text:
using Distributed;
addprocs(20);
@everywhere using PowerModels, Gurobi, Ipopt, Random, Distributions;
@everywhere using CSV;
@everywhere include("solver_run_ac.jl");
@everywhere const GUROBI_ENV = Gurobi.Env();
@everywhere using DataFrames;

list_file = readdir("./data/");
perturbed_file = readdir("./data/perturb/");

perturbResult_dict = Dict();
run_time_dict = Dict();
inf_case = Dict() #the infeasible perturb case

list_file_test = [
    # "case5.m"
    # "case300.m"
    # "case118.m"
    # "case1354pegase.m"
    # "pglib_opf_case10480_goc_api.m"
    "case_ACTIVSg500.m"
]

println("Number of processes: ", nprocs())
println("Number of workers: ", nworkers())

function exist_perturb_file(perturbed_file, itema)
    perturb_data_file = itema[1:end-2] * "_" * "perturb" * ".csv"
    if perturb_data_file in perturbed_file
        return true
    end
    return false
end

for itema in list_file_test
    if occursin(".m", itema)
        if (exist_perturb_file(perturbed_file, itema)) == false
            for σ in [0.05, 0.1, 0.2, 0.3, 0.5]
                time_lim = 10.0
                case_name = AbstractString[]
                obj_val = Float64[]
                time_limit = Float64[]
                solve_time = Float64[]
                Optimal = []
                dual_feasibility = []
                primal_feasibility = []
                esti_obj_val = []
                dual_infeasibility = []
                constraint_violation = []
                overall_nlp = []
                complementarity = []
                var_bound_violation = []
                esti_Optimal = []

                network_data = parse_file("./data/" * itema)
                load_keys = [i for i in keys(network_data["load"])]
                branch_keys = [i for i in keys(network_data["branch"])]
                gen_keys = [i for i in keys(network_data["gen"])]
                bus_keys = [i for i in keys(network_data["bus"])]
                σ_load = Dict()
                load_distrn = Dict()

                for i in 1:length(load_keys)
                    if network_data["load"][load_keys[i]]["pd"] != 0
                        σ_load[load_keys[i]] = abs(network_data["load"][load_keys[i]]["pd"] * σ)
                        load_distrn[load_keys[i]] = Normal(0, σ_load[load_keys[i]])
                    else
                        σ_load[load_keys[i]] = 0
                        load_distrn[load_keys[i]] = Uniform(0, 0.0000001)
                    end
                end

                va_dict = Dict()
                vm_dict = Dict()
                for bus_key in keys(network_data["bus"])
                    va_key = "va" * "_" * bus_key
                    va_value = []

                    vm_key = "vm" * "_" * bus_key
                    vm_value = []

                    va_dict[va_key] = va_value
                    vm_dict[vm_key] = vm_value
                end

                qg_dict = Dict()
                pg_dict = Dict()
                for gen_key in keys(network_data["gen"])
                    qg_key = "qg" * "_" * gen_key
                    qg_value = []

                    pg_key = "pg" * "_" * gen_key
                    pg_value = []

                    qg_dict[qg_key] = qg_value
                    pg_dict[pg_key] = pg_value
                end

                qf_dict = Dict()
                qt_dict = Dict()
                pf_dict = Dict()
                pt_dict = Dict()
                for branch_key in keys(network_data["branch"])
                    qf_key = "qf" * "_" * branch_key
                    qf_value = []

                    qt_key = "qt" * "_" * branch_key
                    qt_value = []

                    pf_key = "pf" * "_" * branch_key
                    pf_value = []

                    pt_key = "pt" * "_" * branch_key
                    pt_value = []

                    pf_dict[pf_key] = pf_value
                    pt_dict[pt_key] = pt_value
                    qt_dict[qt_key] = qt_value
                    qf_dict[qf_key] = qf_value

                end

                pd_dict = Dict()
                qd_dict = Dict()
                for load_key in keys(network_data["load"])
                    pd_key = "pd" * "_" * load_key
                    pd_value = []

                    qd_key = "qd" * "_" * load_key
                    qd_value = []

                    pd_dict[pd_key] = pd_value
                    qd_dict[qd_key] = qd_value

                end

                # perturbResult_dict[itema] = pmap(ω -> run_Ipopt_para_ac(network_data, load_distrn, load_keys, gen_keys, va_dict, vm_dict, qg_dict, pg_dict, qt_dict, qf_dict, pt_dict, pf_dict, pd_dict, qd_dict), 1:1000)

                perturbResult_dict[itema] = pmap(ω -> run_Ipopt_para_ac(itema, network_data, load_distrn, load_keys, gen_keys, time_lim), 1:1000)
                round = 0
                for profile in perturbResult_dict[itema]

                    round = round + 1

                    pd_dict_profile = Dict()

                    case_name_r = itema * "_" * string(round)
                    push!(case_name, case_name_r)
                    push!(obj_val, profile[1])
                    push!(Optimal, profile[2])
                    push!(solve_time, profile[3])
                    push!(primal_feasibility, profile[4])
                    push!(dual_feasibility, profile[5])
                    push!(time_limit, profile[6])
                    push!(esti_obj_val, profile[11])
                    push!(dual_infeasibility, profile[12])
                    push!(constraint_violation, profile[13])
                    push!(overall_nlp, profile[14])
                    push!(complementarity, profile[15])
                    push!(var_bound_violation, profile[16])
                    push!(esti_Optimal, profile[17])

                    for branch_key in [i for i in keys(profile[10])]
                        qf_key = "qf" * "_" * branch_key
                        push!(qf_dict[qf_key], string(profile[10][branch_key]["qf"]))

                        qt_key = "qt" * "_" * branch_key
                        push!(qt_dict[qt_key], string(profile[10][branch_key]["qt"]))

                        pf_key = "pf" * "_" * branch_key
                        push!(pf_dict[pf_key], string(profile[10][branch_key]["pf"]))

                        pt_key = "pt" * "_" * branch_key
                        push!(pt_dict[pt_key], string(profile[10][branch_key]["pt"]))

                    end

                    for gen_key in [i for i in keys(profile[9])]
                        qg_key = "qg" * "_" * gen_key
                        push!(qg_dict[qg_key], string(profile[9][gen_key]["qg"]))

                        pg_key = "pg" * "_" * gen_key
                        push!(pg_dict[pg_key], string(profile[9][gen_key]["pg"]))

                    end

                    for bus_key in [i for i in keys(profile[8])]

                        va_key = "va" * "_" * bus_key
                        push!(va_dict[va_key], string(profile[8][bus_key]["va"]))

                        vm_key = "vm" * "_" * bus_key
                        push!(vm_dict[vm_key], string(profile[8][bus_key]["vm"]))
                    end

                    for load_key in [i for i in keys(profile[7])]
                        pd_key = "pd" * "_" * load_key
                        qd_key = "qd" * "_" * load_key
                        push!(pd_dict[pd_key], string(profile[7][load_key]["pd"]))
                        push!(qd_dict[qd_key], string(profile[7][load_key]["qd"]))

                    end

                end

                result = DataFrame()
                result[!, "case_name"] = case_name
                result[!, "obj_val"] = obj_val
                result[!, "solve_time"] = solve_time
                result[!, "Optimal"] = Optimal
                result[!, "primal_feasibility"] = primal_feasibility
                result[!, "dual_feasibility"] = dual_feasibility
                result[!, "obj_val"] = obj_val
                result[!, "time_limit"] = time_limit
                result[!, "esti_obj_val"] = esti_obj_val
                result[!, "dual_infeasibility"] = dual_infeasibility
                result[!, "constraint_violation"] = constraint_violation
                result[!, "overall_nlp"] = overall_nlp
                result[!, "complementarity"] = complementarity
                result[!, "var_bound_violation"] = var_bound_violation
                result[!, "esti_Optimal"] = esti_obj_val

                for va_key in keys(va_dict)
                    result[!, va_key] = va_dict[va_key]
                end

                for vm_key in keys(vm_dict)
                    result[!, vm_key] = vm_dict[vm_key]
                end

                for pg_key in keys(pg_dict)
                    if length(pg_dict[pg_key]) != 0
                        result[!, pg_key] = pg_dict[pg_key]
                    end
                end

                for qg_key in keys(qg_dict)
                    if length(qg_dict[qg_key]) != 0
                        result[!, qg_key] = qg_dict[qg_key]
                    end
                end

                for pf_key in keys(pf_dict)
                    result[!, pf_key] = pf_dict[pf_key]
                end

                for pt_key in keys(pt_dict)
                    result[!, pt_key] = pt_dict[pt_key]
                end

                for qf_key in keys(qf_dict)
                    result[!, qf_key] = qf_dict[qf_key]
                end

                for qt_key in keys(qt_dict)
                    result[!, qt_key] = qt_dict[qt_key]
                end

                for pd_key in keys(pd_dict)
                    result[!, pd_key] = pd_dict[pd_key]
                end

                for qd_key in keys(qd_dict)
                    result[!, qd_key] = qd_dict[qd_key]
                end

                CSV.write("./data/perturb/" * itema[1:end-2] * "_" * "perturb" * "_" * string(time_lim) * "_" * string(σ) * ".csv", result, bufsize = 2^31)
            end
        end
    end
end

max_stars_repo_name: zhouliang-yu/Grid-ML

max_issues_repo_name: zhouliang-yu/Grid-ML

max_forks_repo_name: zhouliang-yu/Grid-ML


text:
@system StomataBase(Weather, Diffusion) begin
    gs: stomatal_conductance ~ hold
    gb: boundary_layer_conductance ~ hold
    A_net: net_photosynthesis ~ hold
    T: leaf_temperature ~ hold

    drb(Dw, Dc): diffusivity_ratio_boundary_layer => (Dw / Dc)^(2/3) ~ preserve(#= u"H2O/CO2", =# parameter)
    dra(Dw, Dc): diffusivity_ratio_air => (Dw / Dc) ~ preserve(#= u"H2O/CO2", =# parameter)

    Ca(CO2, P_air): co2_air => (CO2 * P_air) ~ track(u"μbar")
    Cs(Ca, A_net, gbc): co2_at_leaf_surface => begin
        Ca - A_net / gbc
    end ~ track(u"μbar")

    gv(gs, gb): total_conductance_h2o => (gs * gb / (gs + gb)) ~ track(u"mol/m^2/s/bar" #= H2O =#)

    rbc(gb, drb): boundary_layer_resistance_co2 => (drb / gb) ~ track(u"m^2*s/mol*bar")
    rsc(gs, dra): stomatal_resistance_co2 => (dra / gs) ~ track(u"m^2*s/mol*bar")
    rvc(rbc, rsc): total_resistance_co2 => (rbc + rsc) ~ track(u"m^2*s/mol*bar")

    gbc(rbc): boundary_layer_conductance_co2 => (1 / rbc) ~ track(u"mol/m^2/s/bar")
    gsc(rsc): stomatal_conductance_co2 => (1 / rsc) ~ track(u"mol/m^2/s/bar")
    gvc(rvc): total_conductance_co2 => (1 / rvc) ~ track(u"mol/m^2/s/bar")
end

@system StomataTuzet begin
    WP_leaf: leaf_water_potential => 0 ~ preserve(u"MPa", parameter)
    Ψv(WP_leaf): bulk_leaf_water_potential ~ track(u"MPa")
    Ψf: reference_leaf_water_potential => -2.0 ~ preserve(u"MPa", parameter)
    sf: stomata_sensitivity_param => 2.3 ~ preserve(u"MPa^-1", parameter)
    fΨv(Ψv, Ψf, sf): stomata_sensitivty => begin
        (1 + exp(sf*Ψf)) / (1 + exp(sf*(Ψf-Ψv)))
    end ~ track
end

@system StomataBallBerry(StomataBase, StomataTuzet) begin
    # Ball-Berry model parameters from Miner and Bauerle 2017, used to be 0.04 and 4.0, respectively (2018-09-04: KDY)
    g0 => 0.017 ~ preserve(u"mol/m^2/s/bar" #= H2O =#, parameter)
    g1 => 4.53 ~ preserve(parameter)

    #HACK: avoid scaling issue with dimensionless unit
    hs(g0, g1, gb, A_net, Cs, fΨv, RH): relative_humidity_at_leaf_surface => begin
        gs = g0 + g1*(A_net*hs/Cs) * fΨv
        (hs - RH)*gb ⩵ (1 - hs)*gs
    end ~ solve(lower=0, upper=1) #, u"percent")
    Ds(D=vp.D, T, hs): vapor_pressure_deficit_at_leaf_surface => begin
        D(T, hs)
    end ~ track(u"kPa")

    gs(g0, g1, A_net, hs, Cs, fΨv): stomatal_conductance => begin
        g0 + g1*(A_net*hs/Cs) * fΨv
    end ~ track(u"mol/m^2/s/bar" #= H2O =#, min=g0)
end

@system StomataMedlyn(StomataBase, StomataTuzet) begin
    g0 => 0.02 ~ preserve(u"mol/m^2/s/bar" #= H2O =#, parameter)
    g1 => 4.0 ~ preserve(u"√kPa", parameter)

    wa(ea=vp.ea, T_air, RH): vapor_pressure_at_air => ea(T_air, RH) ~ track(u"kPa")
    wi(es=vp.es, T): vapor_pressure_at_intercellular_space => es(T) ~ track(u"kPa")
    ws(Ds, wi): vapor_pressure_at_leaf_surface => (wi - Ds) ~ track(u"kPa")
    Ds¹ᐟ²(g0, g1, gb, A_net, Cs, fΨv, wi, wa) => begin
        #HACK: SymPy couldn't extract polynomial coeffs for ps inside √
        gs = g0 + (1 + g1 / Ds¹ᐟ²) * (A_net / Cs) * fΨv
        ws = wi - Ds¹ᐟ²^2
        (ws - wa)*gb ⩵ (wi - ws)*gs
    end ~ solve(lower=0, upper=√wi', u"√kPa")
    Ds(Ds¹ᐟ²): vapor_pressure_deficit_at_leaf_surface => Ds¹ᐟ²^2 ~ track(u"kPa", min=1u"Pa")
    hs(RH=vp.RH, T, Ds): relative_humidity_at_leaf_surface => RH(T, Ds) ~ track

    gs(g0, g1, A_net, Ds, Cs, fΨv): stomatal_conductance => begin
        g0 + (1 + g1/√Ds)*(A_net/Cs) * fΨv
    end ~ track(u"mol/m^2/s/bar" #= H2O =#, min=g0)
end

ext: jl


text:
using Revise

using DrWatson
using GenerativeAD
using StatsBase: fit!, predict, sample
using Statistics
using BSON
using Flux
using DataFrames, CSV
using Random
using EvalMetrics: ConfusionMatrix, recall, precision, f1_score, roccurve, auc_trapezoidal

for dataset in ["ad_cat", "apascal", "bank-additional-full_normalised", "lung-1vs5", "probe", "secom", "u2r"]
    for subsample in 2 .^(3:8)
        parameters = (
            confidence_margin  = 1000.0,
            batchsize          = 256,
            zdim               = 20,
            nlayers            = 1,
            hdim               = 10, # does not matter with one layer
            activation         = "relu",
            ensemble_size      = 50,
            subsample_size     = 16,
            init_seed          = 42
        )

        parameters = merge(parameters, (;subsample_size=subsample))

        # loading data from RDP paper https://github.com/billhhh/RDP/tree/master/RDP-Anomaly/data
        file = datadir("rdp_datasets/$(dataset).csv")
        df = CSV.read(file);

        labels = Int.(df[:class]);
        features = Float32.(copy(Matrix(select(df, Not(:class)))'));

        data = ((features, labels),)
        idim = size(data[1][1], 1)

        @info "Started training REPEN$(parameters) on $(dataset)"
        @info "Number of data points : $(size(data[1][1], 2))"
        @info "Number of features: $(size(data[1][1], 1))"

        model = GenerativeAD.Models.REPEN(;idim=idim, parameters...)
        try
            global info, fit_t, _, _, _ = @timed fit!(model, data; parameters...)
        catch e
            @warn "Training failed due to $e. Continuing."
            continue
        end

        lesinn_trn_scores = GenerativeAD.Models.lesinn(data[1][1]; fit_data=data[1][1], parameters...);
        trn_scores = predict(model, data[1][1]; fit_data=data[1][1], parameters...);

        auc = roccurve(data[1][2], trn_scores) |> auc_trapezoidal
        lesinn_auc = roccurve(data[1][2], lesinn_trn_scores) |> auc_trapezoidal

        @info("", dataset, subsample, auc, lesinn_auc)

        results = Dict(
            :dataset        => dataset,
            :modelname      => "REPEN",
            :history        => info.history,
            :model          => info.model,
            :auc            => auc,
            :lesinn_auc     => lesinn_auc,
            :trn_scores     => trn_scores,
            :trn_labels     => data[1][2],
            :parameters     => parameters,
            :subsample_size => parameters.subsample_size
        )

        wsave(datadir("dumpster/repen/$(dataset)_$(subsample).bson"), results)
    end
end

### analysis
using DataFrames, ValueHistories, CSV.SentinelArrays
files = readdir(datadir("dumpster/repen"), join=true)
results = [load(f) for f in files]

df = DataFrame(
    "dataset" => [r[:dataset] for r in results],
    "subsample" => [r[:subsample_size] for r in results],
    "auc" => [r[:auc] for r in results],
    "lesinn_auc" => [r[:lesinn_auc] for r in results]
)

sort(df, (:dataset, :auc))


text:
# The Computer Language Benchmarks Game
# https://salsa.debian.org/benchmarksgame-team/benchmarksgame/
#
# contributed by David Campbell
# modified by Jarret Revels, Kristoffer Carlsson, Alex Arslan

const revcompdata = Dict(
   'A'=> 'T', 'a'=> 'T',
   'C'=> 'G', 'c'=> 'G',
   'G'=> 'C', 'g'=> 'C',
   'T'=> 'A', 't'=> 'A',
   'U'=> 'A', 'u'=> 'A',
   'M'=> 'K', 'm'=> 'K',
   'R'=> 'Y', 'r'=> 'Y',
   'W'=> 'W', 'w'=> 'W',
   'S'=> 'S', 's'=> 'S',
   'Y'=> 'R', 'y'=> 'R',
   'K'=> 'M', 'k'=> 'M',
   'V'=> 'B', 'v'=> 'B',
   'H'=> 'D', 'h'=> 'D',
   'D'=> 'H', 'd'=> 'H',
   'B'=> 'V', 'b'=> 'V',
   'N'=> 'N', 'n'=> 'N',
)

function print_buff(b)
    br = reverse(b)
    l = length(br)
    for i = 1:60:l
        if i+59 > l
            write(br[i:end]); println()
        else
            write(br[i:i+59]); println()
        end
    end
end

function perf_revcomp()
    buff = UInt8[]
    while true
        line = codeunits(readline())
        if isempty(line)
            print_buff(buff)
            return
        elseif line[1] == UInt8('>')
            print_buff(buff)
            buff = UInt8[]
            write(line)
        else
            l = length(line)-1
            let line = line # Workaround for julia#15276
                append!(buff, [UInt8(revcompdata[Char(line[i])]) for i=1:l])
            end
        end
    end
end

perf_revcomp()


text:
# Sample postprocessing script for virtual satellite trackings.
# Simultaneous virtual satellite trackings at multiple locations.
# Outputs are stored in binary format for sharing within Julia.
#
# Usage:
#   julia -t nthreads demo_pointsextracting.jl
# or
#   JULIA_NUM_THREADS=nthreads julia demo_pointsextracting.jl
#
# Hongyang Zhou, hyzhou@umich.edu

using Glob, Vlasiator
using Vlasiator: RE # Earth radius, [m]
using JLD2: jldsave

function extract_vars(files, locations)
   nfiles = length(files)
   nlocs = length(locations)
   # variables to be extracted
   t   = zeros(Float32, nfiles)
   rho = zeros(Float32, nlocs, nfiles)
   v   = zeros(Float32, 3, nlocs, nfiles)
   p   = zeros(Float32, nlocs, nfiles)
   b   = zeros(Float32, 3, nlocs, nfiles)
   e   = zeros(Float32, 3, nlocs, nfiles)

   ids = Vector{UInt64}(undef, nlocs)
   load(files[1]) do meta
      for iloc in eachindex(locations)
         ids[iloc] = getcell(meta, locations[iloc])
      end
   end

   # Extract data from each frame
   Threads.@threads for i = eachindex(files)
      meta = load(files[i])
      t[i] = meta.time
      rho[:,i] = readvariable(meta, "proton/vg_rho", ids)
      v[:,:,i] = readvariable(meta, "proton/vg_v",   ids)
      p[  :,i] = readvariable(meta, "vg_pressure",   ids)
      b[:,:,i] = readvariable(meta, "vg_b_vol",      ids)
      e[:,:,i] = readvariable(meta, "vg_e_vol",      ids)
   end

   # Save into binary file
   jldsave("satellites.jld2"; locations, t, rho, v, p, b, e)
end

#####

files = glob("bulk*.vlsv", "./")

# virtual satellite locations
locations = [[7RE, 0, 0], [9RE, 0, 0], [11RE, 0, 0], [12RE, 0, 0], [13RE, 0, 0],
   [14RE, 0, 0], [15RE, 0, 0], [16RE, 0, 0], [17RE, 0, 0], [29.3RE, 0, 0]]

println("Number of files: $(length(files))")
println("Number of virtual satellites: $(length(locations))")
println("Running with $(Threads.nthreads()) threads...")

@time extract_vars(files, locations)

println("Virtual satellite extraction done!")


text:
module Diff
using LinearAlgebra
using OffsetArrays
using PaddedViews
export diff_grids
using ..Render, ..Grids

index_overlap(A,B) = length(intersect(CartesianIndices(A),CartesianIndices(B)))

function diff_grids(A::ARCGrid, B::ARCGrid)
    Asz1,Asz2 = size(A)
    Bsz1,Bsz2 = size(B)
    """
    Slide B over A
    Initially B sits entirely inside of A (or vis versa if sizes are swapped, but that doesn't affect anything)
    with [1,1] as their shared top left corner.
    We will be using OffsetArrays here, and specifically B will be the offset array
    and A will always keep its indices (1:Asz1, 1:Asz2).
    We want to start B such that the bottom right pixel of B lines up with the top left pixel of A.
    Since we use the indices of A this is [1,1] wrt A, and so to make the bottom right [Bsz1,Bsz2] pixel
    of B point to [1,1] we basically want all of Bs indices to be <=0 until that one pixel. This ends up
    being `OffsetArray(B, -Bsz1+1, -Bsz1+1)`.
    """
    # todo: precompute these instead of computing them every time
    edges_A = edges_of_grid(A)
    edges_B = edges_of_grid(B)

    oBinit = OffsetArray(B, -Bsz1+1, -Bsz2+1)
    edges_oBinit = OffsetArray(edges_B, -Bsz1+1, -Bsz2+1, 0)

    # highest offests to go to
    d1_end = Asz1+Bsz1-2
    d2_end = Asz2+Bsz2-2

    # sanity checks
    @assert index_overlap(A,oBinit) == 1
    @assert index_overlap(A,OffsetArray(oBinit, d1_end, d2_end)) == 1

    diff_grid = ARCDiffGrid(
        Matrix{ARCDiff}(undef, d1_end+1, d2_end+1),
        A::ARCGrid,
        B::ARCGrid,
    )

    for d1 in 0:d1_end, d2 in 0:d2_end
        # slide B using offset
        oB = OffsetArray(oBinit, d1, d2)
        edges_oB = OffsetArray(edges_oBinit, d1, d2, 0)

        # restrict to the overlapping region
        # shared_indices = intersect(CartesianIndices(A),CartesianIndices(oB))
        # oBint = oB[shared_indices]
        # Aint = A[shared_indices]
        # edges_Aint = edges_A[shared_indices,:]
        # edges_oBint = edges_oB[shared_indices,:]

        # all of the things paddedviews will return will be offsetarrays! so the same elems
        # are still in all the same locations they were in A and oB before.
        Aint, oBint = paddedviews(-1, A, oB)
        edges_Aint, edges_oBint = paddedviews(false,edges_A, edges_oB)

        # nonzero pixels that match. Match will be indexed just like A.
        match = falses(size(A))
        match[findall(Aint .== oBint .> 0)] .= true
        # @show typeof(match)

        # edges that match. These will be indexed just like A.
        edges_match = falses(size(edges_A))
        edges_match[findall(edges_Aint .== edges_oBint .== true)] .= true
        # edges_match,_ = paddedview(0, edges_match,A) # pad edges_match to match A

        diff_grid.grid[d1+1,d2+1] = ARCDiff(A,oB,match,edges_match)

        # helpful debug comment:
        # println("$d1 $d2 -> $(sum(match))/$(index_overlap(A,oB))")
    end
    diff_grid
end

end


text:
using Requires
@init @require CUDAnative = "be33ccc6-a3ff-5ff2-a52e-74243cff1e17" begin
  using .CUDAnative
end

# {{{ FIXME: remove this after we've figure out how to pass through to kernel
const _ξ1x1, _ξ2x1, _ξ3x1 = Grids._ξ1x1, Grids._ξ2x1, Grids._ξ3x1
const _ξ1x2, _ξ2x2, _ξ3x2 = Grids._ξ1x2, Grids._ξ2x2, Grids._ξ3x2
const _ξ1x3, _ξ2x3, _ξ3x3 = Grids._ξ1x3, Grids._ξ2x3, Grids._ξ3x3
const _M, _MI = Grids._M, Grids._MI
const _x1, _x2, _x3 = Grids._x1, Grids._x2, Grids._x3
const _JcV = Grids._JcV

const _n1, _n2, _n3 = Grids._n1, Grids._n2, Grids._n3
const _sM, _vMI = Grids._sM, Grids._vMI
# }}}

"""
    volumerhs!(::Val{dim}, ::Val{N}, ::Val{nstate}, ::Val{nviscstate},
               ::Val{nauxstate}, flux!, source!, rhs, Q, Qvisc, auxstate,
               vgeo, t, D, elems) where {dim, N, nstate, nviscstate,

Computational kernel: Evaluate the volume integrals on right-hand side of a
`DGBalanceLaw` semi-discretization.

See [`odefun!`](@ref) for usage.
"""
function volumerhs!(::Val{dim}, ::Val{N},
                    ::Val{nstate}, ::Val{nviscstate},
                    ::Val{nauxstate},
                    flux!, source!,
                    rhs, Q, Qvisc, auxstate, vgeo, t,
                    ω, D, elems, increment) where {dim, N, nstate, nviscstate,
                                     nauxstate}
  FT = eltype(Q)

  Nq = N + 1

  Nqk = dim == 2 ? 1 : Nq

  s_F = @shmem FT (3, Nq, Nq, Nqk, nstate)
  s_ω = @shmem FT (Nq, )
  s_half_D = @shmem FT (Nq, Nq)
  l_rhs = @scratch FT (nstate, Nq, Nq, Nqk) 3

  source! !== nothing && (l_S = MArray{Tuple{nstate}, FT}(undef))
  l_Q = MArray{Tuple{nstate}, FT}(undef)
  l_Qvisc = MArray{Tuple{nviscstate}, FT}(undef)
  l_aux = MArray{Tuple{nauxstate}, FT}(undef)
  l_F = MArray{Tuple{3, nstate}, FT}(undef)
  l_M = @scratch FT (Nq, Nq, Nqk) 3
  l_ξ1x1 = @scratch FT (Nq, Nq, Nqk) 3
  l_ξ1x2 = @scratch FT (Nq, Nq, Nqk) 3
  l_ξ1x3 = @scratch FT (Nq, Nq, Nqk) 3
  l_ξ2x1 = @scratch FT (Nq, Nq, Nqk) 3
  l_ξ2x2 = @scratch FT (Nq, Nq, Nqk) 3
  l_ξ2x3 = @scratch FT (Nq, Nq, Nqk) 3
  l_ξ3x1 = @scratch FT (Nq, Nq, Nqk) 3
  l_ξ3x2 = @scratch FT (Nq, Nq, Nqk) 3
  l_ξ3x3 = @scratch FT (Nq, Nq, Nqk) 3

  @inbounds @loop for k in (1; threadIdx().z)
    @loop for j in (1:Nq; threadIdx().y)
      s_ω[j] = ω[j]
      @loop for i in (1:Nq; threadIdx().x)
        s_half_D[i, j] = D[i, j] / 2
      end
    end
  end

  @inbounds @loop for e in (elems; blockIdx().x)
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          ijk = i + Nq * ((j-1) + Nq * (k-1))
          l_M[i, j, k] = vgeo[ijk, _M, e]
          l_ξ1x1[i, j, k] = vgeo[ijk, _ξ1x1, e]
          l_ξ1x2[i, j, k] = vgeo[ijk, _ξ1x2, e]
          l_ξ1x3[i, j, k] = vgeo[ijk, _ξ1x3, e]
          l_ξ2x1[i, j, k] = vgeo[ijk, _ξ2x1, e]
          l_ξ2x2[i, j, k] = vgeo[ijk, _ξ2x2, e]
          l_ξ2x3[i, j, k] = vgeo[ijk, _ξ2x3, e]
          l_ξ3x1[i, j, k] = vgeo[ijk, _ξ3x1, e]
          l_ξ3x2[i, j, k] = vgeo[ijk, _ξ3x2, e]
          l_ξ3x3[i, j, k] = vgeo[ijk, _ξ3x3, e]

          @unroll for s = 1:nstate
            l_rhs[s, i, j, k] = increment ? rhs[ijk, s, e] : zero(FT)
          end

          @unroll for s = 1:nstate
            l_Q[s] = Q[ijk, s, e]
          end

          @unroll for s = 1:nviscstate
            l_Qvisc[s] = Qvisc[ijk, s, e]
          end

          @unroll for s = 1:nauxstate
            l_aux[s] = auxstate[ijk, s, e]
          end

          flux!(l_F, l_Q, l_Qvisc, l_aux, t)

          @unroll for s = 1:nstate
            s_F[1,i,j,k,s] = l_F[1,s]
            s_F[2,i,j,k,s] = l_F[2,s]
            s_F[3,i,j,k,s] = l_F[3,s]
          end

          if source! !== nothing
            source!(l_S, l_Q, l_Qvisc, l_aux, t)

            @unroll for s = 1:nstate
              l_rhs[s, i, j, k] += l_S[s]
            end
          end
        end
      end
    end
    @synchronize

    # Weak "outside metrics" derivative
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          @unroll for n = 1:Nq
            @unroll for s = 1:nstate
              Dni = s_half_D[n, i] * s_ω[n] / s_ω[i]
              Dnj = s_half_D[n, j] * s_ω[n] / s_ω[j]
              Nqk > 1 && (Dnk = s_half_D[n, k] * s_ω[n] / s_ω[k])

              # ξ1-grid lines
              l_rhs[s, i, j, k] += l_ξ1x1[i, j, k] * Dni * s_F[1, n, j, k, s]
              l_rhs[s, i, j, k] += l_ξ1x2[i, j, k] * Dni * s_F[2, n, j, k, s]
              l_rhs[s, i, j, k] += l_ξ1x3[i, j, k] * Dni * s_F[3, n, j, k, s]

              # ξ2-grid lines
              l_rhs[s, i, j, k] += l_ξ2x1[i, j, k] * Dnj * s_F[1, i, n, k, s]
              l_rhs[s, i, j, k] += l_ξ2x2[i, j, k] * Dnj * s_F[2, i, n, k, s]
              l_rhs[s, i, j, k] += l_ξ2x3[i, j, k] * Dnj * s_F[3, i, n, k, s]

              # ξ3-grid lines
              if Nqk > 1
                l_rhs[s, i, j, k] += l_ξ3x1[i, j, k] * Dnk * s_F[1, i, j, n, s]
                l_rhs[s, i, j, k] += l_ξ3x2[i, j, k] * Dnk * s_F[2, i, j, n, s]
                l_rhs[s, i, j, k] += l_ξ3x3[i, j, k] * Dnk * s_F[3, i, j, n, s]
              end
            end
          end
        end
      end
    end
    @synchronize

    # Build "inside metrics" flux
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          @unroll for s = 1:nstate
            F1, F2, F3 = s_F[1,i,j,k,s], s_F[2,i,j,k,s], s_F[3,i,j,k,s]
            s_F[1,i,j,k,s] = l_M[i, j, k] * (l_ξ1x1[i, j, k] * F1 +
                                              l_ξ1x2[i, j, k] * F2 +
                                              l_ξ1x3[i, j, k] * F3)
            s_F[2,i,j,k,s] = l_M[i, j, k] * (l_ξ2x1[i, j, k] * F1 +
                                              l_ξ2x2[i, j, k] * F2 +
                                              l_ξ2x3[i, j, k] * F3)
            s_F[3,i,j,k,s] = l_M[i, j, k] * (l_ξ3x1[i, j, k] * F1 +
                                              l_ξ3x2[i, j, k] * F2 +
                                              l_ξ3x3[i, j, k] * F3)
          end
        end
      end
    end
    @synchronize

    # Weak "inside metrics" derivative
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          ijk = i + Nq * ((j-1) + Nq * (k-1))
          MI = vgeo[ijk, _MI, e]
          @unroll for s = 1:nstate
            @unroll for n = 1:Nq
              Dni = s_half_D[n, i]
              Dnj = s_half_D[n, j]
              Nqk > 1 && (Dnk = s_half_D[n, k])
              # ξ1-grid lines
              l_rhs[s, i, j, k] += MI * Dni * s_F[1, n, j, k, s]

              # ξ2-grid lines
              l_rhs[s, i, j, k] += MI * Dnj * s_F[2, i, n, k, s]

              # ξ3-grid lines
              Nqk > 1 && (l_rhs[s, i, j, k] += MI * Dnk * s_F[3, i, j, n, s])
            end
          end
        end
      end
    end
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          ijk = i + Nq * ((j-1) + Nq * (k-1))
          @unroll for s = 1:nstate
            rhs[ijk, s, e] = l_rhs[s, i, j, k]
          end
        end
      end
    end
    @synchronize
  end
  nothing
end

"""
    facerhs!(::Val{dim}, ::Val{N}, ::Val{nstate}, ::Val{nviscstate},
             ::Val{nauxstate}, numerical_flux!,
             numerical_boundary_flux!, rhs, Q, Qvisc, auxstate,
             vgeo, sgeo, t, vmapM, vmapP, elemtobndy,
             elems) where {dim, N, nstate, nviscstate, nauxstate}

Computational kernel: Evaluate the surface integrals on right-hand side of a
`DGBalanceLaw` semi-discretization.

See [`odefun!`](@ref) for usage.
"""
function facerhs!(::Val{dim}, ::Val{N}, ::Val{nstate}, ::Val{nviscstate},
                  ::Val{nauxstate}, numerical_flux!, numerical_boundary_flux!,
                  rhs, Q, Qvisc, auxstate, vgeo, sgeo, t, vmapM, vmapP,
                  elemtobndy, elems) where {dim, N, nstate, nviscstate,
                                            nauxstate}

  FT = eltype(Q)

  if dim == 1
    Np = (N+1)
    Nfp = 1
    nface = 2
  elseif dim == 2
    Np = (N+1) * (N+1)
    Nfp = (N+1)
    nface = 4
  elseif dim == 3
    Np = (N+1) * (N+1) * (N+1)
    Nfp = (N+1) * (N+1)
    nface = 6
  end

  l_QM = MArray{Tuple{nstate}, FT}(undef)
  l_QviscM = MArray{Tuple{nviscstate}, FT}(undef)
  l_auxM = MArray{Tuple{nauxstate}, FT}(undef)

  l_QP = MArray{Tuple{nstate}, FT}(undef)
  l_QviscP = MArray{Tuple{nviscstate}, FT}(undef)
  l_auxP = MArray{Tuple{nauxstate}, FT}(undef)

  l_F = MArray{Tuple{nstate}, FT}(undef)

  @inbounds @loop for e in (elems; blockIdx().x)
    for f = 1:nface
      @loop for n in (1:Nfp; threadIdx().x)
        nM = (sgeo[_n1, n, f, e], sgeo[_n2, n, f, e], sgeo[_n3, n, f, e])
        sM, vMI = sgeo[_sM, n, f, e], sgeo[_vMI, n, f, e]
        idM, idP = vmapM[n, f, e], vmapP[n, f, e]

        eM, eP = e, ((idP - 1) ÷ Np) + 1
        vidM, vidP = ((idM - 1) % Np) + 1,  ((idP - 1) % Np) + 1

        # Load minus side data
        @unroll for s = 1:nstate
          l_QM[s] = Q[vidM, s, eM]
        end

        @unroll for s = 1:nviscstate
          l_QviscM[s] = Qvisc[vidM, s, eM]
        end

        @unroll for s = 1:nauxstate
          l_auxM[s] = auxstate[vidM, s, eM]
        end

        # Load plus side data
        @unroll for s = 1:nstate
          l_QP[s] = Q[vidP, s, eP]
        end

        @unroll for s = 1:nviscstate
          l_QviscP[s] = Qvisc[vidP, s, eP]
        end

        @unroll for s = 1:nauxstate
          l_auxP[s] = auxstate[vidP, s, eP]
        end

        bctype =
            numerical_boundary_flux! === nothing ? 0 : elemtobndy[f, e]
        if bctype == 0
          numerical_flux!(l_F, nM, l_QM, l_QviscM, l_auxM, l_QP, l_QviscP,
                          l_auxP, t)
        else
          numerical_boundary_flux!(l_F, nM, l_QM, l_QviscM, l_auxM, l_QP,
                                   l_QviscP, l_auxP, bctype, t)
        end

        #Update RHS
        @unroll for s = 1:nstate
          # FIXME: Should we pretch these?
          rhs[vidM, s, eM] -= vMI * sM * l_F[s]
        end
      end
      # Need to wait after even faces to avoid race conditions
      f % 2 == 0 && @synchronize
    end
  end
  nothing
end

function volumeviscterms!(::Val{dim}, ::Val{N}, ::Val{nstate},
                          ::Val{ngradstate}, ::Val{nviscstate},
                          ::Val{nauxstate}, viscous_transform!,
                          gradient_transform!, Q, Qvisc, auxstate, vgeo, t, D,
                          elems) where {dim, N, ngradstate,
                                        nviscstate, nstate, nauxstate}
  FT = eltype(Q)

  Nq = N + 1

  Nqk = dim == 2 ? 1 : Nq

  s_G = @shmem FT (Nq, Nq, Nqk, ngradstate)
  s_D = @shmem FT (Nq, Nq)

  l_Q = @scratch FT (nstate, Nq, Nq, Nqk) 3
  l_aux = @scratch FT (nauxstate, Nq, Nq, Nqk) 3
  l_G = MArray{Tuple{ngradstate}, FT}(undef)
  l_Qvisc = MArray{Tuple{nviscstate}, FT}(undef)
  l_gradG = MArray{Tuple{3, ngradstate}, FT}(undef)

  @inbounds @loop for k in (1; threadIdx().z)
    @loop for j in (1:Nq; threadIdx().y)
      @loop for i in (1:Nq; threadIdx().x)
        s_D[i, j] = D[i, j]
      end
    end
  end

  @inbounds @loop for e in (elems; blockIdx().x)
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          ijk = i + Nq * ((j-1) + Nq * (k-1))
          @unroll for s = 1:nstate
            l_Q[s, i, j, k] = Q[ijk, s, e]
          end

          @unroll for s = 1:nauxstate
            l_aux[s, i, j, k] = auxstate[ijk, s, e]
          end

          gradient_transform!(l_G, l_Q[:, i, j, k], l_aux[:, i, j, k], t)
          @unroll for s = 1:ngradstate
            s_G[i, j, k, s] = l_G[s]
          end
        end
      end
    end
    @synchronize

    # Compute gradient of each state
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          ijk = i + Nq * ((j-1) + Nq * (k-1))
          ξ1x1, ξ1x2, ξ1x3 = vgeo[ijk, _ξ1x1, e], vgeo[ijk, _ξ1x2, e], vgeo[ijk, _ξ1x3, e]
          ξ2x1, ξ2x2, ξ2x3 = vgeo[ijk, _ξ2x1, e], vgeo[ijk, _ξ2x2, e], vgeo[ijk, _ξ2x3, e]
          ξ3x1, ξ3x2, ξ3x3 = vgeo[ijk, _ξ3x1, e], vgeo[ijk, _ξ3x2, e], vgeo[ijk, _ξ3x3, e]

          @unroll for s = 1:ngradstate
            Gξ1 = Gξ2 = Gξ3 = zero(FT)
            @unroll for n = 1:Nq
              Din = s_D[i, n]
              Djn = s_D[j, n]
              Nqk > 1 && (Dkn = s_D[k, n])

              Gξ1 += Din * s_G[n, j, k, s]
              Gξ2 += Djn * s_G[i, n, k, s]
              Nqk > 1 && (Gξ3 += Dkn * s_G[i, j, n, s])
            end
            l_gradG[1, s] = ξ1x1 * Gξ1 + ξ2x1 * Gξ2 + ξ3x1 * Gξ3
            l_gradG[2, s] = ξ1x2 * Gξ1 + ξ2x2 * Gξ2 + ξ3x2 * Gξ3
            l_gradG[3, s] = ξ1x3 * Gξ1 + ξ2x3 * Gξ2 + ξ3x3 * Gξ3
          end

          viscous_transform!(l_Qvisc, l_gradG, l_Q[:, i, j, k],
                             l_aux[:, i, j, k], t)

          @unroll for s = 1:nviscstate
            Qvisc[ijk, s, e] = l_Qvisc[s]
          end
        end
      end
    end
    @synchronize
  end
end

function faceviscterms!(::Val{dim}, ::Val{N}, ::Val{nstate},
                        ::Val{ngradstate}, ::Val{nviscstate},
                        ::Val{nauxstate}, viscous_penalty!,
                        viscous_boundary_penalty!, gradient_transform!,
                        Q, Qvisc, auxstate, vgeo, sgeo, t, vmapM, vmapP,
                        elemtobndy, elems) where {dim, N, ngradstate,
                                                  nviscstate, nstate, nauxstate}
  FT = eltype(Q)

  if dim == 1
    Np = (N+1)
    Nfp = 1
    nface = 2
  elseif dim == 2
    Np = (N+1) * (N+1)
    Nfp = (N+1)
    nface = 4
  elseif dim == 3
    Np = (N+1) * (N+1) * (N+1)
    Nfp = (N+1) * (N+1)
    nface = 6
  end

  l_QM = MArray{Tuple{nstate}, FT}(undef)
  l_auxM = MArray{Tuple{nauxstate}, FT}(undef)
  l_GM = MArray{Tuple{ngradstate}, FT}(undef)

  l_QP = MArray{Tuple{nstate}, FT}(undef)
  l_auxP = MArray{Tuple{nauxstate}, FT}(undef)
  l_GP = MArray{Tuple{ngradstate}, FT}(undef)

  l_Qvisc = MArray{Tuple{nviscstate}, FT}(undef)

  @inbounds @loop for e in (elems; blockIdx().x)
    for f = 1:nface
      @loop for n in (1:Nfp; threadIdx().x)
        nM = (sgeo[_n1, n, f, e], sgeo[_n2, n, f, e], sgeo[_n3, n, f, e])
        sM, vMI = sgeo[_sM, n, f, e], sgeo[_vMI, n, f, e]
        idM, idP = vmapM[n, f, e], vmapP[n, f, e]

        eM, eP = e, ((idP - 1) ÷ Np) + 1
        vidM, vidP = ((idM - 1) % Np) + 1,  ((idP - 1) % Np) + 1

        # Load minus side data
        @unroll for s = 1:nstate
          l_QM[s] = Q[vidM, s, eM]
        end

        @unroll for s = 1:nauxstate
          l_auxM[s] = auxstate[vidM, s, eM]
        end

        gradient_transform!(l_GM, l_QM, l_auxM, t)

        # Load plus side data
        @unroll for s = 1:nstate
          l_QP[s] = Q[vidP, s, eP]
        end

        @unroll for s = 1:nauxstate
          l_auxP[s] = auxstate[vidP, s, eP]
        end

        gradient_transform!(l_GP, l_QP, l_auxP, t)

        bctype =
            viscous_boundary_penalty! === nothing ? 0 : elemtobndy[f, e]
        if bctype == 0
          viscous_penalty!(l_Qvisc, nM, l_GM, l_QM, l_auxM, l_GP,
                                  l_QP, l_auxP, t)
        else
          viscous_boundary_penalty!(l_Qvisc, nM, l_GM, l_QM, l_auxM,
                                           l_GP, l_QP, l_auxP, bctype, t)
        end

        @unroll for s = 1:nviscstate
          Qvisc[vidM, s, eM] += vMI * sM * l_Qvisc[s]
        end
      end
      # Need to wait after even faces to avoid race conditions
      f % 2 == 0 && @synchronize
    end
  end
  nothing
end

"""
    initstate!(::Val{dim}, ::Val{N}, ::Val{nvar}, ::Val{nauxstate},
               ic!, Q, auxstate, vgeo, elems) where {dim, N, nvar, nauxstate}

Computational kernel: Initialize the state

See [`DGBalanceLaw`](@ref) for usage.
"""
function initstate!(::Val{dim}, ::Val{N}, ::Val{nvar}, ::Val{nauxstate},
                    ic!, Q, auxstate, vgeo, elems) where {dim, N, nvar, nauxstate}

  FT = eltype(Q)

  Nq = N + 1
  Nqk = dim == 2 ? 1 : Nq
  Np = Nq * Nq * Nqk

  l_Qdof = MArray{Tuple{nvar}, FT}(undef)
  l_auxdof = MArray{Tuple{nauxstate}, FT}(undef)

  @inbounds @loop for e in (elems; blockIdx().x)
    @loop for i in (1:Np; threadIdx().x)
      x1, x2, x3 = vgeo[i, _x1, e], vgeo[i, _x2, e], vgeo[i, _x3, e]

      @unroll for s = 1:nauxstate
        l_auxdof[s] = auxstate[i, s, e]
      end
      ic!(l_Qdof, x1, x2, x3, l_auxdof)
      @unroll for n = 1:nvar
        Q[i, n, e] = l_Qdof[n]
      end
    end
  end
end

"""
    initauxstate!(::Val{dim}, ::Val{N}, ::Val{nauxstate}, auxstatefun!,
                  auxstate, vgeo, elems) where {dim, N, nauxstate}

Computational kernel: Initialize the auxiliary state

See [`DGBalanceLaw`](@ref) for usage.
"""
function initauxstate!(::Val{dim}, ::Val{N}, ::Val{nauxstate}, auxstatefun!,
                       auxstate, vgeo, elems) where {dim, N, nauxstate}

  FT = eltype(auxstate)

  Nq = N + 1
  Nqk = dim == 2 ? 1 : Nq
  Np = Nq * Nq * Nqk

  l_aux = MArray{Tuple{nauxstate}, FT}(undef)

  @inbounds @loop for e in (elems; blockIdx().x)
    @loop for n in (1:Np; threadIdx().x)
      x1, x2, x3 = vgeo[n, _x1, e], vgeo[n, _x2, e], vgeo[n, _x3, e]
      @unroll for s = 1:nauxstate
        l_aux[s] = auxstate[n, s, e]
      end

      auxstatefun!(l_aux, x1, x2, x3)

      @unroll for s = 1:nauxstate
        auxstate[n, s, e] = l_aux[s]
      end
    end
  end
end

"""
    elem_grad_field!(::Val{dim}, ::Val{N}, ::Val{nstate}, Q, vgeo, D, elems, s,
                     sx, sy, sz) where {dim, N, nstate}

Computational kernel: Compute the element gradient of state `s` of `Q` and store
it in `sx`, `sy`, and `sz` of `Q`.

!!! warning

    This does not compute a DG gradient, but only over the element. If ``Q_s``
    is discontinuous you may want to consider another approach.

"""
function elem_grad_field!(::Val{dim}, ::Val{N}, ::Val{nstate}, Q, vgeo,
                          ω, D, elems, s, sx, sy, sz) where {dim, N, nstate}

  FT = eltype(vgeo)

  Nq = N + 1

  Nqk = dim == 2 ? 1 : Nq

  s_f = @shmem FT (3, Nq, Nq, Nqk)
  s_half_D = @shmem FT (Nq, Nq)

  l_f  = @scratch FT (Nq, Nq, Nqk) 3
  l_fd = @scratch FT (3, Nq, Nq, Nqk) 3

  l_J = @scratch FT (Nq, Nq, Nqk) 3
  l_ξ1d = @scratch FT (3, Nq, Nq, Nqk) 3
  l_ξ2d = @scratch FT (3, Nq, Nq, Nqk) 3
  l_ξ3d = @scratch FT (3, Nq, Nq, Nqk) 3

  @inbounds @loop for k in (1; threadIdx().z)
    @loop for j in (1:Nq; threadIdx().y)
      @loop for i in (1:Nq; threadIdx().x)
        s_half_D[i, j] = D[i, j] / 2
      end
    end
  end

  @inbounds @loop for e in (elems; blockIdx().x)
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          ijk = i + Nq * ((j-1) + Nq * (k-1))
          M = dim == 2 ? ω[i] * ω[j] : ω[i] * ω[j] * ω[k]
          l_J[i, j, k] = vgeo[ijk, _M, e] / M
          l_ξ1d[1, i, j, k] = vgeo[ijk, _ξ1x1, e]
          l_ξ1d[2, i, j, k] = vgeo[ijk, _ξ1x2, e]
          l_ξ1d[3, i, j, k] = vgeo[ijk, _ξ1x3, e]
          l_ξ2d[1, i, j, k] = vgeo[ijk, _ξ2x1, e]
          l_ξ2d[2, i, j, k] = vgeo[ijk, _ξ2x2, e]
          l_ξ2d[3, i, j, k] = vgeo[ijk, _ξ2x3, e]
          l_ξ3d[1, i, j, k] = vgeo[ijk, _ξ3x1, e]
          l_ξ3d[2, i, j, k] = vgeo[ijk, _ξ3x2, e]
          l_ξ3d[3, i, j, k] = vgeo[ijk, _ξ3x3, e]
          l_f[i, j, k] = s_f[1, i, j, k] = Q[ijk, s, e]
        end
      end
    end
    @synchronize

    # reference gradient: outside metrics
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          fξ1 = FT(0)
          fξ2 = FT(0)
          fξ3 = FT(0)
          @unroll for n = 1:Nq
            Din = s_half_D[i, n]
            Djn = s_half_D[j, n]
            Nqk > 1 && (Dkn = s_half_D[k, n])

            # ξ1-grid lines
            fξ1 += Din * s_f[1, n, j, k]

            # ξ2-grid lines
            fξ2 += Djn * s_f[1, i, n, k]

            # ξ3-grid lines
            Nqk > 1 && (fξ3 += Dkn * s_f[1, i, j, n])
          end
          @unroll for d = 1:3
            l_fd[d, i, j, k] = l_ξ1d[d, i, j, k] * fξ1 +
                               l_ξ2d[d, i, j, k] * fξ2 +
                               l_ξ3d[d, i, j, k] * fξ3
          end
        end
      end
    end
    @synchronize

    # Build "inside metrics" flux
    for d = 1:3
      @loop for k in (1:Nqk; threadIdx().z)
        @loop for j in (1:Nq; threadIdx().y)
          @loop for i in (1:Nq; threadIdx().x)
            s_f[1,i,j,k] = l_J[i, j, k] * l_ξ1d[d, i, j, k] * l_f[i, j, k]
            s_f[2,i,j,k] = l_J[i, j, k] * l_ξ2d[d, i, j, k] * l_f[i, j, k]
            s_f[3,i,j,k] = l_J[i, j, k] * l_ξ3d[d, i, j, k] * l_f[i, j, k]
          end
        end
      end
      @synchronize
      @loop for k in (1:Nqk; threadIdx().z)
        @loop for j in (1:Nq; threadIdx().y)
          @loop for i in (1:Nq; threadIdx().x)
            fd = FT(0)
            JI = 1 / l_J[i, j, k]
            @unroll for n = 1:Nq
              Din = s_half_D[i, n]
              Djn = s_half_D[j, n]
              Nqk > 1 && (Dkn = s_half_D[k, n])

              l_fd[d, i, j, k] += JI * Din * s_f[1, n, j, k]
              l_fd[d, i, j, k] += JI * Djn * s_f[2, i, n, k]
              Nqk > 1 && (l_fd[d, i, j, k] += JI * Dkn * s_f[3, i, j, n])
            end
          end
        end
      end
      @synchronize
    end

    # Physical gradient
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          ijk = i + Nq * ((j-1) + Nq * (k-1))
          Q[ijk, sx, e] = l_fd[1, i, j, k]
          Q[ijk, sy, e] = l_fd[2, i, j, k]
          Q[ijk, sz, e] = l_fd[3, i, j, k]
        end
      end
    end
    @synchronize
  end
end

"""
    knl_dof_iteration!(::Val{dim}, ::Val{N}, ::Val{nRstate}, ::Val{nstate},
                       ::Val{nviscstate}, ::Val{nauxstate}, dof_fun!, R, Q,
                       QV, auxstate, elems) where {dim, N, nRstate, nstate,
                                                   nviscstate, nauxstate}

Computational kernel: fill postprocessing array

See [`DGBalanceLaw`](@ref) for usage.
"""
function knl_dof_iteration!(::Val{dim}, ::Val{N}, ::Val{nRstate}, ::Val{nstate},
                            ::Val{nviscstate}, ::Val{nauxstate}, dof_fun!, R, Q,
                            QV, auxstate, elems) where {dim, N, nRstate, nstate,
                                                        nviscstate, nauxstate}
  FT = eltype(R)

  Nq = N + 1

  Nqk = dim == 2 ? 1 : Nq

  Np = Nq * Nq * Nqk

  l_R = MArray{Tuple{nRstate}, FT}(undef)
  l_Q = MArray{Tuple{nstate}, FT}(undef)
  l_Qvisc = MArray{Tuple{nviscstate}, FT}(undef)
  l_aux = MArray{Tuple{nauxstate}, FT}(undef)

  @inbounds @loop for e in (elems; blockIdx().x)
    @loop for n in (1:Np; threadIdx().x)
      @unroll for s = 1:nRstate
        l_R[s] = R[n, s, e]
      end

      @unroll for s = 1:nstate
        l_Q[s] = Q[n, s, e]
      end

      @unroll for s = 1:nviscstate
        l_Qvisc[s] = QV[n, s, e]
      end

      @unroll for s = 1:nauxstate
        l_aux[s] = auxstate[n, s, e]
      end

      dof_fun!(l_R, l_Q, l_Qvisc, l_aux)

      @unroll for s = 1:nRstate
        R[n, s, e] = l_R[s]
      end
    end
  end
end

"""
    knl_indefinite_stack_integral!(::Val{dim}, ::Val{N}, ::Val{nstate},
                                            ::Val{nauxstate}, ::Val{nvertelem},
                                            int_knl!, Q, auxstate, vgeo, Imat,
                                            elems, ::Val{outstate}
                                           ) where {dim, N, nstate, nauxstate,
                                                    outstate, nvertelem}

Computational kernel: compute indefinite integral along the vertical stack

See [`DGBalanceLaw`](@ref) for usage.
"""
function knl_indefinite_stack_integral!(::Val{dim}, ::Val{N}, ::Val{nstate},
                                        ::Val{nauxstate}, ::Val{nvertelem},
                                        int_knl!, P, Q, auxstate, vgeo, Imat,
                                        elems, ::Val{outstate}
                                       ) where {dim, N, nstate, nauxstate,
                                                outstate, nvertelem}
  FT = eltype(Q)

  Nq = N + 1
  Nqj = dim == 2 ? 1 : Nq

  nout = length(outstate)

  l_Q = MArray{Tuple{nstate}, FT}(undef)
  l_aux = MArray{Tuple{nauxstate}, FT}(undef)
  l_knl = MArray{Tuple{nout, Nq}, FT}(undef)
  # note that k is the second not 4th index (since this is scratch memory and k
  # needs to be persistent across threads)
  l_int = @scratch FT (nout, Nq, Nq, Nqj) 2

  s_I = @shmem FT (Nq, Nq)

  @inbounds @loop for k in (1; threadIdx().z)
    @loop for i in (1:Nq; threadIdx().x)
      @unroll for n = 1:Nq
        s_I[i, n] = Imat[i, n]
      end
    end
  end
  @synchronize

  @inbounds @loop for eh in (elems; blockIdx().x)
    # Initialize the constant state at zero
    @loop for j in (1:Nqj; threadIdx().y)
      @loop for i in (1:Nq; threadIdx().x)
        @unroll for k in 1:Nq
          @unroll for s = 1:nout
            l_int[s, k, i, j] = 0
          end
        end
      end
    end
    # Loop up the stack of elements
    for ev = 1:nvertelem
      e = ev + (eh - 1) * nvertelem

      # Evaluate the integral kernel at each DOF in the slabk
      @loop for j in (1:Nqj; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          # loop up the pencil
          @unroll for k in 1:Nq
            ijk = i + Nq * ((j-1) + Nqj * (k-1))
            Jc = vgeo[ijk, _JcV, e]
            @unroll for s = 1:nstate
              l_Q[s] = Q[ijk, s, e]
            end

            @unroll for s = 1:nauxstate
              l_aux[s] = auxstate[ijk, s, e]
            end

            int_knl!(view(l_knl, :, k), l_Q, l_aux)

            # multiply in the curve jacobian
            @unroll for s = 1:nout
              l_knl[s, k] *= Jc
            end
          end

          # Evaluate the integral up the element
          @unroll for s = 1:nout
            @unroll for k in 1:Nq
              @unroll for n in 1:Nq
                l_int[s, k, i, j] += s_I[k, n] * l_knl[s, n]
              end
            end
          end

          # Store out to memory and reset the background value for next element
          @unroll for k in 1:Nq
            ijk = i + Nq * ((j-1) + Nqj * (k-1))
            @unroll for ind_out = 1:nout
              s = outstate[ind_out]
              P[ijk, s, e] = l_int[ind_out, k, i, j]
              l_int[ind_out, k, i, j] = l_int[ind_out, Nq, i, j]
            end
          end
        end
      end
    end
  end
  nothing
end

function knl_reverse_indefinite_stack_integral!(::Val{dim}, ::Val{N},
                                                ::Val{nvertelem}, P, elems,
                                                ::Val{outstate},
                                                ::Val{instate}
                                               ) where {dim, N, outstate,
                                                        instate, nvertelem}
  FT = eltype(P)

  Nq = N + 1
  Nqj = dim == 2 ? 1 : Nq

  nout = length(outstate)

  # note that k is the second not 4th index (since this is scratch memory and k
  # needs to be persistent across threads)
  l_T = MArray{Tuple{nout}, FT}(undef)
  l_V = MArray{Tuple{nout}, FT}(undef)

  @inbounds @loop for eh in (elems; blockIdx().x)
    # Initialize the constant state at zero
    @loop for j in (1:Nqj; threadIdx().y)
      @loop for i in (1:Nq; threadIdx().x)
        ijk = i + Nq * ((j-1) + Nqj * (Nq-1))
        et = nvertelem + (eh - 1) * nvertelem
        @unroll for s = 1:nout
          l_T[s] = P[ijk, instate[s], et]
        end

        # Loop up the stack of elements
        for ev = 1:nvertelem
          e = ev + (eh - 1) * nvertelem
          @unroll for k in 1:Nq
            ijk = i + Nq * ((j-1) + Nqj * (k-1))
            @unroll for s = 1:nout
              l_V[s] = P[ijk, instate[s], e]
            end
            @unroll for s = 1:nout
              P[ijk, outstate[s], e] = l_T[s] - l_V[s]
            end
          end
        end
      end
    end
  end
  nothing
end

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
@testset "wrapping_integers_unwrap.jl" begin
    @test unwrap(WrappingInt32(1), WrappingInt32(0), 0) == UInt(1)
    @test unwrap(WrappingInt32(1), WrappingInt32(0), typemax(UInt32)) == (UInt(1) << 32) + 1
    @test unwrap(WrappingInt32(typemax(UInt32) - 1), WrappingInt32(0),
                 3 * (UInt(1) << 32)) == 3 * (UInt(1) << 32) - 2
    @test unwrap(WrappingInt32(typemax(UInt32) - 10), WrappingInt32(0),
                 3 * (UInt(1) << 32)) == 3 * (UInt(1) << 32) - 11
    @test unwrap(WrappingInt32(typemax(UInt32)), WrappingInt32(10), 3 * (UInt(1) << 32)) ==
          3 * (UInt(1) << 32) - 11
    @test unwrap(WrappingInt32(typemax(UInt32)), WrappingInt32(0), 0) ==
          UInt(typemax(UInt32))
    @test unwrap(WrappingInt32(16), WrappingInt32(16), 0) == UInt(0)
    @test unwrap(WrappingInt32(15), WrappingInt32(16), 0) == UInt(typemax(UInt32))
    @test unwrap(WrappingInt32(0), WrappingInt32(typemax(Int32)), 0) ==
           UInt(typemax(Int32)) + 2
    @test unwrap(WrappingInt32(typemax(UInt32)), WrappingInt32(typemax(Int32)), 0) ==
           UInt(1) << 31
    @test unwrap(WrappingInt32(typemax(UInt32)), WrappingInt32(UInt(1) << 31), 0) ==
           UInt(typemax(UInt32)) >> 1
end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
# to be loaded into eFEMpart
# Assembles global linear system from local element matrices
# which are represented in functional form by 'localmat' argument

########################
###### Assembler #######
########################

"""
  assembleScalar(mesh,localmat,parameter)

Assembles sparse matrix for scalar PDEs

ASSUMPTIONS:
- mesh is either Mesh or FluidMesh type
- localmat is functional form of PDE for local element
- parameter can be either number or array (variable coefficient problem)

INPUT:
  mesh:       [AbstractMesh]
                  mesh type that PDE is being solved on
  localmat:   [Function]
                  function dictating what PDE is being solved
  parameter:  [Real or Vector{Real}]
                  parameter for PDE, can be either of type Real (constant coefficient
                    PDE) or Vector{Real} (variable coefficient PDE)

OUTPUT:
  -           [SparseArray]
                  sparse array for linear system to be solved
"""
function assembleScalar(mesh,localmat,parameter)
  # derive parameters
  nNodes = length(mesh.xy)
  nElm   = length(mesh.cm)
  (mesh.order==:Linear ? order=1 : order=2)
  (order==1 ? nNodesPerElm=4 : nNodesPerElm=9)
  nGaussNodes = nNodesPerElm  # possibly not true for axisymmetric quadrature

  #tempArr = Array{Float64}(nNodesPerElm)
  xN      = Array{Float64}(undef,nNodesPerElm)
  yN      = Array{Float64}(undef,nNodesPerElm)
  phi     = Array{Float64}(undef,nNodesPerElm)
  dphidx  = Array{Float64}(undef,nNodesPerElm)
  dphidy  = Array{Float64}(undef,nNodesPerElm)
  dphids  = Array{Float64}(undef,nNodesPerElm)
  dphidt  = Array{Float64}(undef,nNodesPerElm)

  At  = Array{Float64}(undef,nNodesPerElm,nNodesPerElm)

  I = Array{Int64}(undef,nElm,nNodesPerElm,nNodesPerElm)
  J = Array{Int64}(undef,nElm,nNodesPerElm,nNodesPerElm)
  S = Array{Float64}(undef,nElm,nNodesPerElm,nNodesPerElm)

  # compute phi,dphids,dphidt
  w,s,t = GaussQuadPoints2D(order+1)

  # construct local matrices
  for el = 1:nElm
    localmat(mesh,el,xN,yN,w,s,t,nGaussNodes,nNodesPerElm,
             At,phi,dphidx,dphidy,dphids,dphidt,order,parameter)
    for ti = 1:nNodesPerElm, tj = 1:nNodesPerElm
      I[el,ti,tj] = mesh.cm[el].NodeList[ti]
      J[el,ti,tj] = mesh.cm[el].NodeList[tj]
      S[el,ti,tj] = At[ti,tj]
    end
  end

  # assemble global matrix
  return sparse(vec(I),vec(J),vec(S),nNodes,nNodes)
end
# implements dirichlet boundary condition information before neumann
# and affects F simultaneously
# as per https://www.math.colostate.edu/~bangerth/videos.676.21.6.html
function assembleScalarandRHS_2b(mesh,localmat,parameter,prob)
  # derive parameters
  nNodes = length(mesh.xy)
  nElm   = length(mesh.cm)
  (mesh.order==:Linear ? order=1 : order=2)
  (order==1 ? nNodesPerElm=4 : nNodesPerElm=9)
  nGaussNodes = nNodesPerElm  # possibly not true for axisymmetric quadrature

  #tempArr = Array{Float64}(nNodesPerElm)
  xN      = Array{Float64}(undef,nNodesPerElm)
  yN      = Array{Float64}(undef,nNodesPerElm)
  phi     = Array{Float64}(undef,nNodesPerElm)
  dphidx  = Array{Float64}(undef,nNodesPerElm)
  dphidy  = Array{Float64}(undef,nNodesPerElm)
  dphids  = Array{Float64}(undef,nNodesPerElm)
  dphidt  = Array{Float64}(undef,nNodesPerElm)

  At  = Array{Float64}(undef,nNodesPerElm,nNodesPerElm)

  # indice arrays preallocation
  I = Array{Int64}(undef,nElm,nNodesPerElm,nNodesPerElm)
  J = Array{Int64}(undef,nElm,nNodesPerElm,nNodesPerElm)
  S = Array{Float64}(undef,nElm,nNodesPerElm,nNodesPerElm)

  # forcing array preallocation
  F = zeros(Float64,nNodes)

  # for indexing how many times boundary is accessed
  bcCOUNT = zeros(Float64,length(prob.bcid[:dNodes]))

  # compute phi,dphids,dphidt
  w,s,t = GaussQuadPoints2D(order+1)

  # construct local matrices
  # HERE IS WHERE YOU IMPLEMENT BOUNDARY DATA
  for el = 1:nElm
    localmat(mesh,el,xN,yN,w,s,t,nGaussNodes,nNodesPerElm,
             At,phi,dphidx,dphidy,dphids,dphidt,order,parameter)
    # local assembly -- nothing changes
    for ti = 1:nNodesPerElm, tj = 1:nNodesPerElm
      I[el,ti,tj] = mesh.cm[el].NodeList[ti]
      J[el,ti,tj] = mesh.cm[el].NodeList[tj]
      S[el,ti,tj] = At[ti,tj]
    end

    # local F assembly
    c = mesh.cm[el].NodeList
    for gpt=1:nNodesPerElm
      phi,_,_,jac = derivShape2D(s[gpt],t[gpt],xN,yN,order)
	    for i=1:nNodesPerElm
	      fcoefs      = prob.bcval[:forcing][mesh.cm[el].NodeList]
	      fgpt        = shapeEval(fcoefs,phi)
	      F[mesh.cm[el].NodeList[i]] += fgpt*phi[i]*w[gpt]*jac
	    end
	  end

    # check for boundary data & modify local stiffness and forcing arrays accordingly
    for ti=1:nNodesPerElm
      if mesh.cm[el].NodeList[ti] in prob.bcid[:dNodes]
        #find index of bcid
        bcIND = findfirst(x->x==mesh.cm[el].NodeList[ti], prob.bcid[:dNodes])
        bcCOUNT[bcIND] += 1

        # subtract matvec from F
        for tj=1:nNodesPerElm
          F[mesh.cm[el].NodeList[tj]] -= S[el,tj,ti]*prob.bcval[:dBC][bcIND]
          S[el,tj,ti] = 0.0
          S[el,ti,tj] = 0.0
        end

        # place 1 in diagonal of local stiffness
        S[el,ti,ti] = 1.0

        # replace Frocing with dirichlet value
        F[mesh.cm[el].NodeList[ti]] = prob.bcval[:dBC][bcIND]
      end
    end
  end

  F[prob.bcid[:dNodes]] .*= bcCOUNT

  # assemble global matrix
  return sparse(vec(I),vec(J),vec(S),nNodes,nNodes), F
end

"""
  assembleHalfFluid(mesh,localmat,parameter...)

Assembles sparse matrix for Fluids with no cross-terms

ASSUMPTIONS:
- mesh is either Mesh or FluidMesh type
- localmat is functional form of PDE for local element
- parameter can be either number or array (variable coefficient problem),
  also can be a tuple of parameters depending on expectation of localmat

INPUT:
  mesh:       [AbstractMesh]
                  mesh type that PDE is being solved on
  localmat:   [Function]
                  function dictating what PDE is being solved
  parameter:  [Real or Vector{Real}]
                  parameter for PDE, can be either of type Real (constant coefficient
                    PDE) or Vector{Real} (variable coefficient PDE)

OUTPUT:
  -           [SparseArray]
                  sparse array for linear system to be solved
"""
function assembleHalfFluid(mesh,localmat,parameter...)
  nElm = length(mesh.cm)
  nUNodes = length(mesh.xy); nPNodes = length(mesh.xyp)
  nGaussNodes = 9
  nNodesPerElmQ1 = 4
  nNodesPerElmQ2 = 9

  xN = Array{Float64}(undef,nNodesPerElmQ2)
  yN = Array{Float64}(undef,nNodesPerElmQ2)
  phi    = Array{Float64}(undef,nNodesPerElmQ2)
  dphids = Array{Float64}(undef,nNodesPerElmQ2)
  dphidt = Array{Float64}(undef,nNodesPerElmQ2)
  dphidx = Array{Float64}(undef,nNodesPerElmQ2)
  dphidy = Array{Float64}(undef,nNodesPerElmQ2)
  psi    = Array{Float64}(undef,nNodesPerElmQ1)
  dpsids = Array{Float64}(undef,nNodesPerElmQ1)
  dpsidt = Array{Float64}(undef,nNodesPerElmQ1)
  dpsidx = Array{Float64}(undef,nNodesPerElmQ1)
  dpsidy = Array{Float64}(undef,nNodesPerElmQ1)

  Axt = Array{Float64}(undef,nNodesPerElmQ2,nNodesPerElmQ2)
  Ayt = Array{Float64}(undef,nNodesPerElmQ2,nNodesPerElmQ2)
  Bxt = Array{Float64}(undef,nNodesPerElmQ1,nNodesPerElmQ2)
  Byt = Array{Float64}(undef,nNodesPerElmQ1,nNodesPerElmQ2)

  IA1 = zeros(Float64,nElm,nNodesPerElmQ2,nNodesPerElmQ2)
  JA1 = zeros(Float64,nElm,nNodesPerElmQ2,nNodesPerElmQ2)
  SA1 = zeros(Float64,nElm,nNodesPerElmQ2,nNodesPerElmQ2)
  IA2 = zeros(Float64,nElm,nNodesPerElmQ2,nNodesPerElmQ2)
  JA2 = zeros(Float64,nElm,nNodesPerElmQ2,nNodesPerElmQ2)
  SA2 = zeros(Float64,nElm,nNodesPerElmQ2,nNodesPerElmQ2)
  IBx1 = zeros(Float64,nElm,nNodesPerElmQ1,nNodesPerElmQ2)
  JBx1 = zeros(Float64,nElm,nNodesPerElmQ1,nNodesPerElmQ2)
  SBx1 = zeros(Float64,nElm,nNodesPerElmQ1,nNodesPerElmQ2)
  IBy1 = zeros(Float64,nElm,nNodesPerElmQ1,nNodesPerElmQ2)
  JBy1 = zeros(Float64,nElm,nNodesPerElmQ1,nNodesPerElmQ2)
  SBy1 = zeros(Float64,nElm,nNodesPerElmQ1,nNodesPerElmQ2)
  IBx2 = zeros(Float64,nElm,nNodesPerElmQ1,nNodesPerElmQ2)
  JBx2 = zeros(Float64,nElm,nNodesPerElmQ1,nNodesPerElmQ2)
  SBx2 = zeros(Float64,nElm,nNodesPerElmQ1,nNodesPerElmQ2)
  IBy2 = zeros(Float64,nElm,nNodesPerElmQ1,nNodesPerElmQ2)
  JBy2 = zeros(Float64,nElm,nNodesPerElmQ1,nNodesPerElmQ2)
  SBy2 = zeros(Float64,nElm,nNodesPerElmQ1,nNodesPerElmQ2)

  # for axisymmetric, better to use 4x4 gauss?
  w,s,t = GaussQuadPoints2D(3)

  for el = 1:nElm
    localmat(mesh,el,xN,yN,w,s,t,Axt,Ayt,Bxt,Byt,
              phi,dphidx,dphidy,dphids,dphidt,
              psi,dpsidx,dpsidy,dpsids,dpsidt,
              parameter...)

    for ti = 1:nNodesPerElmQ2, tj = 1:nNodesPerElmQ2
      # Axt
      IA1[el,ti,tj] = mesh.cm[el].NodeList[ti]
      JA1[el,ti,tj] = mesh.cm[el].NodeList[tj]
      SA1[el,ti,tj] = Axt[ti,tj]
      # Ayt
      IA2[el,ti,tj] = mesh.cm[el].NodeList[ti] + nUNodes
      JA2[el,ti,tj] = mesh.cm[el].NodeList[tj] + nUNodes
      SA2[el,ti,tj] = Ayt[ti,tj]
    end
    for tk = 1:nNodesPerElmQ1, ti=1:nNodesPerElmQ2
      # Bx
      IBx1[el,tk,ti] = mesh.cmp[el].NodeList[tk] + 2*nUNodes
      JBx1[el,tk,ti] = mesh.cm[el].NodeList[ti]
      SBx1[el,tk,ti] = Bxt[tk,ti]
      # By
      IBy1[el,tk,ti] = mesh.cmp[el].NodeList[tk] + 2*nUNodes
      JBy1[el,tk,ti] = mesh.cm[el].NodeList[ti]  + nUNodes
      SBy1[el,tk,ti] = Byt[tk,ti]
      # Bx'
      IBx2[el,tk,ti] = mesh.cm[el].NodeList[ti]
      JBx2[el,tk,ti] = mesh.cmp[el].NodeList[tk] + 2*nUNodes
      SBx2[el,tk,ti] = Bxt[tk,ti]
      # By'
      IBy2[el,tk,ti] = mesh.cm[el].NodeList[ti]  + nUNodes
      JBy2[el,tk,ti] = mesh.cmp[el].NodeList[tk] + 2*nUNodes
      SBy2[el,tk,ti] = Byt[tk,ti]
    end
  end

  I = vcat(vec(IA1),vec(IA2),vec(IBx1),vec(IBx2),vec(IBy1),vec(IBy2))
  J = vcat(vec(JA1),vec(JA2),vec(JBx1),vec(JBx2),vec(JBy1),vec(JBy2))
  S = vcat(vec(SA1),vec(SA2),vec(SBx1),vec(SBx2),vec(SBy1),vec(SBy2))

  return sparse(I,J,S,2*nUNodes+nPNodes,2*nUNodes+nPNodes)
end
"""
  assembleFullFluid_WRONG(mesh,localmat,parameter...)

actually not wrong... ?????
"""
function assembleFullFluid_WRONG(mesh,localmat,parameter...)
  nElm = length(mesh.cm)
  nUNodes = length(mesh.xy); nPNodes = length(mesh.xyp)
  nGaussNodes = 9
  nNodesPerElmQ1 = 4
  nNodesPerElmQ2 = 9

  w,s,t = GaussQuadPoints2D(3)

  tempArr4 = Array{Float64}(undef,nNodesPerElmQ1)
  tempArr9 = Array{Float64}(undef,nNodesPerElmQ2)

  xN     = copy(tempArr9); yN     = copy(tempArr9)

  phi    = copy(tempArr9); dphids = copy(tempArr9); dphidt=copy(tempArr9)
  dphidx = copy(tempArr9); dphidy = copy(tempArr9)

  psi    = copy(tempArr4); dpsids = copy(tempArr4); dpsidt = copy(tempArr4)
  dpsidx = copy(tempArr4); dpsidy = copy(tempArr4)

  # top left
  Axt  = Array{Float64}(undef,nNodesPerElmQ2,nNodesPerElmQ2)
  # top right
  A12t = Array{Float64}(undef,nNodesPerElmQ2,nNodesPerElmQ2)
  # bottom left
  A21t = Array{Float64}(undef,nNodesPerElmQ2,nNodesPerElmQ2)
  # bottom right
  Ayt  = Array{Float64}(undef,nNodesPerElmQ2,nNodesPerElmQ2)
  Bxt = Array{Float64}(undef,nNodesPerElmQ1,nNodesPerElmQ2)
  Byt = Array{Float64}(undef,nNodesPerElmQ1,nNodesPerElmQ2)

  IA11 = zeros(Float64,nElm,nNodesPerElmQ2,nNodesPerElmQ2)
  JA11 = copy(IA11); SA11 = copy(IA11);
  IA12 = copy(IA11); JA12 = copy(IA11); SA12 = copy(IA11)
  IA21 = copy(IA11); JA21 = copy(IA11); SA21 = copy(IA11)
  IA22 = copy(IA11); JA22 = copy(IA11); SA22 = copy(IA11)

  IBx1 = zeros(Float64,nElm,nNodesPerElmQ1,nNodesPerElmQ2)
  JBx1 = copy(IBx1); SBx1 = copy(IBx1)
  IBy1 = copy(IBx1); JBy1 = copy(IBx1); SBy1 = copy(IBx1)
  IBx2 = copy(IBx1); JBx2 = copy(IBx1); SBx2 = copy(IBx1)
  IBy2 = copy(IBx1); JBy2 = copy(IBx1); SBy2 = copy(IBx1)

  for el = 1:nElm
    localmat(mesh,el,xN,yN,w,s,t,Axt,Ayt,A12t,A21t,Bxt,Byt,
              phi,dphidx,dphidy,dphids,dphidt,
              psi,dpsidx,dpsidy,dpsids,dpsidt,
              parameter...)

    for ti = 1:nNodesPerElmQ2, tj = 1:nNodesPerElmQ2
      # Axt
      IA11[el,ti,tj] = mesh.cm[el].NodeList[ti]
      JA11[el,ti,tj] = mesh.cm[el].NodeList[tj]
      SA11[el,ti,tj] = Axt[ti,tj]
      # A12t
      IA12[el,ti,tj] = mesh.cm[el].NodeList[ti]
      JA12[el,ti,tj] = mesh.cm[el].NodeList[tj] + nUNodes
      SA12[el,ti,tj] = A12t[ti,tj]
      # A21t
      IA21[el,ti,tj] = mesh.cm[el].NodeList[ti] + nUNodes
      JA21[el,ti,tj] = mesh.cm[el].NodeList[tj]
      SA21[el,ti,tj] = A21t[ti,tj]
      # Ayt
      IA22[el,ti,tj] = mesh.cm[el].NodeList[ti] + nUNodes
      JA22[el,ti,tj] = mesh.cm[el].NodeList[tj] + nUNodes
      SA22[el,ti,tj] = Ayt[ti,tj]
    end
    for tk = 1:nNodesPerElmQ1, ti=1:nNodesPerElmQ2
      # Bx
      IBx1[el,tk,ti] = mesh.cmp[el].NodeList[tk] + 2*nUNodes
      JBx1[el,tk,ti] = mesh.cm[el].NodeList[ti]
      SBx1[el,tk,ti] = Bxt[tk,ti]
      # By
      IBy1[el,tk,ti] = mesh.cmp[el].NodeList[tk] + 2*nUNodes
      JBy1[el,tk,ti] = mesh.cm[el].NodeList[ti]  + nUNodes
      SBy1[el,tk,ti] = Byt[tk,ti]
      # Bx'
      IBx2[el,tk,ti] = mesh.cm[el].NodeList[ti]
      JBx2[el,tk,ti] = mesh.cmp[el].NodeList[tk] + 2*nUNodes
      SBx2[el,tk,ti] = Bxt[tk,ti]
      # By'
      IBy2[el,tk,ti] = mesh.cm[el].NodeList[ti]  + nUNodes
      JBy2[el,tk,ti] = mesh.cmp[el].NodeList[tk] + 2*nUNodes
      SBy2[el,tk,ti] = Byt[tk,ti]
    end
  end

  I = vcat(vec(IA11),vec(IA12),vec(IA21),vec(IA22),vec(IBx1),vec(IBx2),vec(IBy1),vec(IBy2))
  J = vcat(vec(JA11),vec(JA12),vec(JA21),vec(JA22),vec(JBx1),vec(JBx2),vec(JBy1),vec(JBy2))
  S = vcat(vec(SA11),vec(SA12),vec(SA21),vec(SA22),vec(SBx1),vec(SBx2),vec(SBy1),vec(SBy2))

  return sparse(I,J,S,2*nUNodes+nPNodes,2*nUNodes+nPNodes)
end
"""
  assembleFullFluid(mesh,localmat,parameter...)

Assembles full 3x3 block matrix for fluid problems
"""
# Assembles full 3x3 block matrix for fluid problems
function assembleFullFluid(mesh,localmat,parameter...)
  nElm = length(mesh.cm)
  nUNodes = length(mesh.xy); nPNodes = length(mesh.xyp)
  nGaussNodes = 9
  nNodesPerElmQ1 = 4
  nNodesPerElmQ2 = 9

  w,s,t = GaussQuadPoints2D(3)

  tempArr4 = Array{Float64}(undef,nNodesPerElmQ1)
  tempArr9 = Array{Float64}(undef,nNodesPerElmQ2)

  xN     = copy(tempArr9); yN     = copy(tempArr9)

  phi    = copy(tempArr9); dphids = copy(tempArr9); dphidt=copy(tempArr9)
  dphidx = copy(tempArr9); dphidy = copy(tempArr9)

  psi    = copy(tempArr4); dpsids = copy(tempArr4); dpsidt = copy(tempArr4)
  dpsidx = copy(tempArr4); dpsidy = copy(tempArr4)

  # generate temporary element block matrices
  At = Array{Float64}(undef,nNodesPerElmQ2,nNodesPerElmQ2)
  Bt = Array{Float64}(undef,nNodesPerElmQ2,nNodesPerElmQ2)
  Ct = Array{Float64}(undef,nNodesPerElmQ2,nNodesPerElmQ1)
  Dt = Array{Float64}(undef,nNodesPerElmQ2,nNodesPerElmQ2)
  Et = Array{Float64}(undef,nNodesPerElmQ2,nNodesPerElmQ2)
  Ft = Array{Float64}(undef,nNodesPerElmQ2,nNodesPerElmQ1)
  Gt = Array{Float64}(undef,nNodesPerElmQ1,nNodesPerElmQ2)
  Ht = Array{Float64}(undef,nNodesPerElmQ1,nNodesPerElmQ2)

  IA = zeros(Float64,nElm,nNodesPerElmQ2,nNodesPerElmQ2)
  JA = copy(IA); SA = copy(IA)
  IB = copy(IA); JB = copy(IA); SB = copy(IA)
  ID = copy(IA); JD = copy(IA); SD = copy(IA)
  IE = copy(IA); JE = copy(IA); SE = copy(IA)

  IG = zeros(Float64,nElm,nNodesPerElmQ1,nNodesPerElmQ2)
  JG = copy(IG); SG = copy(IG)
  IH = copy(IG); JH = copy(IG); SH = copy(IG)

  # transpose!!!!
  IC = zeros(Float64,nElm,nNodesPerElmQ1,nNodesPerElmQ2)
  JC = copy(IC); SC = copy(IC)
  IF = copy(IC); JF = copy(IC); SF = copy(IC)

  BigMat = spzeros(2*nUNodes+nPNodes,2*nUNodes+nPNodes)

  for el = 1:nElm
    localmat(mesh,el,xN,yN,w,s,t,
             At,Bt,Ct,Dt,Et,Ft,Gt,Ht,
             phi,dphidx,dphidy,dphids,dphidt,
             psi,dpsidx,dpsidy,dpsids,dpsidt,
             parameter...)
    # big block matrices -- A, B, D, E
    for ti = 1:nNodesPerElmQ2, tj = 1:nNodesPerElmQ2
      # At
      iii = mesh.cm[el].NodeList[ti]
      jjj = mesh.cm[el].NodeList[tj]
      BigMat[iii,jjj] += At[ti,tj]

      # Bt
      iii = mesh.cm[el].NodeList[ti]
      jjj = mesh.cm[el].NodeList[tj] + nUNodes
      BigMat[iii,jjj] += Bt[ti,tj]

      # Dt
      iii = mesh.cm[el].NodeList[ti] + nUNodes
      jjj = mesh.cm[el].NodeList[tj]
      BigMat[iii,jjj] += Dt[ti,tj]

      # Et
      iii = mesh.cm[el].NodeList[ti] + nUNodes
      jjj = mesh.cm[el].NodeList[tj] + nUNodes
      BigMat[iii,jjj] += Et[ti,tj]

      # At
      #IA[el,ti,tj] = mesh.cm[el].NodeList[ti]
      #JA[el,ti,tj] = mesh.cm[el].NodeList[tj]
      #SA[el,ti,tj] = At[ti,tj]
      # Bt
      #IB[el,ti,tj] = mesh.cm[el].NodeList[ti]
      #JB[el,ti,tj] = mesh.cm[el].NodeList[tj] + nUNodes
      #SB[el,ti,tj] = Bt[ti,tj]
      # Dt
      #ID[el,ti,tj] = mesh.cm[el].NodeList[ti] + nUNodes
      #JD[el,ti,tj] = mesh.cm[el].NodeList[tj]
      #SD[el,ti,tj] = Dt[ti,tj]
      # Et
      #IE[el,ti,tj] = mesh.cm[el].NodeList[ti] + nUNodes
      #JE[el,ti,tj] = mesh.cm[el].NodeList[tj] + nUNodes
      #SE[el,ti,tj] = Et[ti,tj]
    end
    # bottom block matrices -- G, H
    for ti = 1:nNodesPerElmQ1, tj = 1:nNodesPerElmQ2
      # Gt
      iii = mesh.cmp[el].NodeList[ti] + 2*nUNodes
      jjj = mesh.cm[el].NodeList[tj]
      BigMat[iii,jjj] += Gt[ti,tj]

      # Ht
      iii = mesh.cmp[el].NodeList[ti] + 2*nUNodes
      jjj = mesh.cm[el].NodeList[tj]  + nUNodes
      BigMat[iii,jjj] += Ht[ti,tj]

      # Ct
      iii = mesh.cm[el].NodeList[tj]
      jjj = mesh.cmp[el].NodeList[ti] + 2*nUNodes
      BigMat[iii,jjj] += Ct[tj,ti]

      # Ft
      iii = mesh.cm[el].NodeList[tj]  + nUNodes
      jjj = mesh.cmp[el].NodeList[ti] + 2*nUNodes
      BigMat[iii,jjj] += Ft[tj,ti]

      # Gt
      #IG[el,ti,tj] = mesh.cmp[el].NodeList[ti] + 2*nUNodes
      #JG[el,ti,tj] = mesh.cm[el].NodeList[tj]
      #SG[el,ti,tj] = Gt[ti,tj]
      # Ht
      #IH[el,ti,tj] = mesh.cmp[el].NodeList[ti] + 2*nUNodes
      #JH[el,ti,tj] = mesh.cm[el].NodeList[tj]  + nUNodes
      #SH[el,ti,tj] = Ht[ti,tj]
      # Ct
      #IC[el,ti,tj] = mesh.cm[el].NodeList[tj]
      #JC[el,ti,tj] = mesh.cmp[el].NodeList[ti] + 2*nUNodes
      #SC[el,ti,tj] = Ct[tj,ti]
      # Ft
      #IF[el,ti,tj] = mesh.cm[el].NodeList[tj]  + nUNodes
      #JF[el,ti,tj] = mesh.cmp[el].NodeList[ti] + 2*nUNodes
      #SF[el,ti,tj] = Ft[tj,ti]
    end
  end

  # bug could be in how this is vectorized...
  #I = vcat(vec(IA),vec(IB),vec(IC),vec(ID),vec(IE),vec(IF),vec(IG),vec(IH))
  #J = vcat(vec(JA),vec(JB),vec(JC),vec(JD),vec(JE),vec(JF),vec(JG),vec(JH))
  #S = vcat(vec(SA),vec(SB),vec(SC),vec(SD),vec(SE),vec(SF),vec(SG),vec(SH))

  #return sparse(I,J,S,2*nUNodes+nPNodes,2*nUNodes+nPNodes)
  return BigMat
end
"""
  assembleMPB(mesh,localmat,parameter...)

Assembles full 3x3 block matrix for fluid problems. Not sure how it differs
  from the above.....
"""
function assembleMPB(mesh,localmat,parameter...)
  nElm = length(mesh.cm)
  nUNodes = length(mesh.xy); nPNodes = length(mesh.xyp)
  nGaussNodes = 9
  nNodesPerElmQ1 = 4
  nNodesPerElmQ2 = 9

  w,s,t = GaussQuadPoints2D(3)

  tempArr4 = Array{Float64}(undef,nNodesPerElmQ1)
  tempArr9 = Array{Float64}(undef,nNodesPerElmQ2)

  xN     = copy(tempArr9); yN     = copy(tempArr9)

  phi    = copy(tempArr9); dphids = copy(tempArr9); dphidt=copy(tempArr9)
  dphidx = copy(tempArr9); dphidy = copy(tempArr9)

  psi    = copy(tempArr4); dpsids = copy(tempArr4); dpsidt = copy(tempArr4)
  dpsidx = copy(tempArr4); dpsidy = copy(tempArr4)

  # generate temporary element block matrices
  At = Array{Float64}(undef,nNodesPerElmQ2,nNodesPerElmQ2)
  Ct = Array{Float64}(undef,nNodesPerElmQ2,nNodesPerElmQ1)
  Et = Array{Float64}(undef,nNodesPerElmQ2,nNodesPerElmQ2)
  Ft = Array{Float64}(undef,nNodesPerElmQ2,nNodesPerElmQ1)
  Gt = Array{Float64}(undef,nNodesPerElmQ1,nNodesPerElmQ2)
  Ht = Array{Float64}(undef,nNodesPerElmQ1,nNodesPerElmQ2)

  IA = zeros(Float64,nElm,nNodesPerElmQ2,nNodesPerElmQ2)
  JA = copy(IA); SA = copy(IA)
  IE = copy(IA); JE = copy(IA); SE = copy(IA)

  IG = zeros(Float64,nElm,nNodesPerElmQ1,nNodesPerElmQ2)
  JG = copy(IG); SG = copy(IG)
  IH = copy(IG); JH = copy(IG); SH = copy(IG)

  # transpose!!!!
  IC = zeros(Float64,nElm,nNodesPerElmQ1,nNodesPerElmQ2)
  JC = copy(IC); SC = copy(IC)
  IF = copy(IC); JF = copy(IC); SF = copy(IC)

  for el = 1:nElm
    localmat(mesh,el,xN,yN,w,s,t,
             At,Ct,Et,Ft,Gt,Ht,
             phi,dphidx,dphidy,dphids,dphidt,
             psi,dpsidx,dpsidy,dpsids,dpsidt,
             parameter...)
    # big block matrices -- A, B, D, E
    for ti = 1:nNodesPerElmQ2, tj = 1:nNodesPerElmQ2
      # At
      IA[el,ti,tj] = mesh.cm[el].NodeList[ti]
      JA[el,ti,tj] = mesh.cm[el].NodeList[tj]
      SA[el,ti,tj] = At[ti,tj]
      # Et
      IE[el,ti,tj] = mesh.cm[el].NodeList[ti] + nUNodes
      JE[el,ti,tj] = mesh.cm[el].NodeList[tj] + nUNodes
      SE[el,ti,tj] = Et[ti,tj]
    end
    # bottom block matrices -- G, H
    for ti = 1:nNodesPerElmQ1, tj = 1:nNodesPerElmQ2
      # Gt
      IG[el,ti,tj] = mesh.cmp[el].NodeList[ti] + 2*nUNodes
      JG[el,ti,tj] = mesh.cm[el].NodeList[tj]
      SG[el,ti,tj] = Gt[ti,tj]
      # Ht
      IH[el,ti,tj] = mesh.cmp[el].NodeList[ti] + 2*nUNodes
      JH[el,ti,tj] = mesh.cm[el].NodeList[tj]  + nUNodes
      SH[el,ti,tj] = Ht[ti,tj]
      # Ct
      IC[el,ti,tj] = mesh.cm[el].NodeList[tj]
      JC[el,ti,tj] = mesh.cmp[el].NodeList[ti] + 2*nUNodes
      SC[el,ti,tj] = Ct[tj,ti]
      # Ft
      IF[el,ti,tj] = mesh.cm[el].NodeList[tj]  + nUNodes
      JF[el,ti,tj] = mesh.cmp[el].NodeList[ti] + 2*nUNodes
      SF[el,ti,tj] = Ft[tj,ti]
    end
  end

  # bug could be in how this is vectorized...
  I = vcat(vec(IA),vec(IC),vec(IE),vec(IF),vec(IG),vec(IH))
  J = vcat(vec(JA),vec(JC),vec(JE),vec(JF),vec(JG),vec(JH))
  S = vcat(vec(SA),vec(SC),vec(SE),vec(SF),vec(SG),vec(SH))

  return sparse(I,J,S,2*nUNodes+nPNodes,2*nUNodes+nPNodes)
end

"""
WeakScalar2D(mesh,farr)

Computes the weak form of some function farr using cartesian integral

ASSUMPTIONS:
- mesh is either Mesh or FluidMesh type

INPUT:
  mesh:       [AbstractMesh]
                  mesh type that PDE is being solved on
  farr:       [Vector{Real}]
                  forcing values at mesh nodes

OUTPUT:
  F           [Vector{Real}]
                  RHS vector in form of RHS of linear Finite Element system
"""
function WeakScalar2D(mesh,farr)
  Nnodes     = length(mesh.xy)
  Nel        = length(mesh.cm)
  NumElNodes = length(mesh.cm[1].NodeList)
  F          = zeros(Nnodes)

  gaussorder = 0
  NumGaussNodes = 0
  if mesh.order == :Linear
    gaussorder = 2
    NumGaussNodes = Int(gaussorder^2)
  elseif mesh.order == :Quadratic
    gaussorder = 3
    NumGaussNodes = Int(gaussorder^2)
  end

  xNodes = zeros(Float64,NumGaussNodes)
  yNodes = zeros(Float64,NumGaussNodes)

  # Generate local F by elements
  # generate basis values at gauss points
  w,s,t = GaussQuadPoints2D(gaussorder)
  for el=1:Nel
    c = mesh.cm[el].NodeList
    for i=1:NumGaussNodes
      xNodes[i] = mesh.xy[mesh.cm[el].NodeList[i]].x
      yNodes[i] = mesh.xy[mesh.cm[el].NodeList[i]].y
    end

    for gpt=1:NumElNodes
      phi,_,_,jac = derivShape2D(s[gpt],t[gpt],xNodes,yNodes,gaussorder-1)

	    for i=1:NumElNodes
	      fcoefs      = farr[c]
	      fgpt        = shapeEval(fcoefs,phi)
	      F[mesh.cm[el].NodeList[i]] += fgpt*phi[i]*w[gpt]*jac
	    end
	  end
  end

  return F
end

"""
  WeakScalarAS(mesh,farr)

Computes the weak form of some function farr using axisymmetric integral

ASSUMPTIONS:
- mesh is either Mesh or FluidMesh type
- axisymmetric coordinates are accounted for.
  Expects farr without modification by e.g. radial coordinate

INPUT:
  mesh:       [AbstractMesh]
                  mesh type that PDE is being solved on
  farr:       [Vector{Real}]
                  forcing values at mesh nodes

OUTPUT:
  F           [Vector{Real}]
                  RHS vector in form of RHS of linear Finite Element system
                    in axisymmetric coordinates
"""
function WeakScalarAS(mesh,farr)
  Nnodes     = length(mesh.xy)
  Nel        = length(mesh.cm)
  NumElNodes = length(mesh.cm[1].NodeList)
  F          = zeros(Nnodes)

  gaussorder = 0
  NumGaussNodes = 0
  if mesh.order == :Linear
    gaussorder = 2
    NumGaussNodes = Int(gaussorder^2)
  elseif mesh.order == :Quadratic
    gaussorder = 3
    NumGaussNodes = Int(gaussorder^2)
  end

  xNodes = zeros(Float64,NumGaussNodes)
  yNodes = zeros(Float64,NumGaussNodes)

  # Generate local F by elements
  # generate basis values at gauss points
  w,s,t = GaussQuadPoints2D(gaussorder)
  for el=1:Nel
    c = mesh.cm[el].NodeList
    for i=1:NumGaussNodes
      xNodes[i] = mesh.xy[mesh.cm[el].NodeList[i]].x
      yNodes[i] = mesh.xy[mesh.cm[el].NodeList[i]].y
    end

    for gpt=1:NumElNodes
      phi,_,_,jac = derivShape2D(s[gpt],t[gpt],xNodes,yNodes,gaussorder-1)

	    for i=1:NumElNodes
	      fcoefs = farr[c]
        fgpt   = shapeEval(fcoefs,phi)
        yg     = shapeEval(yNodes,phi)
	      F[mesh.cm[el].NodeList[i]] += fgpt*phi[i]*w[gpt]*jac*yg
	    end
	  end
  end

  return F
end

# generates SUPG adjustment for axisymmetric cases
# parameter is of type for advection-diffusion, includes velocity field
# expects parameter to be a constant
function WeakSUPG2D(mesh,farr,parameter::T) where T<:AbstractConstantParameter
  Nnodes     = length(mesh.xy)
  Nel        = length(mesh.cm)
  NumElNodes = length(mesh.cm[1].NodeList)
  F          = zeros(Nnodes)

  # define 'wind'
  U = parameter.u
  V = parameter.v

  gaussorder = 0
  NumGaussNodes = 0
  if mesh.order == :Linear
    gaussorder = 2
    NumGaussNodes = Int(gaussorder^2)
  elseif mesh.order == :Quadratic
    gaussorder = 3
    NumGaussNodes = Int(gaussorder^2)
  end

  xNodes  = zeros(Float64,NumGaussNodes)
  yNodes  = zeros(Float64,NumGaussNodes)
  UNodes  = zeros(Float64,NumGaussNodes)
  VNodes  = zeros(Float64,NumGaussNodes)
  phi     = Array{Float64}(undef,NumElNodes)
  dphidx  = Array{Float64}(undef,NumElNodes)
  dphidy  = Array{Float64}(undef,NumElNodes)
  dphids  = Array{Float64}(undef,NumElNodes)
  dphidt  = Array{Float64}(undef,NumElNodes)

  # Generate local F by elements
  # generate basis values at gauss points
  w,s,t = GaussQuadPoints2D(gaussorder)
  for el=1:Nel
    c = mesh.cm[el].NodeList
    for i=1:NumGaussNodes
      xNodes[i] = mesh.xy[mesh.cm[el].NodeList[i]].x
      yNodes[i] = mesh.xy[mesh.cm[el].NodeList[i]].y
      UNodes[i] = U[c[i]]
      VNodes[i] = V[c[i]]
    end

    # compute h (max element length)
    hk = sqrt(areaCalc(xNodes,yNodes))

    # compute element Peclet number
    wmag = 0.0
    elPe = 0.0
    if mesh.order == :Linear
      uavg = sum(U[c])/4
      vavg = sum(V[c])/4
      wmag = sqrt(uavg^2 + vavg^2)

      elPe = 0.5*wmag*hk/parameter.κ
    elseif mesh.order == :Quadratic
      wmag = sqrt(U[c[9]]^2 + V[c[9]]^2)

      elPe = 0.5*wmag*hk/parameter.κ
    else
      error("mesh.order not specified correctly")
    end

    # compute delta for SUPG
    δ = 0.5*hk/wmag*(1.0 - 1.0/elPe)

    for gpt=1:NumElNodes
      shape2D!(phi,dphids,dphidt,s[gpt],t[gpt],gaussorder-1)
      derivShape2D!(phi,dphidx,dphidy,dphids,dphidt,s[gpt],t[gpt],xNodes,yNodes,gaussorder-1)
      ug  = shapeEval(UNodes,phi)
      vg  = shapeEval(VNodes,phi)

	    for i=1:NumElNodes
	      fcoefs = farr[c]
        fgpt   = shapeEval(fcoefs,phi)
        yg     = shapeEval(yNodes,phi)
	      F[mesh.cm[el].NodeList[i]] += δ*fgpt*(ug*dphidx[i] + vg*dphidy[i])*w[gpt]
	    end
	  end
  end

  return F
end

# generates SUPG adjustment for axisymmetric cases
# parameter is of type for advection-diffusion, includes velocity field
# expects parameter to be a constant
function WeakSUPGAS(mesh,farr,parameter::T) where T<:AbstractConstantParameter
  Nnodes     = length(mesh.xy)
  Nel        = length(mesh.cm)
  NumElNodes = length(mesh.cm[1].NodeList)
  F          = zeros(Nnodes)

  # define 'wind'
  U = parameter.u
  V = parameter.v

  gaussorder = 0
  NumGaussNodes = 0
  if mesh.order == :Linear
    gaussorder = 2
    NumGaussNodes = Int(gaussorder^2)
  elseif mesh.order == :Quadratic
    gaussorder = 3
    NumGaussNodes = Int(gaussorder^2)
  end

  xNodes  = zeros(Float64,NumGaussNodes)
  yNodes  = zeros(Float64,NumGaussNodes)
  UNodes  = zeros(Float64,NumGaussNodes)
  VNodes  = zeros(Float64,NumGaussNodes)
  phi     = Array{Float64}(undef,NumElNodes)
  dphidx  = Array{Float64}(undef,NumElNodes)
  dphidy  = Array{Float64}(undef,NumElNodes)
  dphids  = Array{Float64}(undef,NumElNodes)
  dphidt  = Array{Float64}(undef,NumElNodes)

  # Generate local F by elements
  # generate basis values at gauss points
  w,s,t = GaussQuadPoints2D(gaussorder)
  for el=1:Nel
    c = mesh.cm[el].NodeList
    for i=1:NumGaussNodes
      xNodes[i] = mesh.xy[mesh.cm[el].NodeList[i]].x
      yNodes[i] = mesh.xy[mesh.cm[el].NodeList[i]].y
      UNodes[i] = U[c[i]]
      VNodes[i] = V[c[i]]
    end

    # compute h (max element length)
    hk = sqrt(areaCalc(xNodes,yNodes))

    # compute element Peclet number
    wmag = 0.0
    elPe = 0.0
    if mesh.order == :Linear
      uavg = sum(U[c])/4
      vavg = sum(V[c])/4
      wmag = sqrt(uavg^2 + vavg^2)

      elPe = 0.5*wmag*hk/parameter.κ
    elseif mesh.order == :Quadratic
      wmag = sqrt(U[c[9]]^2 + V[c[9]]^2)

      elPe = 0.5*wmag*hk/parameter.κ
    else
      error("mesh.order not specified correctly")
    end

    # compute delta for SUPG
    δ = 0.5*hk/wmag*(1.0 - 1.0/elPe)

    for gpt=1:NumElNodes
      shape2D!(phi,dphids,dphidt,s[gpt],t[gpt],gaussorder-1)
      derivShape2D!(phi,dphidx,dphidy,dphids,dphidt,s[gpt],t[gpt],xNodes,yNodes,gaussorder-1)
      ug  = shapeEval(UNodes,phi)
      vg  = shapeEval(VNodes,phi)

	    for i=1:NumElNodes
	      fcoefs = farr[c]
        fgpt   = shapeEval(fcoefs,phi)
        yg     = shapeEval(yNodes,phi)
	      F[mesh.cm[el].NodeList[i]] += δ*fgpt*(ug*dphidx[i] + vg*dphidy[i])*w[gpt]*yg
	    end
	  end
  end

  return F
end

# written in BilinearForms.jl
# given 4 nodes, computes area
#function areaCalc(xNodes,yNodes)
#  area1 = xNodes[1]*(yNodes[2] - yNodes[3]) +
#          xNodes[2]*(yNodes[3] - yNodes[1]) +
#          xNodes[3]*(yNodes[1] - yNodes[2])
#  area2 = xNodes[1]*(yNodes[4] - yNodes[3]) +
#          xNodes[4]*(yNodes[3] - yNodes[1]) +
#          xNodes[3]*(yNodes[1] - yNodes[4])
#
#  return 0.5*(abs(area1) + abs(area2))
#end


text:
"""
    function alleles2gt(snp)
---
Merge alleles to genotypes.

- SNP haplotypes are SNP × (Nid × 2)
- SNP alleles are 0 & 1.
- 00 → 0; 01/10 → 1; 11 → 2
"""
function alleles2gt(snp)
    nlc, nid = size(snp)
    nid ÷= 2
    g = Array{Int8, 2}(undef, nlc, nid)
    @inbounds for i in 1:nid
        g[:, i] = snp[:, 2i-1] + snp[:, 2i]
    end
    return g
end

"""
    function hap2gt(hap; GBLUP = false)
Note `hap` here are ID majored.  The two haplotypes of one ID are in two columns.
This is for ease of gene-dropping.
"""
function hap2gt(hap; GBLUP = false)
    nlc, nid = size(hap)
    nid ÷= 2
    if GBLUP
        gt = Matrix{Int8}(undef, nlc, nid)
        for i in 1:nid
            gt[:, i] = hap[:, 2i-1] + hap[:, 2i]
        end
    else
        gt = Matrix{Int8}(undef, nid, nlc)
        for i in 1:nid
            gt[i, :] = hap[:, 2i-1] + hap[:, 2i]
        end
    end
    gt
end

"""
    function breeding_value(snp, qtl)
Given SNP genotypes `snp`, and QTL loci and effects in `qtl`, this function return
true breeding values of each individual.
"""
function breeding_value(snp, qtl)
    l = qtl.pos
    e = qtl.effect             # shorthands
    n = size(snp)[2] ÷ 2       # number of individuals
    g = alleles2gt(snp[l, :])
    return g'e
end

"""
    function phenotype(bv, h²)
---
Given breeding values `bv` and `h²`, and assuming `σₐ²` = 1,
this function returns a phenotype array.
"""
function phenotype(bv, h²)
    σₑ = √(1/h² - 1)            # σₐ² = 1
    ni = length(bv)             # number of individuals
    er = randn(ni) .* σₑ
    return bv, bv + er
end

"""
    function phenotype(snp, qtl, h²)
---
This function returns phenotypes of a trait.
`σₑ²` is determined assuming `σₐ²` = 1 in the base population.
"""
function phenotype(snp, qtl, h²)
    σₑ = √(1/h² - 1)            # σₐ² = 1
    bv = breeding_value(snp, qtl)
    ni = size(bv)
    er = randn(ni) .* σₑ
    return bv, bv + er
end

"""
    function phenotype(snp, qtl, h², threshold)
---
This function returns phenotypes of a binary trait.
`σₑ²` is determined assuming `σₐ² = 1` in the base population.
This result can also serve as a challenge result.
Assuming the liability is for resisitance.
ID with a `phenotype < qtl.mean + threshold` are incident/dead.
That is an ID is alive and phenotype is 1.
"""
function phenotype(snp, qtl, h², percentage)
    bv, ph = phenotype(snp, qtl, h²)
    th = quantile(ph, percentage)
    ni = length(ph)
    bn = zeros(Int8, ni)
    @inbounds for i in 1:ni
        (ph[i] > th) && (bn[i] = 1)
    end
    return bv, bn
end


text:
@testset "Axes" begin

    using Dates
    using YAXArrays.Cubes.Axes:
        axVal2Index, axVal2Index_lb, axVal2Index_ub, findAxis, axname, axsym, renameaxis
    using YAXArrayBase: dimname, dimvals, iscontdim

    axestotest = [
        (CategoricalAxis, "CatAxis", ["One", "Two", "Three"], nothing),
        (RangeAxis, "IntRange", 1:10, 1),
        (RangeAxis, "FloatRange", 1.0:-0.1:0.1, 0.1),
        (
            RangeAxis,
            "MonthTimeRange",
            Date(2001, 1, 1):Month(1):Date(2002, 12, 31),
            Day(30),
        ),
        (
            RangeAxis,
            "DayTimeRange",
            DateTime(2002, 1, 1):Day(1):DateTime(2002, 1, 31),
            Hour(24),
        ),
        (RangeAxis, "PureArray", [0.0, 0.2, 0.4, 0.6], 0.2),
    ]

    for (axt, axn, axv, axstep) in axestotest
        ax = axt(axn, axv)
        @test size(ax) == size(axv)
        @test size(ax, 1) == size(axv, 1)
        @test ndims(ax) == 1
        @test length(ax) == length(axv)
        ax2 = YAXArrays.Cubes.Axes.axcopy(ax, axv)
        @test typeof(ax2) == typeof(ax)
        @test ax2.values == ax.values
        b = IOBuffer()
        show(b, ax)
        @test caxes(ax) == [ax]
        @test YAXArrays.Cubes.Axes.axname(ax) == axn
        @test YAXArrays.Cubes.axsym(ax) == Symbol(axn)
        @test dimname(ax, 1) == axn
        @test dimvals(ax, 1) == axv
        #Test axVal2Index separately
        if ax isa CategoricalAxis
            for (i, v) in enumerate(axv)
                @test axVal2Index(ax, v) == i
                @test axVal2Index(ax, v[1:2], fuzzy = true) == i
            end
            @test iscontdim(ax, 1) == false
        else
            for (i, v) in enumerate(axv)
                @test axVal2Index(ax, v) == i
                @test axVal2Index_ub(ax, v + axstep / 2, fuzzy = true) == i
                @test axVal2Index_lb(ax, v - axstep / 2, fuzzy = true) == i
            end
            @test iscontdim(ax, 1) == true
        end
        ax3 = renameaxis(ax, "Test")
        ax4 = renameaxis(ax, :Test)
        @test axname(ax3) == "Test"
        @test axsym(ax3) == :Test
    end
    axlist = map(i -> i[1](i[2], i[3]), axestotest)
    @test findAxis("Int", axlist) == 2
    @test findAxis(RangeAxis("FloatRange", 1.0:-0.1:0.1), axlist) == 3
    @test getAxis(RangeAxis("FloatRange", 1.0:-0.1:0.1), axlist) ==
          RangeAxis("FloatRange", 1.0:-0.1:0.1)
end


text:
using Tinkering
using Test

@testset "Tinkering.jl" begin
    # Write your tests here.
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
    inflate_ir!(ci::CodeInfo, linfo::MethodInstance) -> ir::IRCode
    inflate_ir!(ci::CodeInfo, sptypes::Vector{Any}, argtypes::Vector{Any}) -> ir::IRCode

Inflates `ci::CodeInfo`-IR to `ir::IRCode`-format.
This should be used with caution as it is a in-place transformation where the fields of
the original `ci::CodeInfo` are modified.
"""
function inflate_ir!(ci::CodeInfo, linfo::MethodInstance)
    sptypes = sptypes_from_meth_instance(linfo)
    if ci.inferred
        argtypes, _ = matching_cache_argtypes(linfo, nothing)
    else
        argtypes = Any[ Any for i = 1:length(ci.slotflags) ]
    end
    return inflate_ir!(ci, sptypes, argtypes)
end
function inflate_ir!(ci::CodeInfo, sptypes::Vector{Any}, argtypes::Vector{Any})
    code = ci.code
    cfg = compute_basic_blocks(code)
    for i = 1:length(code)
        stmt = code[i]
        # Translate statement edges to bb_edges
        if isa(stmt, GotoNode)
            code[i] = GotoNode(block_for_inst(cfg, stmt.label))
        elseif isa(stmt, GotoIfNot)
            code[i] = GotoIfNot(stmt.cond, block_for_inst(cfg, stmt.dest))
        elseif isa(stmt, PhiNode)
            code[i] = PhiNode(Int32[block_for_inst(cfg, Int(edge)) for edge in stmt.edges], stmt.values)
        elseif isexpr(stmt, :enter)
            stmt.args[1] = block_for_inst(cfg, stmt.args[1]::Int)
            code[i] = stmt
        end
    end
    nstmts = length(code)
    ssavaluetypes = ci.ssavaluetypes
    if !isa(ssavaluetypes, Vector{Any})
        ssavaluetypes = Any[ Any for i = 1:ssavaluetypes::Int ]
    end
    info = Any[nothing for i = 1:nstmts]
    stmts = InstructionStream(code, ssavaluetypes, info, ci.codelocs, ci.ssaflags)
    linetable = ci.linetable
    if !isa(linetable, Vector{LineInfoNode})
        linetable = collect(LineInfoNode, linetable::Vector{Any})::Vector{LineInfoNode}
    end
    meta = Expr[]
    return IRCode(stmts, cfg, linetable, argtypes, meta, sptypes)
end

"""
    inflate_ir(ci::CodeInfo, linfo::MethodInstance) -> ir::IRCode
    inflate_ir(ci::CodeInfo, sptypes::Vector{Any}, argtypes::Vector{Any}) -> ir::IRCode
    inflate_ir(ci::CodeInfo) -> ir::IRCode

Non-destructive version of `inflate_ir!`.
Mainly used for testing or interactive use.
"""
inflate_ir(ci::CodeInfo, linfo::MethodInstance) = inflate_ir!(copy(ci), linfo)
inflate_ir(ci::CodeInfo, sptypes::Vector{Any}, argtypes::Vector{Any}) = inflate_ir!(copy(ci), sptypes, argtypes)
inflate_ir(ci::CodeInfo) = inflate_ir(ci, Any[], Any[ Any for i = 1:length(ci.slotflags) ])

function replace_code_newstyle!(ci::CodeInfo, ir::IRCode, nargs::Int)
    @assert isempty(ir.new_nodes)
    # All but the first `nargs` slots will now be unused
    resize!(ci.slotflags, nargs)
    stmts = ir.stmts
    code = ci.code = stmts.inst
    ssavaluetypes = ci.ssavaluetypes = stmts.type
    codelocs = ci.codelocs = stmts.line
    ssaflags = ci.ssaflags = stmts.flag
    linetable = ci.linetable = ir.linetable
    for metanode in ir.meta
        push!(code, metanode)
        push!(codelocs, 1)
        push!(ssavaluetypes, Any)
        push!(ssaflags, IR_FLAG_NULL)
    end
    # Translate BB Edges to statement edges
    # (and undo normalization for now)
    for i = 1:length(code)
        stmt = code[i]
        if isa(stmt, GotoNode)
            code[i] = GotoNode(first(ir.cfg.blocks[stmt.label].stmts))
        elseif isa(stmt, GotoIfNot)
            code[i] = GotoIfNot(stmt.cond, first(ir.cfg.blocks[stmt.dest].stmts))
        elseif isa(stmt, PhiNode)
            code[i] = PhiNode(Int32[last(ir.cfg.blocks[edge].stmts) for edge in stmt.edges], stmt.values)
        elseif isexpr(stmt, :enter)
            stmt.args[1] = first(ir.cfg.blocks[stmt.args[1]::Int].stmts)
            code[i] = stmt
        end
    end
end

max_stars_repo_name: Physolia/julia

max_issues_repo_name: Eurus-Holmes/julia

max_forks_repo_name: Nexuscompute/julia


text:
# Support for CUDA

if Pkg.installed("CUDArt") == nothing
    macro cuda(ex)
    end
else
    using CUDArt

    macro cuda(ex)
        ex
    end

    immutable CudaDevice <: Device end

    export CudaDevice

end


text:
workspace()
using TimeWindowStats
using Base.Test, Dates
reload("TimeWindowStats")
dt1 = DateTime(2015,6,4,7,1,1)
dt2 = dt1 + Second(2)
dt3 = dt1 + Second(4)
dt4 = dt1 + Second(5)
dt5 = dt1 + Second(8)

c = Sum(5000,4)
update!(c,1., dt1)
@test c.head == 1
@test c.tail == 1
@test c.len == 1
@test c.v == 1.
@test typeof(c.xs) == Vector{Float64}
@test typeof(c.v) == Float64

update!(c,2., dt2)
@test c.head == 2
@test c.tail == 1
@test c.len == 2
@test c.v == 3.

update!(c,3., dt3)
@test c.head == 3
@test c.tail == 1
@test c.len == 3
@test c.v == 6.

update!(c,4., dt4)
@test c.head == 4
@test c.tail == 1
@test c.len == 4
@test c.size == 5
@test c.v == 10.

update!(c,5., dt4)
@test c.head == 5
@test c.tail == 1
@test c.len == 5
@test c.size == 6
@test c.v == 15.

c.head
c.xs

update!(c,6., dt5)
@test c.head == 6
@test c.tail == 3
@test c.len == 4
@test c.size == 6
@test c.v == 18.

update!(c,7., dt5)
@test c.head == 1
@test c.tail == 3
@test c.len == 5
@test c.size == 6
@test c.v == 25.

update!(c,8., dt5)
@test c.head == 6
@test c.tail == 1
@test c.len == 6
@test c.size == 7
@test c.v == 33.

update!(c,9., dt5)
@test c.head == 7
@test c.tail == 1
@test c.len == 7
@test c.size == 8
@test c.v == 42.

c = Sum(5000)
update!(c,1, dt1)
@test c.v == 1
@test typeof(c.xs) == Vector{Int64}
@test typeof(c.v) == Int64

c = Pct(5000)
update!(c,1,1, dt1)
@test c.v == 1.
update!(c,0,1, dt2)
@test c.v == .5

c = Mean(5000)
update!(c,1, dt1)
@test c.v == 1.
update!(c,2, dt2)
@test c.v == 1.5
update!(c,3, dt3)
@test c.v == 2.
update!(c,4, dt4)
@test c.v == 2.5
update!(c,5, dt5)
@test c.v == 4.

max_stars_repo_name: Andy-P/StreamingUtil

max_issues_repo_name: Andy-P/StreamingUtil

max_forks_repo_name: Andy-P/StreamingUtil


text:
module TestLogParsing
include("../../../examples/logparsing.jl")
end  # module


text:
# advent of code 2020 day 3

function load_data(path:: String):: AbstractArray
    file = open(path)
    data_array = readlines(file)
    close(file)
    return data_array
end

function walker_acc(data_array:: AbstractArray, x_step=1, y_step=3):: Int
    size_x, size_y = size(data_array)[1], length(data_array[1])
    numb_tree, pos_x, pos_y = 0, 1 + x_step, 1 + y_step
    while (pos_x <= size_x)
        if data_array[pos_x][(pos_y - 1)%end + 1] == '#'
            numb_tree += 1
        end
        pos_x += x_step
        pos_y += y_step
    end

    return numb_tree
end

data_array = load_data("./data/day3.txt")
numb_tree = walker_acc(data_array)
println("Numb of tree is: $numb_tree")

slope_to_test = [[1, 1], [3, 1], [5, 1], [7, 1], [1, 2]]
let prod_path = 1
    for slope in slope_to_test
        prod_path *= walker_acc(data_array, slope[2], slope[1])
    end
    println("product of multipel slopes is: $prod_path")
end