text:
using CurrencyAmounts
using Statistics
using Test

@currencies USD, EUR, GBP, CHF

@testset "CurrencyAmounts.jl" begin

    @test EUR == Currency{:EUR}()

    @test 4USD == CurrencyAmounts.CurrencyAmount{Int64,Currency{:USD}}(4)
    @test 4USD + 5.5USD == 9.5USD
    @test 4USD - 5.5USD == -1.5USD
    @test 4.5USD * 2 == 9.0USD
    @test 4*4.5USD * 2 == 36.0USD
    @test 8USD / 2 == 4.0USD
    @test 8USD / 2 ≈ 4.0USD
    @test 6EUR / 3EUR == 2.0
    @test 4USD/USD == 4.0
    @test -4USD == CurrencyAmounts.CurrencyAmount{Int64,Currency{:USD}}(-4)
    @test zero(3EUR) == 0EUR
    @test zero(typeof(3.5EUR)) == 0.0EUR
    @test 3EUR < 3.1EUR
    @test 6.8USD > 6.7USD
    @test_throws MethodError 2EUR < 5USD
    @test isfinite(3USD)
    @test !isfinite(Inf * USD)

    @test Statistics.middle(4EUR, 5EUR) == 4.5EUR
    @test median([2EUR, 4EUR, 5EUR, 8EUR]) == 4.5EUR
    @test isnan(4EUR) == false
    @test isnan(NaN*EUR) == true

    # invalid operations
    @test_throws MethodError 2 / 8USD
    @test_throws MethodError 3USD * 4USD
    @test_throws MethodError 4USD + 3EUR
    @test_throws MethodError 4USD - 3EUR

    # exchange rates
    rate_eur_usd = 1.2USD / 1EUR
    @test rate_eur_usd == ExchangeRate{Float64,Currency{:EUR},Currency{:USD}}(1.2) == ExchangeRate(1.2, EUR, USD)
    @test 2.4USD / 2EUR == rate_eur_usd
    @test rate_eur_usd * 2 == 2.4USD / 1EUR
    @test rate_eur_usd/2 == ExchangeRate{Float64,Currency{:EUR},Currency{:USD}}(0.6)
    @test 1/rate_eur_usd == ExchangeRate{Float64,Currency{:USD},Currency{:EUR}}(1/1.2)
    @test rate_eur_usd + 2*rate_eur_usd == 3*rate_eur_usd
    @test ExchangeRate(2.2, EUR, USD) - rate_eur_usd ≈ ExchangeRate(1.0, EUR, USD)
    @test_throws ErrorException ExchangeRate(1, EUR, EUR)

    @test rate_eur_usd / rate_eur_usd == 1
    @test (2.4USD / EUR) / rate_eur_usd == 2

    # triangulation
    @test (1.2USD/EUR) / (0.8USD/GBP) ≈ (1.2/0.8)*GBP/EUR
    @test (1.2USD/EUR) / (1.1GBP/EUR) ≈ (1.2/1.1)*USD/GBP
    @test (1.2USD/EUR) * (2.2EUR/GBP) ≈ (1.2*2.2)*USD/GBP
    @test (1.2USD/EUR) * (1.1GBP/USD) ≈ (1.2*1.1)*GBP/EUR

    @test Statistics.middle(rate_eur_usd, 2*rate_eur_usd) == 1.5rate_eur_usd
    @test isnan(rate_eur_usd) == false

    # currency conversions
    @test 3EUR * rate_eur_usd ≈ 3.6USD
    @test rate_eur_usd * 3EUR == 3EUR * rate_eur_usd
    @test_throws MethodError 3USD * rate_eur_usd # dimensions do not match
    @test 3USD / rate_eur_usd == 2.5EUR

    # test pretty printing
    io = IOBuffer(append=true)
    print(io, EUR)
    @test read(io, String) == "EUR"
    print(io, 100.5EUR)
    @test read(io, String) == "100.5 EUR"
    print(io, 1.2USD/EUR)
    @test read(io, String) == "1.2 USD/EUR"

    # conversions
    @test convert(EUR, 2.4USD, 1.2USD/EUR) == 2.0EUR
    @test convert(USD, 2EUR, 1.2USD/EUR) == 2.4USD
    @test convert(USD, 4USD, 1.2USD/EUR) == 4USD
    @test convert(GBP, 4USD, 1.2USD/EUR) === missing

    exchange_rates = (1.2USD/EUR, 1.3EUR/GBP, 1.56USD/GBP)
    @test convert(GBP, 5.2EUR, exchange_rates) == 4.0GBP
    @test convert(CHF, 5.2USD, exchange_rates) === missing # this exchange rate is not given, no automatic triangulation performed

end


text:
using OmegaCore
using Distributions

x = 1 ~ Normal(0, 1)
y = 2 ~ Normal(0, 3)

function f(id, ω)
  a = ([id; 3] ~ Normal(0, 1))(ω)
  for i = 1:10
    a += ([id; i] ~ Normal(0, 1))(ω)
  end
  a + x(ω) + y(ω)
end

f1 = 3 ~ f
f2 = 4 ~ f

# Need to make sure that all ids are globally unique
# Need to compose id with everything

ext: jl

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
function output_data(name,funcs,pvals,sol,tsol)
    Nend = sol.t[end]
    bNs = 0.0:0.01:Nend
    d = length(sol.u[1])÷2
    #Mσσ(N) = funcs["Mσσ"]([sol(N);pvals]...)
    #Mss(N) = funcs["Mss"]([sol(N);pvals]...)
    Eh(N) = funcs["Eh"]([sol(N)[1:2*d];pvals]...)
    Ev(N) = funcs["Ev"](sol(N)...,pvals...)
    Om(N) = funcs["Om"]([sol(N)[1:2*d];pvals]...)
    if haskey(funcs,"Mab")
        Mab(N) = funcs["Mab"]([sol(N)[1:2*d];pvals]...)
        mabs = Mab.(bNs)
    end
    #Mσσs = Mσσ.(bNs)
    #Msss = Mss.(bNs)
    ehs = Eh.(bNs)
    evs = Ev.(bNs)
    oms = Om.(bNs)
    phis = [u[1:d] for u in sol.(bNs)]
    pis = [u[d+1:end] for u in sol.(bNs)]
    jldopen(name,"w") do f
        write(f,"eh",ehs)
        write(f,"ev",evs)
        write(f,"omega",oms)
        if haskey(funcs,"Mab")
            write(f,"Mab",mabs)
        end
        write(f,"Phi",phis)
        write(f,"Pi",pis)
        write(f,"pvalues",pvals)
        write(f,"back_Ne",bNs)
        if !isnothing(tsol)
            for (key,val) in tsol
                write(f,key,val)
            end
        end
    end
end

function generate_default_plots(output_fname)
    # write to a pluto notebook
end


text:
# adaptors for converting abstract arrays to Base.Array

## arrays

Adapt.adapt_storage(::Type{Array}, xs::AT) where {AT<:AbstractArray} =
  convert(Array, xs)

# if an element type is specified, convert to it
Adapt.adapt_storage(::Type{<:Array{T}}, xs::AT) where {T, AT<:AbstractArray} =
  convert(Array{T}, xs)

# NOTE: this flattens all <:AbstractArray leaves, e.g., Base.Slice(1:2) -> [1,2].


text:
function filter_cells!(dataset::AbstractVTKUnstructuredData, cell_types_to_remove::Vector{Int})
    inds = findall(i -> i in cell_types_to_remove, dataset.cell_types)
    deleteat!(dataset.cell_types, inds)
    deleteat!(dataset.cell_connectivity, inds)
    for m in keys(dataset.cell_data)
        _var_dim = var_dim(dataset, m, "Cell")
        if _var_dim == 1
            deleteat!(dataset.cell_data[m], inds)
        else
            dataset.cell_data[m] = dataset.cell_data[m][:,setdiff(1:num_of_cells(dataset), inds)]
        end
    end
    return dataset
end

function keep_volume_cells_only!(dataset::VTKUnstructuredData)
    return filter_cells!(dataset, [POINT_CELLS; LINE_CELLS; FACE_CELLS])
end
function keep_face_cells_only!(dataset::VTKPolyData)
    return filter_cells!(dataset, [POINT_CELLS; LINE_CELLS])
end
function add_new_cell!( dataset::AbstractVTKUnstructuredData,
                        _cell_type::Int,
                        _cell_connectivity,
                        _filter::Bool = true
                      )
    TData = typeof(dataset)
    begin
        out = true
        for i in 1:length(dataset.cell_connectivity)
            (out = !similar_cells(dataset.cell_connectivity[i], cell_type(dataset, i),
                _cell_connectivity, _cell_type)) || break
        end
        out
    end || _filter && (return num_of_cells(dataset)) || throw("$(TData): Repeat cells are not allowed.")

    push!(dataset.cell_connectivity, _cell_connectivity)
    push!(dataset.cell_types, _cell_type)
    for m in keys(dataset.cell_data)
        _var_dim = var_dim(dataset, m, "Cell")
        if _var_dim == 1
            dataset.cell_data[m] = [dataset.cell_data[m]; [0]]
        else
            dataset.cell_data[m] = [dataset.cell_data[m] zeros(size(dataset.cell_data[m], 1))]
        end
    end
    return num_of_cells(dataset)
end

function remove_cell!(dataset::AbstractVTKUnstructuredData, cell_ind::Int)
    deleteat!(dataset.cell_connectivity, cell_ind)
    deleteat!(dataset.cell_types, cell_ind)
    for m in keys(dataset.cell_data)
        _var_dim = var_dim(dataset, m, "Cell")
        if _var_dim == 1
            dataset.cell_data[m] = vcat(dataset.cell_data[m][1:cell_ind-1],
                                            dataset.cell_data[m][cell_ind+1:end])
        else
            dataset.cell_data[m] = hcat(dataset.cell_data[m][:, 1:cell_ind-1],
                                            dataset.cell_data[m][:, cell_ind+1:end])
        end
    end
    return num_of_cells(dataset)
end

function add_point_id_offset!(cell_connectivity, offset::Int)
    for i in 1:length(cell_connectivity)
        cell_connectivity[i] = cell_connectivity[i] .+ offset
    end
    return cell_connectivity
end

function append(datasets::AbstractVTKUnstructuredData...)
    dataset1 = deepcopy(datasets[1])
    for i in 2:length(datasets)
        append!(dataset1, datasets[i])
    end
    return dataset1
end

insert_new_block!(dataset::AbstractVTKMultiblockData, block) = push!(dataset.blocks, block)
remove_block!(dataset::AbstractVTKMultiblockData, ind::Int) = deleteat!(dataset.blocks, ind)
function remove_block!(dataset::AbstractVTKMultiblockData, block::AbstractVTKMultiblockData)
    return deleteat!(dataset.blocks, findin(dataset.blocks, block))
end
function insert_timed_data!(dataset::VTKTimeSeriesData{TTime, TData},
                            time::TTime,
                            data::TData
                           ) where {TTime<:Real, TData<:AbstractStaticVTKData}
    loc = searchsorted(dataset.timemarkers, time)
    if loc.start == length(dataset) + 1 && loc.stop == length(dataset)
        push!(dataset.timemarkers, time)
        push!(dataset.data, data)
    elseif loc.start == 1 && loc.stop == 0
        unshift!(dataset.timemarkers, time)
        unshift!(dataset.data, data)
    else
        insert!(dataset.timemarkers, loc.start, time)
        insert!(dataset.data, loc.start, data)
    end
end

function remove_timed_data!(dataset::VTKTimeSeriesData{TTime, TData},
                            _time::TTime
                           ) where {TTime<:Real, TData<:AbstractStaticVTKData}
    ind = findin(dataset.timemarkers, _time)
    if typeof(ind) == Int
        deleteat!(dataset.timemarkers, ind)
        deleteat!(dataset.data, ind)
    else
        throw("No data at time $_time.")
    end
end

function remove_timed_data!(dataset::VTKTimeSeriesData, ind::Int)
    if ind >= 1 && ind <= length(dataset)
        deleteat!(dataset.timemarkers, ind)
        deleteat!(dataset.data, ind)
    else
        throw("Index out of bound!")
    end
end

function increase_resolution!(dataset::VTKTimeSeriesData, resolution::Int)
    _times = deepcopy(collect(zip(dataset.timemarkers[1:end-1], dataset.timemarkers[2:end])))
    for (t1, t2) in _times
        _step = (t2-t1)/resolution
        for j in 1:resolution-1
            VTKDataTypes.insert_timed_data!(dataset, Float64(t1+j*_step), dataset[t1+j*_step])
        end
    end
    return
end

@resumable function simple_block_generator(multiblock::VTKMultiblockData)
    for block in multiblock
        if typeof(block) <: VTKMultiblockData
            for b in simple_block_generator(block)
                @yield b
            end
        else
            @yield block
        end
    end
end

function dim3!(a::AbstractVTKUnstructuredData)
    if dim(a) == 3
        return
    elseif dim(a) == 2
        a.point_coords = [a.point_coords; zeros(1,num_of_points(a))]
        return
    else
        throw("Invalid dimension.")
    end
end

function dim3!(a::VTKUniformRectilinearData)
    if dim(a) == 3
        return
    elseif dim(a) == 2
        a.origin = [a.origin; 0]
        a.spacing = [a.spacing; 0]
        a.extents = [a.extents; 1]
        for m in keys(a.point_data)
            a.point_data[m] = reshape(a.point_data[m], (size(a.point_data[m])..., 1))
        end
        for m in keys(a.cell_data)
            a.cell_data[m] = reshape(a.cell_data[m], (size(a.cell_data[m])..., 1))
        end

        return
    else
        throw("Invalid dimension.")
    end
end

function dim3!(a::VTKRectilinearData)
    if dim(a) == 3
        return
    elseif dim(a) == 2
        a.point_coords = [a.point_coords; [0]]
        for m in keys(a.point_data)
            a.point_data[m] = reshape(a.point_data[m], (size(a.point_data[m])..., 1))
        end
        for m in keys(a.cell_data)
            a.cell_data[m] = reshape(a.cell_data[m], (size(a.cell_data[m])..., 1))
        end

        return
    else
        throw("Invalid dimension.")
    end
end

function dim3!(a::VTKStructuredData)
    if dim(a) == 3
        return
    elseif dim(a) == 2
        a.point_coords = reshape(cat(1, a.point_coords, zeros(1,extents(a)...)), (3, extents(a), 1))
        for m in keys(a.point_data)
            a.point_data[m] = reshape(a.point_data[m], (size(a.point_data[m])..., 1))
        end
        for m in keys(a.cell_data)
            a.cell_data[m] = reshape(a.cell_data[m], (size(a.cell_data[m])..., 1))
        end

        return
    else
        throw("Invalid dimension.")
    end
end

function dim3!(a::VTKMultiblockData)
    if dim(a) == 3
        return
    elseif dim(a) == 2
        for block in simple_block_generator(a)
            dim3!(block)
        end
        return
    else
        throw("Invalid dimension.")
    end
end

function dim3!(a::VTKTimeSeriesData)
    if dim(a) == 3
        return
    elseif dim(a) == 2
        for timedblock in a
            dim3!(timedblock)
        end
        return
    else
        throw("Invalid dimension.")
    end
end

function celldata_to_pointdata!(dataset::AbstractVTKUnstructuredData)
    point_coords = dataset.point_coords
    point_data = dataset.point_data

    for m in keys(dataset.cell_data)
        _var_dim = var_dim(dataset, m, "Cell")
        if _var_dim == 1
            point_data[m] = zeros(num_of_points(dataset))
        else
            point_data[m] = zeros(_var_dim, num_of_points(dataset))
        end
    end

    point_counter = zeros(num_of_points(dataset))
    for i in 1:num_of_cells(dataset), j in dataset.cell_connectivity[i]
        point_counter[j] += 1
        for m in keys(dataset.cell_data)
            _var_dim = var_dim(dataset, m, "Cell")
            if _var_dim == 1
                point_data[m][j] += dataset.cell_data[m][i]
            else
                @views point_data[m][:,j] += dataset.cell_data[m][:,i]
            end
        end
    end

    for m in keys(dataset.cell_data)
        _var_dim = var_dim(dataset, m, "Cell")
        if _var_dim == 1
            point_data[m] = point_data[m] ./ max.(point_counter, 1)
        else
            point_data[m] = point_data[m] ./ max.(point_counter', 1)
        end
    end
    empty!(dataset.cell_data)

    return
end


text:
using SafeTestsets
@safetestset "Calc Tests" begin include("CalcTests.jl") end


text:
export solicitation_table_parser

function solicitation_table_parser(df,target_rho,total_vol)

    lbbbl_to_kgm3=2.857142857142857u"kg/m^3"

    ppg_to_kgm3=119.83u"kg/m^3"

    galbbl_to_mLmL=0.0238

    mass=Float64[]
    mass_unit=Unitful.Quantity[]

    QSP_indexes=Integer[]
    vol_indexes=Integer[]

    water_index=Integer[]
    water_content=Fluid[]
    oil_index=Integer[]
    oil_content=Fluid[]
    solids_index=Integer[]
    solid_content=Solid[]

    aditive_indexes=Integer[]

    has_QSP=false
    has_pure_water=false
    pure_water_index=0
    len_components=size(df)[1]

    #check if QSP of water and salt
    comp_df=copy(df)
    if any(df.Componente .== :H2O) && any(df.Componente .== :NaCl)
        w_i=findall(df.Componente .== :H2O)
        Na_i=findall(df.Componente .== :NaCl)

        w_QSP=any()

    end

    for i in range(1,len_components)

        row=df[i,:]

        # println(row.Componente)
        if row.unit == "lb/bbl"
            convert=lbbbl_to_kgm3
            m=upreferred(row.conc*convert*total_vol)
            if row.Componente == :Obiturante
                append!(solid_content,[CaCO3()])
                append!(solids_index,i)
            elseif row.Componente == :Reologico
                ben=Bentonite()
                append!(solid_content,[ben])
                append!(solids_index,i)
            end
        elseif row.unit == "bbl/bbl" || row.unit == "v/v"
            vol=row.conc*total_vol
            append!(vol_indexes,i)
            rho=1.0
            if row.Componente == :Olefina || row.Componente == :Hexadeceno
                fluid=Hexadecene()
                rho=0.728*1000.0
                # rho=fluid.ρ
                append!(oil_index,i)
                append!(oil_content,[fluid])
            elseif row.Componente == :Brine_NaCl && !ismissing(row.conc)
                fluid=BrineNaCl(9.8)
                rho=fluid.ρ
                append!(water_index,i)
                append!(water_content,[fluid])
            elseif row.Componente == :Obiturante
                append!(solid_content,[CaCO3()])
                rho=CaCO3().ρ
                append!(solids_index,i)
            end

            m=upreferred(vol*rho*u"kg/m^3")
        elseif row.unit == "% v/v"
            vol=row.conc*total_vol/100

            if row.Componente == :Liovac4260
                rho=700.0
            end
            println("Liovac")
            m=upreferred(vol*rho*u"kg/m^3")
        elseif row.unit == "gal/bbl"
            vol=galbbl_to_mLmL*total_vol

            if row.Componente == :Polifoan
                rho=1.01u"g/cm^3"
            end

            m=upreferred(vol*rho)
        elseif ismissing(row.conc) && occursin(row.unit,"QSP")
            has_QSP=true
            append!(QSP_indexes,i)
            m=0.0u"kg"

            if row.Componente == :Barita
                append!(solids_index,i)
                append!(solid_content,[Barite()])
            elseif row.Componente == :H2O
                has_pure_water=true
                pure_water_index=i
            end

        end
        println((i=i,))
        append!(mass_unit,m)
        append!(mass,ustrip(m))
    end
    println(mass_unit)

    ## Encontra as massas dos QSP...
    fluid_rho=target_rho*ppg_to_kgm3
    total_mass=upreferred(fluid_rho*total_vol)
    if !isempty(QSP_indexes)
        if length(QSP_indexes) == 1
            m_a=0.0u"kg"
            m_B=0.0u"kg"
            for i in range(1,len_components)
                if df[i,:].aditive
                    append!(aditive_indexes,i)
                    m_a+=mass_unit[i]
                end
                if i in vol_indexes
                    m_B+=mass_unit[i]
                end
            end
            m_s=0.0u"kg"
            for j in solids_index
                if j ∉ QSP_indexes
                    m_s+=mass_unit[j]
                end
            end
            mass_unit[QSP_indexes[1]]=total_mass-m_a-m_B
        else
            # if has_pure_water && any(findall(.Componente .== :NaCl))
            #     println("hello")
            #     @error "Not Implemented"
            # end
        end
    end

    ## Encontrar as frações volumétricas considerando aditivos dissolvidos na água.

    fP=[0.0,0.0,0.0]
    fcomp=[
        [],
        [],
        []
    ]

    vol_sol=Unitful.Quantity[]
    vol_oil=0.0u"m^3"
    vol_water=0.0u"m^3"

    for (i,f) in zip(oil_index,oil_content)
        rho=f.ρ*u"kg/m^3"
        vol=upreferred(mass_unit[i]/rho)
        vol_oil+=vol
        append!(fcomp[2],vol_oil)
    end

    for (i,s) in zip(solids_index,solid_content)
        rho=s.ρ*u"kg/m^3"
        vol=upreferred(mass_unit[i]/rho)
        append!(vol_sol,vol)
        append!(fcomp[3],vol)
    end

    vol_water=total_vol-sum(vol_oil)-sum(vol_sol)

    fP[1]=upreferred(vol_water/total_vol)
    fP[2]=upreferred(sum(vol_oil)/total_vol)
    fP[3]=upreferred(sum(vol_sol)/total_vol)

    fcomp[1] = [upreferred(vol_water/total_vol)]
    fcomp[2] .= upreferred.(fcomp[2]./total_vol)
    fcomp[3] = upreferred.(fcomp[3]./total_vol)

    mass_frac=mass_unit./total_mass

    df.mass=mass_unit
    (
        total_mass=total_mass,
        total_mass_no_unit=ustrip(total_mass),
        total_vol=total_vol,
        total_vol_no_unit=ustrip(upreferred(total_vol)),
        components_mass=mass_unit,
        components_mass_no_unit=mass,
        aditive_mass=m_a,
        mass_frac=mass_frac,
        fP=fP,
        fcomp=fcomp,
        rho_std=target_rho*ppg_to_kgm3,
        solid_indexes=solids_index,
        aditive_indexes=aditive_indexes,
        df=df,
        # sol_indexes=solids_index
    )

end


text:
using ElasticAnisotropy
using Test

O = OrthoModel(8e9,4e9, 0.35, 0.5, 1.48148e9)
n = [-0.573576436, 0.0, 0.819152044]
R = getRmatrix(n)
C_voigt = generateVoigtMatrix(O)
c_tensor = VoigtMatrix2Tensor(C_voigt)
c_tensor_R = rotateTensor(c_tensor,n)
c_tensor_R2 = rotateTensor(c_tensor,R)
C_voigt_R = Tensor2VoigtMatrix(c_tensor_R)
C_voigt_R2 = Tensor2VoigtMatrix(c_tensor_R2)

@testset "ElasticAnisotropy.jl" begin
    @test C_voigt_R == C_voigt_R2
end


text:
using Test, Statistics

pdir = joinpath(@__DIR__, "../examples/m12.6sl")
cd(pdir) do

  isdir("tmp") && rm("tmp", recursive=true);

  include(joinpath(pdir, "m12.6sl.jl"))

  @test 1.0 <  mean(m12_6sl[:alpha].value)< 1.2

  isdir("tmp") && rm("tmp", recursive=true);

end


text:
actions = getallactions!(client)

@testset "List all Actions" begin
    @test actions[1].id == 36804636
    @test actions[1].status == "completed"
    @test actions[1].actiontype == "create"
    @test actions[1].started_at == DateTime("2014-11-14T16:29:21")
    @test actions[1].completed_at == DateTime("2014-11-14T16:30:06")
    @test actions[1].resourceid == 3164444
    @test actions[1].resourcetype == "droplet"
    @test actions[1].region.slug == "nyc3"
    @test actions[1].regionslug == "nyc3"
    @test "s-16vcpu-64gb" in actions[1].region.sizes
    @test "private_networking" in actions[1].region.features
end;

action = getaction!(client, 36804636)

@testset "Retrieve an existing Action" begin
    @test action.id == 36804636
    @test action.status == "completed"
    @test action.actiontype == "create"
    @test action.started_at == DateTime("2014-11-14T16:29:21")
    @test action.completed_at == DateTime("2014-11-14T16:30:06")
    @test action.resourceid == 3164444
    @test action.resourcetype == "droplet"
    @test action.region.slug == "nyc3"
    @test action.regionslug == "nyc3"
    @test "s-4vcpu-8gb" in action.region.sizes
    @test "metadata" in action.region.features
end;


text:
#Scales
#xlim!/ylim!
for fun in [:x, :y]

    fname = Symbol("$(fun)lim!")

    @eval begin function ($fname)(v::VegaVisualization; min::Real = 0, max::Real = 1000, reverse::Bool = false, round::Bool = false, _type::AbstractString = "", exponent::Real = 0.3)

                    if min > max
                        error("Min must be less than Max")
                    end

                    s = v.scales[findfirst([z.name == $"$fun" for z in v.scales])]
                    s.domainMin = min == 0? s.domainMin : min
                    s.domainMax = max == 1000? s.domainMax : max
                    s.reverse = reverse
                    s.round = round
                    min != 0? s.zero = false : s.zero = true

                    #Allow for changing scale type; if type == pow, also need to set exponent. Defaulting to 0.3 instead of d3 value of 1, since 1 returns identity
                    #Returning identity is safer, but also makes user believe function is broken
                    _type != ""? s._type = _type: nothing
                    _type == "pow"? s.exponent = exponent: nothing
                    return v
                end
          end
end

max_stars_repo_name: jpthor/vega


text:
# ===============================
# Written by AAE
# TU Delft, Spring 2014
# simulkade.com
# ===============================

# ===============================================================
# Changes
#    2014-12-29 Import all the convection terms from matlab code
#    2015-01-10 extended to accept nonuniform grids
# ===============================================================

function convectionTerm(u::FaceValue)
d = u.domain.dimension
if d==1
  M = convectionTerm1D(u)
elseif d==1.5
  M = convectionTermCylindrical1D(u)
elseif d==2
  M, Mx, My = convectionTerm2D(u)
elseif d==2.5
  M, Mx, My = convectionTermCylindrical2D(u)
elseif d==2.8
  M, Mx, My = convectionTermRadial2D(u)
elseif d==3
  M, Mx, My, Mz = convectionTerm3D(u)
elseif d==3.2
  M, Mx, My, Mz = convectionTermCylindrical3D(u)
end
M
end

function convectionUpwindTerm(u::FaceValue)
d = u.domain.dimension
if d==1
  M = convectionUpwindTerm1D(u)
elseif d==1.5
  M = convectionUpwindTermCylindrical1D(u)
elseif d==2
  M, Mx, My = convectionUpwindTerm2D(u)
elseif d==2.5
  M, Mx, My = convectionUpwindTermCylindrical2D(u)
elseif d==2.8
  M, Mx, My = convectionUpwindTermRadial2D(u)
elseif d==3
  M, Mx, My, Mz = convectionUpwindTerm3D(u)
elseif d==3.2
  M, Mx, My, Mz = convectionUpwindTermCylindrical3D(u)
end
M
end

function convectionTvdTerm(u::FaceValue, phi::CellValue, FL::Function)
d = u.domain.dimension
if d==1
  M, RHS = convectionTvdTerm1D(u::FaceValue, phi::CellValue, FL::Function)
elseif d==1.5
  M, RHS = convectionTvdTermCylindrical1D(u::FaceValue, phi::CellValue, FL::Function)
elseif d==2
  M, RHS, Mx, My, RHSx, RHSy = convectionTvdTerm2D(u::FaceValue, phi::CellValue, FL::Function)
elseif d==2.5
  M, RHS, Mx, My, RHSx, RHSy = convectionTvdTermCylindrical2D(u::FaceValue, phi::CellValue, FL::Function)
elseif d==2.8
  M, RHS, Mx, My, RHSx, RHSy = convectionTvdTermRadial2D(u::FaceValue, phi::CellValue, FL::Function)
elseif d==3
  M, RHS, Mx, My, Mz, RHSx, RHSy, RHSz = convectionTvdTerm3D(u::FaceValue, phi::CellValue, FL::Function)
elseif d==3.2
  M, RHS, Mx, My, Mz, RHSx, RHSy, RHSz = convectionTvdTermCylindrical3D(u::FaceValue, phi::CellValue, FL::Function)
end
(M, RHS)
end

# =================== 1D Convection Terms =======================
# =============== 1D Convection Terms Cartesian Central =================
function convectionTerm1D(u::FaceValue)
# u is a face variable

# extract data from the mesh structure
Nx = u.domain.dims[1]
G = [1:Nx+2;]
#DX = u.domain.cellsize.x

DXe = u.domain.cellsize.x[3:end]
DXw = u.domain.cellsize.x[1:end-2]
DXp = u.domain.cellsize.x[2:end-1]

# define the vectors to store the sparse matrix data
iix = zeros(Int64, 3*(Nx+2))
jjx = zeros(Int64, 3*(Nx+2))
sx = zeros(Float64, 3*(Nx+2))

# reassign the east, west for code readability
ue = u.xvalue[2:Nx+1]./(DXp+DXe)
uw = u.xvalue[1:Nx]./(DXp+DXw)

# calculate the coefficients for the internal cells
AE = reshape(ue,Nx)
AW = reshape(-uw,Nx)
APx = reshape((ue.*DXe-uw.*DXw)./DXp,Nx)

# build the sparse matrix based on the numbering system
rowx_index = reshape(G[2:Nx+1],Nx) # main diagonal x
iix[1:3*Nx] = repmat(rowx_index,3)
jjx[1:3*Nx] = [reshape(G[1:Nx],Nx); reshape(G[2:Nx+1],Nx); reshape(G[3:Nx+2],Nx)]
sx[1:3*Nx] = [AW; APx; AE]

# build the sparse matrix
kx = 3*Nx
M = sparse(iix[1:kx], jjx[1:kx], sx[1:kx], Nx+2, Nx+2)

end

# =============== 1D Convection Terms Cylindrical Central =================
function convectionTermCylindrical1D(u::FaceValue)
# u is a face variable

# extract data from the mesh structure
Nx = u.domain.dims[1]
G = [1:Nx+2;]
DXe = u.domain.cellsize.x[3:end]
DXw = u.domain.cellsize.x[1:end-2]
DXp = u.domain.cellsize.x[2:end-1]
rp = u.domain.cellcenters.x
rf = u.domain.facecenters.x

# define the vectors to store the sparse matrix data
iix = zeros(Int64, 3*(Nx+2))
jjx = zeros(Int64, 3*(Nx+2))
sx = zeros(Float64, 3*(Nx+2))

# reassign the east, west for code readability
ue = rf[2:Nx+1].*u.xvalue[2:Nx+1]./(rp.*(DXp+DXe))
uw = rf[1:Nx].*u.xvalue[1:Nx]./(rp.*(DXp+DXw))

# calculate the coefficients for the internal cells
AE = reshape(ue,Nx)
AW = reshape(-uw,Nx)
APx = reshape((ue.*DXe-uw.*DXw)./DXp,Nx)

# build the sparse matrix based on the numbering system
rowx_index = reshape(G[2:Nx+1],Nx) # main diagonal x
iix[1:3*Nx] = repmat(rowx_index,3)
jjx[1:3*Nx] = [reshape(G[1:Nx],Nx); reshape(G[2:Nx+1],Nx); reshape(G[3:Nx+2],Nx)]
sx[1:3*Nx] = [AW; APx; AE]

# build the sparse matrix
kx = 3*Nx
M = sparse(iix[1:kx], jjx[1:kx], sx[1:kx], Nx+2, Nx+2)

end

# =============== 1D Convection Term Cartesian Upwind =================
function convectionUpwindTerm1D(u::FaceValue)
# u is a face variable

# extract data from the mesh structure
Nx = u.domain.dims[1]
G = [1:Nx+2;]
DXp = u.domain.cellsize.x[2:end-1]

# define the vectors to store the sparse matrix data
iix = zeros(Int64, 3*(Nx+2))
jjx = zeros(Int64, 3*(Nx+2))
sx = zeros(Float64, 3*(Nx+2))

# find the velocity direction for the upwind scheme
ue_min = min(u.xvalue[2:Nx+1],0.0)
ue_max = max(u.xvalue[2:Nx+1],0.0)
uw_min = min(u.xvalue[1:Nx],0.0)
uw_max = max(u.xvalue[1:Nx],0.0)

# calculate the coefficients for the internal cells
AE = reshape(ue_min./DXp,Nx)
AW = reshape(-uw_max./DXp,Nx)
APx = reshape((ue_max-uw_min)./DXp,Nx)

# correct for the cells next to the boundary
# Left boundary:
APx[1] = APx[1]-uw_max[1]/(2.0*DXp[1])
AW[1] = AW[1]/2.0
# Right boundary:
AE[end] = AE[end]/2.0
APx[end] = APx[end] + ue_min[end]/(2.0*DXp[end])

# build the sparse matrix based on the numbering system
rowx_index = reshape(G[2:Nx+1],Nx) # main diagonal x
iix[1:3*Nx] = repmat(rowx_index,3)
jjx[1:3*Nx] = [reshape(G[1:Nx],Nx); reshape(G[2:Nx+1],Nx); reshape(G[3:Nx+2],Nx)]
sx[1:3*Nx] = [AW; APx; AE]

# build the sparse matrix
kx = 3*Nx
M = sparse(iix[1:kx], jjx[1:kx], sx[1:kx], Nx+2, Nx+2)

end

# =================== 1D Convection Terms Cylindrical Upwind =====================
function convectionUpwindTermCylindrical1D(u::FaceValue)
# u is a face variable

# extract data from the mesh structure
Nx = u.domain.dims[1]
G = [1:Nx+2;]
DXp = u.domain.cellsize.x[2:end-1]
rp = u.domain.cellcenters.x
rf = u.domain.facecenters.x

# define the vectors to store the sparse matrix data
iix = zeros(Int64, 3*(Nx+2))
jjx = zeros(Int64, 3*(Nx+2))
sx = zeros(Float64, 3*(Nx+2))

# reassign the east, west for code readability
re = rf[2:Nx+1]
rw = rf[1:Nx]

# find the velocity direction for the upwind scheme
ue_min = min(u.xvalue[2:Nx+1],0.0)
ue_max = max(u.xvalue[2:Nx+1],0.0)
uw_min = min(u.xvalue[1:Nx],0.0)
uw_max = max(u.xvalue[1:Nx],0.0)

# calculate the coefficients for the internal cells
AE = reshape(re.*ue_min./(DXp.*rp),Nx)
AW = reshape(-rw.*uw_max./(DXp.*rp),Nx)
APx = reshape((re.*ue_max-rw.*uw_min)./(DXp.*rp),Nx)

# correct for the cells next to the boundary
# Left boundary:
APx[1] = APx[1]-rw[1]*uw_max[1]/(2.0*DXp[1]*rp[1])
AW[1] = AW[1]/2.0
# Right boundary:
AE[end] = AE[end]/2.0
APx[end] = APx[end] + re[end]*ue_min[end]/(2.0*DXp[end]*rp[end])

# build the sparse matrix based on the numbering system
rowx_index = reshape(G[2:Nx+1],Nx) # main diagonal x
iix[1:3*Nx] = repmat(rowx_index,3)
jjx[1:3*Nx] = [reshape(G[1:Nx],Nx); reshape(G[2:Nx+1],Nx); reshape(G[3:Nx+2],Nx)]
sx[1:3*Nx] = [AW; APx; AE]

# build the sparse matrix
kx = 3*Nx
M = sparse(iix[1:kx], jjx[1:kx], sx[1:kx], Nx+2, Nx+2)

end

# ================= 1D Convection Terms Cylindrical TVD =====================
function convectionTvdTermCylindrical1D(u::FaceValue, phi::CellValue, FL::Function)
# u is a face variable
# phi is a cell variable

# a function to avoid division by zero
eps1 = 1.0e-20
fsign(phi_in) = (abs(phi_in).>=eps1).*phi_in+eps1*(phi_in.==0.0)+eps1*(abs(phi_in).<eps1).*sign(phi_in)

# extract data from the mesh structure
Nx = u.domain.dims[1]
G = [1:Nx+2;]
DXp = u.domain.cellsize.x[2:end-1]
dx = 0.5*(u.domain.cellsize.x[1:end-1]+u.domain.cellsize.x[2:end])
RHS = zeros(Float64, Nx+2)
psi_p = zeros(Float64, Nx+1)
psi_m = zeros(Float64, Nx+1)
r = u.domain.cellcenters.x
rf = u.domain.facecenters.x

# define the vectors to store the sparse matrix data
iix = zeros(Int64, 3*(Nx+2))
jjx = zeros(Int64, 3*(Nx+2))
sx = zeros(Float64, 3*(Nx+2))

# calculate the upstream to downstream gradient ratios for u>0 (+ ratio)
dphi_p = (phi.value[2:Nx+2]-phi.value[1:Nx+1])./dx
rp = dphi_p[1:end-1]./fsign(dphi_p[2:end])
psi_p[2:Nx+1] = 0.5*FL(rp).*(phi.value[3:Nx+2]-phi.value[2:Nx+1])
psi_p[1] = 0.0 # left boundary will be handled explicitly

# calculate the upstream to downstream gradient ratios for u<0 (- ratio)
rm = dphi_p[2:end]./fsign(dphi_p[1:end-1])
psi_m[1:Nx] = 0.5*FL(rm).*(phi.value[1:Nx]-phi.value[2:Nx+1])
psi_m[Nx+1] = 0.0 # right boundary will be handled explicitly

# reassign the east, west for code readability
re = rf[2:Nx+1]
rw = rf[1:Nx]

# find the velocity direction for the upwind scheme
ue_min = min(u.xvalue[2:Nx+1],0.0)
ue_max = max(u.xvalue[2:Nx+1],0.0)
uw_min = min(u.xvalue[1:Nx],0.0)
uw_max = max(u.xvalue[1:Nx],0.0)

# calculate the TVD correction term
RHS[2:Nx+1] = -(1.0./(DXp.*r)).*(re.*(ue_max.*psi_p[2:Nx+1]+ue_min.*psi_m[2:Nx+1])-
              rw.*(uw_max.*psi_p[1:Nx]+uw_min.*psi_m[1:Nx]))

# calculate the coefficients for the internal cells
AE = reshape(re.*ue_min./(DXp.*r),Nx)
AW = reshape(-rw.*uw_max./(DXp.*r),Nx)
APx = reshape((re.*ue_max-rw.*uw_min)./(DXp.*r),Nx)

# correct for the cells next to the boundary
# Left boundary:
APx[1] = APx[1]-rw[1]*uw_max[1]/(2.0*DXp[1]*r[1])
AW[1] = AW[1]/2.0
# Right boundary:
AE[end] = AE[end]/2.0
APx[end] = APx[end] + re[end]*ue_min[end]/(2.0*DXp[end]*r[end])

# build the sparse matrix based on the numbering system
rowx_index = reshape(G[2:Nx+1],Nx) # main diagonal x
iix[1:3*Nx] = repmat(rowx_index,3)
jjx[1:3*Nx] = [reshape(G[1:Nx],Nx); reshape(G[2:Nx+1],Nx); reshape(G[3:Nx+2],Nx)]
sx[1:3*Nx] = [AW; APx; AE]

# build the sparse matrix
kx = 3*Nx
M = sparse(iix[1:kx], jjx[1:kx], sx[1:kx], Nx+2, Nx+2)

(M, RHS)
end

# ================= 1D Convection Terms Cartesian TVD =====================
function convectionTvdTerm1D(u::FaceValue, phi::CellValue, FL::Function)
# u is a face variable
# phi is a cell variable

# a function to avoid division by zero
eps1 = 1.0e-20
fsign(phi_in) = (abs(phi_in).>=eps1).*phi_in+eps1*(phi_in.==0.0)+eps1*(abs(phi_in).<eps1).*sign(phi_in)

# extract data from the mesh structure
Nx = u.domain.dims[1]
G = [1:Nx+2;]
DXp = u.domain.cellsize.x[2:end-1]
dx = 0.5*(u.domain.cellsize.x[1:end-1]+u.domain.cellsize.x[2:end])
RHS = zeros(Float64, Nx+2)
psi_p = zeros(Float64, Nx+1)
psi_m = zeros(Float64, Nx+1)

# define the vectors to store the sparse matrix data
iix = zeros(Int64, 3*(Nx+2))
jjx = zeros(Int64, 3*(Nx+2))
sx = zeros(Float64, 3*(Nx+2))

# calculate the upstream to downstream gradient ratios for u>0 (+ ratio)
dphi_p = (phi.value[2:Nx+2]-phi.value[1:Nx+1])./dx
rp = dphi_p[1:end-1]./fsign(dphi_p[2:end])
psi_p[2:Nx+1] = 0.5*FL(rp).*(phi.value[3:Nx+2]-phi.value[2:Nx+1])
psi_p[1] = 0.0 # left boundary will be handled explicitly

# calculate the upstream to downstream gradient ratios for u<0 (- ratio)
rm = dphi_p[2:end]./fsign(dphi_p[1:end-1])
psi_m[1:Nx] = 0.5*FL(rm).*(phi.value[1:Nx]-phi.value[2:Nx+1])
psi_m[Nx+1] = 0.0 # right boundary will be handled explicitly

# find the velocity direction for the upwind scheme
ue_min = min(u.xvalue[2:Nx+1],0.0)
ue_max = max(u.xvalue[2:Nx+1],0.0)
uw_min = min(u.xvalue[1:Nx],0.0)
uw_max = max(u.xvalue[1:Nx],0.0)

# calculate the TVD correction term
RHS[2:Nx+1] = -(1.0./DXp).*((ue_max.*psi_p[2:Nx+1]+ue_min.*psi_m[2:Nx+1])-
              (uw_max.*psi_p[1:Nx]+uw_min.*psi_m[1:Nx]))

# calculate the coefficients for the internal cells
AE = reshape(ue_min./DXp,Nx)
AW = reshape(-uw_max./DXp,Nx)
APx = reshape((ue_max-uw_min)./DXp,Nx)

# correct for the cells next to the boundary
# Left boundary:
APx[1] = APx[1]-uw_max[1]/(2.0*DXp[1])
AW[1] = AW[1]/2.0
# Right boundary:
AE[end] = AE[end]/2.0
APx[end] = APx[end] + ue_min[end]/(2.0*DXp[end])

# build the sparse matrix based on the numbering system
rowx_index = reshape(G[2:Nx+1],Nx) # main diagonal x
iix[1:3*Nx] = repmat(rowx_index,3)
jjx[1:3*Nx] = [reshape(G[1:Nx],Nx); reshape(G[2:Nx+1],Nx); reshape(G[3:Nx+2],Nx)]
sx[1:3*Nx] = [AW; APx; AE]

# build the sparse matrix
kx = 3*Nx
M = sparse(iix[1:kx], jjx[1:kx], sx[1:kx], Nx+2, Nx+2)

(M, RHS)
end

# =================== 2D Convection Terms =======================
# ================ 2D Convection Term Central ===================
function convectionTerm2D(u::FaceValue)
# u is a face variable
# extract data from the mesh structure
Nx = u.domain.dims[1]
Ny = u.domain.dims[2]
G=reshape([1:(Nx+2)*(Ny+2);], Nx+2, Ny+2)
DXe = u.domain.cellsize.x[3:end]
DXw = u.domain.cellsize.x[1:end-2]
DXp = u.domain.cellsize.x[2:end-1]
DY = Array(Float64, 1, Ny+2)
DY[:] = u.domain.cellsize.y
DYn = DY[1,3:end]
DYs = DY[1,1:end-2]
DYp = DY[1,2:end-1]

# define the vectors to store the sparse matrix data
iix = zeros(Int64, 3*(Nx+2)*(Ny+2))
iiy = zeros(Int64, 3*(Nx+2)*(Ny+2))
jjx = zeros(Int64, 3*(Nx+2)*(Ny+2))
jjy = zeros(Int64, 3*(Nx+2)*(Ny+2))
sx = zeros(Float64, 3*(Nx+2)*(Ny+2))
sy = zeros(Float64, 3*(Nx+2)*(Ny+2))
mnx = Nx*Ny
mny = Nx*Ny

# reassign the east, west for code readability
ue = u.xvalue[2:Nx+1,:]./(DXp+DXe)
uw = u.xvalue[1:Nx,:]./(DXp+DXw)
vn = u.yvalue[:,2:Ny+1]./(DYp+DYn)
vs = u.yvalue[:,1:Ny]./(DYp+DYs)

# calculate the coefficients for the internal cells
AE = reshape(ue,mnx)
AW = reshape(-uw,mnx)
AN = reshape(vn,mny)
AS = reshape(-vs,mny)
APx = reshape((ue.*DXe-uw.*DXw)./DXp,mnx)
APy = reshape((vn.*DYn-vs.*DYs)./DYp,mny)

# build the sparse matrix based on the numbering system
rowx_index = reshape(G[2:Nx+1,2:Ny+1],mnx) # main diagonal x
iix[1:3*mnx] = repmat(rowx_index,3)
rowy_index = reshape(G[2:Nx+1,2:Ny+1],mny) # main diagonal y
iiy[1:3*mny] = repmat(rowy_index,3)
jjx[1:3*mnx] = [reshape(G[1:Nx,2:Ny+1],mnx); reshape(G[2:Nx+1,2:Ny+1],mnx); reshape(G[3:Nx+2,2:Ny+1],mnx)]
jjy[1:3*mny] = [reshape(G[2:Nx+1,1:Ny],mny); reshape(G[2:Nx+1,2:Ny+1],mny); reshape(G[2:Nx+1,3:Ny+2],mny)]
sx[1:3*mnx] = [AW; APx; AE]
sy[1:3*mny] = [AS; APy; AN]

# build the sparse matrix
kx = 3*mnx
ky = 3*mny
Mx = sparse(iix[1:kx], jjx[1:kx], sx[1:kx], (Nx+2)*(Ny+2), (Nx+2)*(Ny+2))
My = sparse(iiy[1:ky], jjy[1:ky], sy[1:ky], (Nx+2)*(Ny+2), (Nx+2)*(Ny+2))
M = Mx + My
(M, Mx, My)
end

# ================ 2D Convection Term Upwind ===================
function convectionUpwindTerm2D(u::FaceValue)
# u is a face variable
# extract data from the mesh structure
Nx = u.domain.dims[1]
Ny = u.domain.dims[2]
G=reshape([1:(Nx+2)*(Ny+2);], Nx+2, Ny+2)
DXp = u.domain.cellsize.x[2:end-1]
DYp = Array(Float64, 1, Ny)
DYp[:] = u.domain.cellsize.y[2:end-1]

# define the vectors to store the sparse matrix data
iix = zeros(Int64, 3*(Nx+2)*(Ny+2))
iiy = zeros(Int64, 3*(Nx+2)*(Ny+2))
jjx = zeros(Int64, 3*(Nx+2)*(Ny+2))
jjy = zeros(Int64, 3*(Nx+2)*(Ny+2))
sx = zeros(Float64, 3*(Nx+2)*(Ny+2))
sy = zeros(Float64, 3*(Nx+2)*(Ny+2))
mnx = Nx*Ny
mny = Nx*Ny

# find the velocity direction for the upwind scheme
ue_min = min(u.xvalue[2:Nx+1,:],0.0)
ue_max = max(u.xvalue[2:Nx+1,:],0.0)
uw_min = min(u.xvalue[1:Nx,:],0.0)
uw_max = max(u.xvalue[1:Nx,:],0.0)
vn_min = min(u.yvalue[:,2:Ny+1],0.0)
vn_max = max(u.yvalue[:,2:Ny+1],0.0)
vs_min = min(u.yvalue[:,1:Ny],0.0)
vs_max = max(u.yvalue[:,1:Ny],0.0)

# calculate the coefficients for the internal cells, not reshape
AE = ue_min./DXp
AW = -uw_max./DXp
AN = vn_min./DYp
AS = -vs_max./DYp
APx = (ue_max-uw_min)./DXp
APy = (vn_max-vs_min)./DYp

# Also correct for the boundary cells (not the ghost cells)
# Left boundary:
APx[1,:] = APx[1,:]-uw_max[1,:]/(2.0*DXp[1])
AW[1,:] = AW[1,:]/2.0
# Right boundary:
AE[end,:] = AE[end,:]/2.0
APx[end,:] = APx[end,:]+ue_min[end,:]/(2.0*DXp[end])
# Bottom boundary:
APy[:,1] = APy[:,1]-vs_max[:,1]/(2.0*DYp[1])
AS[:,1] = AS[:,1]/2.0
# Top boundary:
AN[:,end] = AN[:,end]/2.0
APy[:,end] = APy[:,end]+vn_min[:,end]/(2.0*DYp[end])

# now reshape
AE = reshape(AE,mnx)
AW = reshape(AW,mnx)
AN = reshape(AN,mny)
AS = reshape(AS,mny)
APx = reshape(APx,mnx)
APy = reshape(APy,mny)

# build the sparse matrix based on the numbering system
rowx_index = reshape(G[2:Nx+1,2:Ny+1],mnx) # main diagonal x
iix[1:3*mnx] = repmat(rowx_index,3)
rowy_index = reshape(G[2:Nx+1,2:Ny+1],mny) # main diagonal y
iiy[1:3*mny] = repmat(rowy_index,3)
jjx[1:3*mnx] = [reshape(G[1:Nx,2:Ny+1],mnx); reshape(G[2:Nx+1,2:Ny+1],mnx); reshape(G[3:Nx+2,2:Ny+1],mnx)]
jjy[1:3*mny] = [reshape(G[2:Nx+1,1:Ny],mny); reshape(G[2:Nx+1,2:Ny+1],mny); reshape(G[2:Nx+1,3:Ny+2],mny)]
sx[1:3*mnx] = [AW; APx; AE]
sy[1:3*mny] = [AS; APy; AN]

# build the sparse matrix
kx = 3*mnx
ky = 3*mny
Mx = sparse(iix[1:kx], jjx[1:kx], sx[1:kx], (Nx+2)*(Ny+2), (Nx+2)*(Ny+2))
My = sparse(iiy[1:ky], jjy[1:ky], sy[1:ky], (Nx+2)*(Ny+2), (Nx+2)*(Ny+2))
M = Mx + My
(M, Mx, My)
end

# ================ 2D Convection Term TVD =======================
function convectionTvdTerm2D(u::FaceValue, phi::CellValue, FL::Function)
# u is a face variable
# phi is a cell variable

# a function to avoid division by zero
eps1 = 1.0e-20
fsign(phi_in) = (abs(phi_in).>=eps1).*phi_in+eps1*(phi_in.==0.0)+eps1*(abs(phi_in).<eps1).*sign(phi_in)

# extract data from the mesh structure
Nx = u.domain.dims[1]
Ny = u.domain.dims[2]
G=reshape([1:(Nx+2)*(Ny+2);], Nx+2, Ny+2)
DXp = u.domain.cellsize.x[2:end-1]
DYp = Array(Float64, 1, Ny)
DYp[:] = u.domain.cellsize.y[2:end-1]
dx=0.5*(u.domain.cellsize.x[1:end-1]+u.domain.cellsize.x[2:end])
dy=Array(Float64, 1, Ny+1)
dy[:]=0.5*(u.domain.cellsize.y[1:end-1]+u.domain.cellsize.y[2:end])

psiX_p = zeros(Nx+1,Ny)
psiX_m = zeros(Nx+1,Ny)
psiY_p = zeros(Nx,Ny+1)
psiY_m = zeros(Nx,Ny+1)

# define the vectors to store the sparse matrix data
iix = zeros(Int64, 3*(Nx+2)*(Ny+2))
iiy = zeros(Int64, 3*(Nx+2)*(Ny+2))
jjx = zeros(Int64, 3*(Nx+2)*(Ny+2))
jjy = zeros(Int64, 3*(Nx+2)*(Ny+2))
sx = zeros(Float64, 3*(Nx+2)*(Ny+2))
sy = zeros(Float64, 3*(Nx+2)*(Ny+2))
mnx = Nx*Ny
mny = Nx*Ny

# calculate the upstream to downstream gradient ratios for u>0 (+ ratio)
# x direction
dphiX_p = (phi.value[2:Nx+2, 2:Ny+1]-phi.value[1:Nx+1, 2:Ny+1])./dx
rX_p = dphiX_p[1:end-1,:]./fsign(dphiX_p[2:end,:])
psiX_p[2:Nx+1,:] = 0.5*FL(rX_p).*(phi.value[3:Nx+2,2:Ny+1]-
		    phi.value[2:Nx+1, 2:Ny+1])
psiX_p[1, :] = 0.0 # left boundary will be handled in the main matrix
# y direction
dphiY_p = (phi.value[2:Nx+1, 2:Ny+2]-phi.value[2:Nx+1, 1:Ny+1])./dy
rY_p = dphiY_p[:,1:end-1]./fsign(dphiY_p[:,2:end])
psiY_p[:,2:Ny+1] = 0.5*FL(rY_p).*(phi.value[2:Nx+1,3:Ny+2]-
		  phi.value[2:Nx+1, 2:Ny+1])
psiY_p[:,1] = 0.0 # Bottom boundary will be handled in the main matrix

# calculate the upstream to downstream gradient ratios for u<0 (- ratio)
# x direction
rX_m = dphiX_p[2:end,:]./fsign(dphiX_p[1:end-1,:])
psiX_m[1:Nx,:] = 0.5*FL(rX_m).*(phi.value[1:Nx, 2:Ny+1]-
		phi.value[2:Nx+1, 2:Ny+1])
psiX_m[Nx+1,:] = 0.0 # right boundary
# y direction
rY_m = dphiY_p[:,2:end]./fsign(dphiY_p[:,1:end-1])
psiY_m[:,1:Ny] = 0.5*FL(rY_m).*(phi.value[2:Nx+1, 1:Ny]-
	      phi.value[2:Nx+1, 2:Ny+1])
psiY_m[:, Ny+1] = 0.0 # top boundary will be handled in the main matrix

# find the velocity direction for the upwind scheme
ue_min = min(u.xvalue[2:Nx+1,:],0.0)
ue_max = max(u.xvalue[2:Nx+1,:],0.0)
uw_min = min(u.xvalue[1:Nx,:],0.0)
uw_max = max(u.xvalue[1:Nx,:],0.0)
vn_min = min(u.yvalue[:,2:Ny+1],0.0)
vn_max = max(u.yvalue[:,2:Ny+1],0.0)
vs_min = min(u.yvalue[:,1:Ny],0.0)
vs_max = max(u.yvalue[:,1:Ny],0.0)

# calculate the coefficients for the internal cells, not reshape
AE = ue_min./DXp
AW = -uw_max./DXp
AN = vn_min./DYp
AS = -vs_max./DYp
APx = (ue_max-uw_min)./DXp
APy = (vn_max-vs_min)./DYp

# Also correct for the boundary cells (not the ghost cells)
# Left boundary:
APx[1,:] = APx[1,:]-uw_max[1,:]/(2.0*DXp[1])
AW[1,:] = AW[1,:]/2.0
# Right boundary:
AE[end,:] = AE[end,:]/2.0
APx[end,:] = APx[end,:]+ue_min[end,:]/(2.0*DXp[end])
# Bottom boundary:
APy[:,1] = APy[:,1]-vs_max[:,1]/(2.0*DYp[1])
AS[:,1] = AS[:,1]/2.0
# Top boundary:
AN[:,end] = AN[:,end]/2.0
APy[:,end] = APy[:,end]+vn_min[:,end]/(2.0*DYp[end])

# now reshape
AE = reshape(AE,mnx)
AW = reshape(AW,mnx)
AN = reshape(AN,mny)
AS = reshape(AS,mny)
APx = reshape(APx,mnx)
APy = reshape(APy,mny)

# build the sparse matrix based on the numbering system
rowx_index = reshape(G[2:Nx+1,2:Ny+1],mnx) # main diagonal x
iix[1:3*mnx] = repmat(rowx_index,3)
rowy_index = reshape(G[2:Nx+1,2:Ny+1],mny) # main diagonal y
iiy[1:3*mny] = repmat(rowy_index,3)
jjx[1:3*mnx] = [reshape(G[1:Nx,2:Ny+1],mnx); reshape(G[2:Nx+1,2:Ny+1],mnx); reshape(G[3:Nx+2,2:Ny+1],mnx)]
jjy[1:3*mny] = [reshape(G[2:Nx+1,1:Ny],mny); reshape(G[2:Nx+1,2:Ny+1],mny); reshape(G[2:Nx+1,3:Ny+2],mny)]
sx[1:3*mnx] = [AW; APx; AE]
sy[1:3*mny] = [AS; APy; AN]

# calculate the TVD correction term
div_x = -(1./DXp).*((ue_max.*psiX_p[2:Nx+1,:]+ue_min.*psiX_m[2:Nx+1,:])-
              (uw_max.*psiX_p[1:Nx,:]+uw_min.*psiX_m[1:Nx,:]))
div_y = -(1./DYp).*((vn_max.*psiY_p[:,2:Ny+1]+vn_min.*psiY_m[:,2:Ny+1])-
              (vs_max.*psiY_p[:,1:Ny]+vs_min.*psiY_m[:,1:Ny]))

# define the RHS Vector
RHS = zeros((Nx+2)*(Ny+2))
RHSx = zeros((Nx+2)*(Ny+2))
RHSy = zeros((Nx+2)*(Ny+2))

# assign the values of the RHS vector
RHS[rowx_index] = reshape(div_x+div_y,Nx*Ny)
RHSx[rowx_index] = reshape(div_x,Nx*Ny)
RHSy[rowy_index] = reshape(div_y,Nx*Ny)

# build the sparse matrix
kx = 3*mnx
ky = 3*mny
Mx = sparse(iix[1:kx], jjx[1:kx], sx[1:kx], (Nx+2)*(Ny+2), (Nx+2)*(Ny+2))
My = sparse(iiy[1:ky], jjy[1:ky], sy[1:ky], (Nx+2)*(Ny+2), (Nx+2)*(Ny+2))
M = Mx + My
(M, RHS, Mx, My, RHSx, RHSy)
end

# ================ 2D Convection Term Cylindrical Central ===================
function convectionTermCylindrical2D(u::FaceValue)
# u is a face variable
# extract data from the mesh structure
Nr = u.domain.dims[1]
Nz = u.domain.dims[2]
G=reshape([1:(Nr+2)*(Nz+2);], Nr+2, Nz+2)
DXe = u.domain.cellsize.x[3:end]
DXw = u.domain.cellsize.x[1:end-2]
DXp = u.domain.cellsize.x[2:end-1]
DZ = Array(Float64, 1, Nz+2)
DZ[:] = u.domain.cellsize.y
DYn = DZ[1,3:end]
DYs = DZ[1,1:end-2]
DYp = DZ[1,2:end-1]
rp = repmat(u.domain.cellcenters.x, 1, Nz)
rf = repmat(u.domain.facecenters.x, 1, Nz)

# define the vectors to store the sparse matrix data
iix = zeros(Int64, 3*(Nr+2)*(Nz+2))
iiy = zeros(Int64, 3*(Nr+2)*(Nz+2))
jjx = zeros(Int64, 3*(Nr+2)*(Nz+2))
jjy = zeros(Int64, 3*(Nr+2)*(Nz+2))
sx = zeros(Float64, 3*(Nr+2)*(Nz+2))
sy = zeros(Float64, 3*(Nr+2)*(Nz+2))
mnx = Nr*Nz
mny = Nr*Nz

# reassign the east, west for code readability
ue = rf[2:Nr+1,:].*u.xvalue[2:Nr+1,:]./(rp.*(DXp+DXe))
uw = rf[1:Nr,:].*u.xvalue[1:Nr,:]./(rp.*(DXp+DXw))
vn = u.yvalue[:,2:Nz+1]./(DYp+DYn)
vs = u.yvalue[:,1:Nz]./(DYp+DYs)
re = rf[2:Nr+1,:]
rw = rf[1:Nr,:]

# calculate the coefficients for the internal cells
AE = reshape(ue,mnx)
AW = reshape(-uw,mnx)
AN = reshape(vn,mny)
AS = reshape(-vs,mny)
APx = reshape((ue.*DXe-uw.*DXw)./DXp,mnx)
APy = reshape((vn.*DYn-vs.*DYs)./DYp,mny)

# build the sparse matrix based on the numbering system
rowx_index = reshape(G[2:Nr+1,2:Nz+1],mnx) # main diagonal x
iix[1:3*mnx] = repmat(rowx_index,3)
rowy_index = reshape(G[2:Nr+1,2:Nz+1],mny) # main diagonal y
iiy[1:3*mny] = repmat(rowy_index,3)
jjx[1:3*mnx] = [reshape(G[1:Nr,2:Nz+1],mnx); reshape(G[2:Nr+1,2:Nz+1],mnx); reshape(G[3:Nr+2,2:Nz+1],mnx)]
jjy[1:3*mny] = [reshape(G[2:Nr+1,1:Nz],mny); reshape(G[2:Nr+1,2:Nz+1],mny); reshape(G[2:Nr+1,3:Nz+2],mny)]
sx[1:3*mnx] = [AW; APx; AE]
sy[1:3*mny] = [AS; APy; AN]

# build the sparse matrix
kx = 3*mnx
ky = 3*mny
Mx = sparse(iix[1:kx], jjx[1:kx], sx[1:kx], (Nr+2)*(Nz+2), (Nr+2)*(Nz+2))
My = sparse(iiy[1:ky], jjy[1:ky], sy[1:ky], (Nr+2)*(Nz+2), (Nr+2)*(Nz+2))
M = Mx + My
(M, Mx, My)
end

# ================ 2D Convection Term Cylindrical Upwind ==============================
function convectionUpwindTermCylindrical2D(u::FaceValue)
# u is a face variable
# extract data from the mesh structure
Nr = u.domain.dims[1]
Nz = u.domain.dims[2]
G=reshape([1:(Nr+2)*(Nz+2);], Nr+2, Nz+2)
DRp = u.domain.cellsize.x[2:end-1]
DZp = Array(Float64, 1, Nz)
DZp[:] = u.domain.cellsize.y[2:end-1]
rp = repmat(u.domain.cellcenters.x, 1, Nz)
rf = repmat(u.domain.facecenters.x, 1, Nz)

# define the vectors to store the sparse matrix data
iix = zeros(Int64, 3*(Nr+2)*(Nz+2))
iiy = zeros(Int64, 3*(Nr+2)*(Nz+2))
jjx = zeros(Int64, 3*(Nr+2)*(Nz+2))
jjy = zeros(Int64, 3*(Nr+2)*(Nz+2))
sx = zeros(Float64, 3*(Nr+2)*(Nz+2))
sy = zeros(Float64, 3*(Nr+2)*(Nz+2))
mnx = Nr*Nz
mny = Nr*Nz

re = rf[2:Nr+1,:]
rw = rf[1:Nr,:]

# find the velocity direction for the upwind scheme
ue_min = min(u.xvalue[2:Nr+1,:],0.0)
ue_max = max(u.xvalue[2:Nr+1,:],0.0)
uw_min = min(u.xvalue[1:Nr,:],0.0)
uw_max = max(u.xvalue[1:Nr,:],0.0)
vn_min = min(u.yvalue[:,2:Nz+1],0.0)
vn_max = max(u.yvalue[:,2:Nz+1],0.0)
vs_min = min(u.yvalue[:,1:Nz],0.0)
vs_max = max(u.yvalue[:,1:Nz],0.0)

# calculate the coefficients for the internal cells, do not reshape yet
AE = re.*ue_min./(DRp.*rp)
AW = -rw.*uw_max./(DRp.*rp)
AN = vn_min./DZp
AS = -vs_max./DZp
APx = (re.*ue_max-rw.*uw_min)./(DRp.*rp)
APy = (vn_max-vs_min)./DZp

# Also correct for the boundary cells (not the ghost cells)
# Left boundary:
APx[1,:] = APx[1,:]-rw[1,:].*uw_max[1,:]./(2.0*DRp[1]*rp[1,:])
AW[1,:] = AW[1,:]/2.0
# Right boundary:
AE[end,:] = AE[end,:]/2.0
APx[end,:] = APx[end,:]+re[end,:].*ue_min[end,:]./(2.0*DRp[end]*rp[end,:])
# Bottom boundary:
APy[:,1] = APy[:,1]-vs_max[:,1]/(2.0*DZp[1])
AS[:,1] = AS[:,1]/2.0
# Top boundary:
AN[:,end] = AN[:,end]/2.0
APy[:,end] = APy[:,end]+vn_min[:,end]/(2.0*DZp[end])

# now reshape
AE = reshape(AE,mnx)
AW = reshape(AW,mnx)
AN = reshape(AN,mny)
AS = reshape(AS,mny)
APx = reshape(APx,mnx)
APy = reshape(APy,mny)

# build the sparse matrix based on the numbering system
rowx_index = reshape(G[2:Nr+1,2:Nz+1],mnx) # main diagonal x
iix[1:3*mnx] = repmat(rowx_index,3)
rowy_index = reshape(G[2:Nr+1,2:Nz+1],mny) # main diagonal y
iiy[1:3*mny] = repmat(rowy_index,3)
jjx[1:3*mnx] = [reshape(G[1:Nr,2:Nz+1],mnx); reshape(G[2:Nr+1,2:Nz+1],mnx); reshape(G[3:Nr+2,2:Nz+1],mnx)]
jjy[1:3*mny] = [reshape(G[2:Nr+1,1:Nz],mny); reshape(G[2:Nr+1,2:Nz+1],mny); reshape(G[2:Nr+1,3:Nz+2],mny)]
sx[1:3*mnx] = [AW; APx; AE]
sy[1:3*mny] = [AS; APy; AN]

# build the sparse matrix
kx = 3*mnx
ky = 3*mny
Mx = sparse(iix[1:kx], jjx[1:kx], sx[1:kx], (Nr+2)*(Nz+2), (Nr+2)*(Nz+2))
My = sparse(iiy[1:ky], jjy[1:ky], sy[1:ky], (Nr+2)*(Nz+2), (Nr+2)*(Nz+2))
M = Mx + My
(M, Mx, My)
end

# ==================== 2D Convection Term Cylindrical TVD ===========================
function convectionTvdTermCylindrical2D(u::FaceValue, phi::CellValue, FL::Function)
# u is a face variable
# phi is a cell variable

# a function to avoid division by zero
eps1 = 1.0e-20
fsign(phi_in) = (abs(phi_in).>=eps1).*phi_in+eps1*(phi_in.==0.0)+eps1*(abs(phi_in).<eps1).*sign(phi_in)

# extract data from the mesh structure
Nr = u.domain.dims[1]
Nz = u.domain.dims[2]
G=reshape([1:(Nr+2)*(Nz+2);], Nr+2, Nz+2)
DRp = u.domain.cellsize.x[2:end-1]
DZp = Array(Float64, 1, Nz)
DZp[:] = u.domain.cellsize.y[2:end-1]
dr=0.5*(u.domain.cellsize.x[1:end-1]+u.domain.cellsize.x[2:end])
dz=Array(Float64, 1, Nz+1)
dz[:]=0.5*(u.domain.cellsize.y[1:end-1]+u.domain.cellsize.y[2:end])
rp = repmat(u.domain.cellcenters.x, 1, Nz)
rf = repmat(u.domain.facecenters.x, 1, Nz)
psiX_p = zeros(Nr+1,Nz)
psiX_m = zeros(Nr+1,Nz)
psiY_p = zeros(Nr,Nz+1)
psiY_m = zeros(Nr,Nz+1)

# define the vectors to store the sparse matrix data
iix = zeros(Int64, 3*(Nr+2)*(Nz+2))
iiy = zeros(Int64, 3*(Nr+2)*(Nz+2))
jjx = zeros(Int64, 3*(Nr+2)*(Nz+2))
jjy = zeros(Int64, 3*(Nr+2)*(Nz+2))
sx = zeros(Float64, 3*(Nr+2)*(Nz+2))
sy = zeros(Float64, 3*(Nr+2)*(Nz+2))
mnx = Nr*Nz
mny = Nr*Nz

# calculate the upstream to downstream gradient ratios for u>0 (+ ratio)
# x direction
dphiX_p = (phi.value[2:Nr+2, 2:Nz+1]-phi.value[1:Nr+1, 2:Nz+1])./dr
rX_p = dphiX_p[1:end-1,:]./fsign(dphiX_p[2:end,:])
psiX_p[2:Nr+1,:] = 0.5*FL(rX_p).*(phi.value[3:Nr+2,2:Nz+1]-
		    phi.value[2:Nr+1, 2:Nz+1])
psiX_p[1, :] = 0.0 # left boundary will be handled in the main matrix
# y direction
dphiY_p = (phi.value[2:Nr+1, 2:Nz+2]-phi.value[2:Nr+1, 1:Nz+1])./dz
rY_p = dphiY_p[:,1:end-1]./fsign(dphiY_p[:,2:end])
psiY_p[:,2:Nz+1] = 0.5*FL(rY_p).*(phi.value[2:Nr+1,3:Nz+2]-
		  phi.value[2:Nr+1, 2:Nz+1])
psiY_p[:,1] = 0.0 # Bottom boundary will be handled in the main matrix

# calculate the upstream to downstream gradient ratios for u<0 (- ratio)
# x direction
rX_m = dphiX_p[2:end,:]./fsign(dphiX_p[1:end-1,:])
psiX_m[1:Nr,:] = 0.5*FL(rX_m).*(phi.value[1:Nr, 2:Nz+1]-
		phi.value[2:Nr+1, 2:Nz+1])
psiX_m[Nr+1,:] = 0.0 # right boundary
# y direction
rY_m = dphiY_p[:,2:end]./fsign(dphiY_p[:,1:end-1])
psiY_m[:,1:Nz] = 0.5*FL(rY_m).*(phi.value[2:Nr+1, 1:Nz]-
	      phi.value[2:Nr+1, 2:Nz+1])
psiY_m[:, Nz+1] = 0.0 # top boundary will be handled in the main matrix

re = rf[2:Nr+1,:]
rw = rf[1:Nr,:]

# find the velocity direction for the upwind scheme
ue_min = min(u.xvalue[2:Nr+1,:],0.0)
ue_max = max(u.xvalue[2:Nr+1,:],0.0)
uw_min = min(u.xvalue[1:Nr,:],0.0)
uw_max = max(u.xvalue[1:Nr,:],0.0)
vn_min = min(u.yvalue[:,2:Nz+1],0.0)
vn_max = max(u.yvalue[:,2:Nz+1],0.0)
vs_min = min(u.yvalue[:,1:Nz],0.0)
vs_max = max(u.yvalue[:,1:Nz],0.0)

# calculate the coefficients for the internal cells, do not reshape yet
AE = re.*ue_min./(DRp.*rp)
AW = -rw.*uw_max./(DRp.*rp)
AN = vn_min./DZp
AS = -vs_max./DZp
APx = (re.*ue_max-rw.*uw_min)./(DRp.*rp)
APy = (vn_max-vs_min)./DZp

# Also correct for the boundary cells (not the ghost cells)
# Left boundary:
APx[1,:] = APx[1,:]-rw[1,:].*uw_max[1,:]./(2.0*DRp[1]*rp[1,:])
AW[1,:] = AW[1,:]/2.0
# Right boundary:
AE[end,:] = AE[end,:]/2.0
APx[end,:] = APx[end,:]+re[end,:].*ue_min[end,:]./(2.0*DRp[end]*rp[end,:])
# Bottom boundary:
APy[:,1] = APy[:,1]-vs_max[:,1]/(2.0*DZp[1])
AS[:,1] = AS[:,1]/2.0
# Top boundary:
AN[:,end] = AN[:,end]/2.0
APy[:,end] = APy[:,end]+vn_min[:,end]/(2.0*DZp[end])

# now reshape
AE = reshape(AE,mnx)
AW = reshape(AW,mnx)
AN = reshape(AN,mny)
AS = reshape(AS,mny)
APx = reshape(APx,mnx)
APy = reshape(APy,mny)

# build the sparse matrix based on the numbering system
rowx_index = reshape(G[2:Nr+1,2:Nz+1],mnx) # main diagonal x
iix[1:3*mnx] = repmat(rowx_index,3)
rowy_index = reshape(G[2:Nr+1,2:Nz+1],mny) # main diagonal y
iiy[1:3*mny] = repmat(rowy_index,3)
jjx[1:3*mnx] = [reshape(G[1:Nr,2:Nz+1],mnx); reshape(G[2:Nr+1,2:Nz+1],mnx); reshape(G[3:Nr+2,2:Nz+1],mnx)]
jjy[1:3*mny] = [reshape(G[2:Nr+1,1:Nz],mny); reshape(G[2:Nr+1,2:Nz+1],mny); reshape(G[2:Nr+1,3:Nz+2],mny)]
sx[1:3*mnx] = [AW; APx; AE]
sy[1:3*mny] = [AS; APy; AN]

# calculate the TVD correction term
div_x = -(1./(DRp.*rp)).*(re.*(ue_max.*psiX_p[2:Nr+1,:]+ue_min.*psiX_m[2:Nr+1,:])-
              rw.*(uw_max.*psiX_p[1:Nr,:]+uw_min.*psiX_m[1:Nr,:]))
div_y = -(1./DZp).*((vn_max.*psiY_p[:,2:Nz+1]+vn_min.*psiY_m[:,2:Nz+1])-
              (vs_max.*psiY_p[:,1:Nz]+vs_min.*psiY_m[:,1:Nz]))

# define the RHS Vector
RHS = zeros((Nr+2)*(Nz+2))
RHSx = zeros((Nr+2)*(Nz+2))
RHSy = zeros((Nr+2)*(Nz+2))

# assign the values of the RHS vector
RHS[rowx_index] = reshape(div_x+div_y,Nr*Nz)
RHSx[rowx_index] = reshape(div_x,Nr*Nz)
RHSy[rowy_index] = reshape(div_y,Nr*Nz)

# build the sparse matrix
kx = 3*mnx
ky = 3*mny
Mx = sparse(iix[1:kx], jjx[1:kx], sx[1:kx], (Nr+2)*(Nz+2), (Nr+2)*(Nz+2))
My = sparse(iiy[1:ky], jjy[1:ky], sy[1:ky], (Nr+2)*(Nz+2), (Nr+2)*(Nz+2))
M = Mx + My
(M, RHS, Mx, My, RHSx, RHSy)
end

# ================ 2D Convection Term Radial Central ===================
function convectionTermRadial2D(u::FaceValue)
# u is a face variable
# extract data from the mesh structure
Nr = u.domain.dims[1]
Ntheta = u.domain.dims[2]
G=reshape([1:(Nr+2)*(Ntheta+2);], Nr+2, Ntheta+2)
DRe = u.domain.cellsize.x[3:end]
DRw = u.domain.cellsize.x[1:end-2]
DRp = u.domain.cellsize.x[2:end-1]
DTHETA = Array(Float64, 1, Ntheta+2)
DTHETA[:] = u.domain.cellsize.y
DTHETAn = DTHETA[1,3:end]
DTHETAs = DTHETA[1,1:end-2]
DTHETAp = DTHETA[1,2:end-1]
rp = u.domain.cellcenters.x
rf = u.domain.facecenters.x

# define the vectors to store the sparse matrix data
iix = zeros(Int64, 3*(Nr+2)*(Ntheta+2))
iiy = zeros(Int64, 3*(Nr+2)*(Ntheta+2))
jjx = zeros(Int64, 3*(Nr+2)*(Ntheta+2))
jjy = zeros(Int64, 3*(Nr+2)*(Ntheta+2))
sx = zeros(Float64, 3*(Nr+2)*(Ntheta+2))
sy = zeros(Float64, 3*(Nr+2)*(Ntheta+2))
mnx = Nr*Ntheta
mny = Nr*Ntheta

re = rf[2:Nr+1,:]
rw = rf[1:Nr,:]

# reassign the east, west for code readability
ue = re.*u.xvalue[2:Nr+1,:]./(DRp+DRe)
uw = rw.*u.xvalue[1:Nr,:]./(DRp+DRw)
vn = u.yvalue[:,2:Ntheta+1]./(rp.*(DTHETAp+DTHETAn))
vs = u.yvalue[:,1:Ntheta]./(rp.*(DTHETAp+DTHETAs))

# calculate the coefficients for the internal cells
AE = reshape(ue,mnx)
AW = reshape(-uw,mnx)
AN = reshape(vn,mny)
AS = reshape(-vs,mny)
APx = reshape((ue.*DRe-uw.*DRw)./DRp,mnx)
APy = reshape((vn.*DTHETAn-vs.*DTHETAs)./DTHETAp,mny)

# build the sparse matrix based on the numbering system
rowx_index = reshape(G[2:Nr+1,2:Ntheta+1],mnx) # main diagonal x
iix[1:3*mnx] = repmat(rowx_index,3)
rowy_index = reshape(G[2:Nr+1,2:Ntheta+1],mny) # main diagonal y
iiy[1:3*mny] = repmat(rowy_index,3)
jjx[1:3*mnx] = [reshape(G[1:Nr,2:Ntheta+1],mnx); reshape(G[2:Nr+1,2:Ntheta+1],mnx); reshape(G[3:Nr+2,2:Ntheta+1],mnx)]
jjy[1:3*mny] = [reshape(G[2:Nr+1,1:Ntheta],mny); reshape(G[2:Nr+1,2:Ntheta+1],mny); reshape(G[2:Nr+1,3:Ntheta+2],mny)]
sx[1:3*mnx] = [AW; APx; AE]
sy[1:3*mny] = [AS; APy; AN]

# build the sparse matrix
kx = 3*mnx
ky = 3*mny
Mx = sparse(iix[1:kx], jjx[1:kx], sx[1:kx], (Nr+2)*(Ntheta+2), (Nr+2)*(Ntheta+2))
My = sparse(iiy[1:ky], jjy[1:ky], sy[1:ky], (Nr+2)*(Ntheta+2), (Nr+2)*(Ntheta+2))
M = Mx + My
(M, Mx, My)
end

# ================ 2D Convection Term Radial Upwind ==============================
function convectionUpwindTermRadial2D(u::FaceValue)
# u is a face variable
# extract data from the mesh structure
Nr = u.domain.dims[1]
Ntheta = u.domain.dims[2]
G=reshape([1:(Nr+2)*(Ntheta+2);], Nr+2, Ntheta+2)
DRp = u.domain.cellsize.x[2:end-1]
DTHETAp = Array(Float64, 1, Ntheta)
DTHETAp[:] = u.domain.cellsize.y[2:end-1]
rp = u.domain.cellcenters.x
rf = u.domain.facecenters.x

# define the vectors to store the sparse matrix data
iix = zeros(Int64, 3*(Nr+2)*(Ntheta+2))
iiy = zeros(Int64, 3*(Nr+2)*(Ntheta+2))
jjx = zeros(Int64, 3*(Nr+2)*(Ntheta+2))
jjy = zeros(Int64, 3*(Nr+2)*(Ntheta+2))
sx = zeros(Float64, 3*(Nr+2)*(Ntheta+2))
sy = zeros(Float64, 3*(Nr+2)*(Ntheta+2))
mnx = Nr*Ntheta
mny = Nr*Ntheta

re = rf[2:Nr+1,:]
rw = rf[1:Nr,:]

# find the velocity direction for the upwind scheme
ue_min = min(u.xvalue[2:Nr+1,:],0.0)
ue_max = max(u.xvalue[2:Nr+1,:],0.0)
uw_min = min(u.xvalue[1:Nr,:],0.0)
uw_max = max(u.xvalue[1:Nr,:],0.0)
vn_min = min(u.yvalue[:,2:Ntheta+1],0.0)
vn_max = max(u.yvalue[:,2:Ntheta+1],0.0)
vs_min = min(u.yvalue[:,1:Ntheta],0.0)
vs_max = max(u.yvalue[:,1:Ntheta],0.0)

# calculate the coefficients for the internal cells, do not reshape yet
AE = re.*ue_min./(DRp.*rp)
AW = -rw.*uw_max./(DRp.*rp)
AN = vn_min./(DTHETAp.*rp)
AS = -vs_max./(DTHETAp.*rp)
APx = (re.*ue_max-rw.*uw_min)./(DRp.*rp)
APy = (vn_max-vs_min)./(DTHETAp.*rp)

# Also correct for the boundary cells (not the ghost cells)
# Left boundary:
APx[1,:] = APx[1,:]-rw[1,:].*uw_max[1,:]./(2.0*DRp[1]*rp[1,:])
AW[1,:] = AW[1,:]/2.0
# Right boundary:
AE[end,:] = AE[end,:]/2.0
APx[end,:] = APx[end,:]+re[end,:].*ue_min[end,:]./(2.0*DRp[end]*rp[end,:])
# Bottom boundary:
APy[:,1] = APy[:,1]-vs_max[:,1]./(2.0*DTHETAp[1]*rp[:,1])
AS[:,1] = AS[:,1]/2.0
# Top boundary:
AN[:,end] = AN[:,end]/2.0
APy[:,end] = APy[:,end]+vn_min[:,end]./(2.0*DTHETAp[end]*rp[:,end])

# now reshape
AE = reshape(AE,mnx)
AW = reshape(AW,mnx)
AN = reshape(AN,mny)
AS = reshape(AS,mny)
APx = reshape(APx,mnx)
APy = reshape(APy,mny)

# build the sparse matrix based on the numbering system
rowx_index = reshape(G[2:Nr+1,2:Ntheta+1],mnx) # main diagonal x
iix[1:3*mnx] = repmat(rowx_index,3)
rowy_index = reshape(G[2:Nr+1,2:Ntheta+1],mny) # main diagonal y
iiy[1:3*mny] = repmat(rowy_index,3)
jjx[1:3*mnx] = [reshape(G[1:Nr,2:Ntheta+1],mnx); reshape(G[2:Nr+1,2:Ntheta+1],mnx); reshape(G[3:Nr+2,2:Ntheta+1],mnx)]
jjy[1:3*mny] = [reshape(G[2:Nr+1,1:Ntheta],mny); reshape(G[2:Nr+1,2:Ntheta+1],mny); reshape(G[2:Nr+1,3:Ntheta+2],mny)]
sx[1:3*mnx] = [AW; APx; AE]
sy[1:3*mny] = [AS; APy; AN]

# build the sparse matrix
kx = 3*mnx
ky = 3*mny
Mx = sparse(iix[1:kx], jjx[1:kx], sx[1:kx], (Nr+2)*(Ntheta+2), (Nr+2)*(Ntheta+2))
My = sparse(iiy[1:ky], jjy[1:ky], sy[1:ky], (Nr+2)*(Ntheta+2), (Nr+2)*(Ntheta+2))
M = Mx + My
(M, Mx, My)
end

# ================ 2D Convection Term Radial TVD =======================
function convectionTvdTermRadial2D(u::FaceValue, phi::CellValue, FL::Function)
# u is a face variable
# phi is a cell variable

# a function to avoid division by zero
eps1 = 1.0e-20
fsign(phi_in) = (abs(phi_in).>=eps1).*phi_in+eps1*(phi_in.==0.0)+eps1*(abs(phi_in).<eps1).*sign(phi_in)

# extract data from the mesh structure
Nr = u.domain.dims[1]
Ntheta = u.domain.dims[2]
G=reshape([1:(Nr+2)*(Ntheta+2);], Nr+2, Ntheta+2)
DRp = u.domain.cellsize.x[2:end-1]
DTHETAp = Array(Float64, 1, Ntheta)
DTHETAp[:] = u.domain.cellsize.y[2:end-1]
dr = 0.5*(u.domain.cellsize.x[1:end-1]+u.domain.cellsize.x[2:end])
dtheta = Array(Float64, 1, 1+Ntheta)
dtheta[:] = 0.5*(u.domain.cellsize.y[1:end-1]+u.domain.cellsize.y[2:end])
rp = u.domain.cellcenters.x
rf = u.domain.facecenters.x
psiX_p = zeros(Nr+1,Ntheta)
psiX_m = zeros(Nr+1,Ntheta)
psiY_p = zeros(Nr,Ntheta+1)
psiY_m = zeros(Nr,Ntheta+1)

# define the vectors to store the sparse matrix data
iix = zeros(Int64, 3*(Nr+2)*(Ntheta+2))
iiy = zeros(Int64, 3*(Nr+2)*(Ntheta+2))
jjx = zeros(Int64, 3*(Nr+2)*(Ntheta+2))
jjy = zeros(Int64, 3*(Nr+2)*(Ntheta+2))
sx = zeros(Float64, 3*(Nr+2)*(Ntheta+2))
sy = zeros(Float64, 3*(Nr+2)*(Ntheta+2))
mnx = Nr*Ntheta
mny = Nr*Ntheta

# calculate the upstream to downstream gradient ratios for u>0 (+ ratio)
# x direction
dphiX_p = (phi.value[2:Nr+2, 2:Ntheta+1]-phi.value[1:Nr+1, 2:Ntheta+1])./dr
rX_p = dphiX_p[1:end-1,:]./fsign(dphiX_p[2:end,:])
psiX_p[2:Nr+1,:] = 0.5*FL(rX_p).*(phi.value[3:Nr+2,2:Ntheta+1]-
		    phi.value[2:Nr+1, 2:Ntheta+1])
psiX_p[1, :] = 0.0 # left boundary will be handled in the main matrix
# y direction
dphiY_p = (phi.value[2:Nr+1, 2:Ntheta+2]-phi.value[2:Nr+1, 1:Ntheta+1])./dtheta
rY_p = dphiY_p[:,1:end-1]./fsign(dphiY_p[:,2:end])
psiY_p[:,2:Ntheta+1] = 0.5*FL(rY_p).*(phi.value[2:Nr+1,3:Ntheta+2]-
		  phi.value[2:Nr+1, 2:Ntheta+1])
psiY_p[:,1] = 0.0 # Bottom boundary will be handled in the main matrix

# calculate the upstream to downstream gradient ratios for u<0 (- ratio)
# x direction
rX_m = dphiX_p[2:end,:]./fsign(dphiX_p[1:end-1,:])
psiX_m[1:Nr,:] = 0.5*FL(rX_m).*(phi.value[1:Nr, 2:Ntheta+1]-
		phi.value[2:Nr+1, 2:Ntheta+1])
psiX_m[Nr+1,:] = 0.0 # right boundary
# y direction
rY_m = dphiY_p[:,2:end]./fsign(dphiY_p[:,1:end-1])
psiY_m[:,1:Ntheta] = 0.5*FL(rY_m).*(phi.value[2:Nr+1, 1:Ntheta]-
	      phi.value[2:Nr+1, 2:Ntheta+1])
psiY_m[:, Ntheta+1] = 0.0 # top boundary will be handled in the main matrix

re = rf[2:Nr+1,:]
rw = rf[1:Nr,:]

# find the velocity direction for the upwind scheme
ue_min = min(u.xvalue[2:Nr+1,:],0.0)
ue_max = max(u.xvalue[2:Nr+1,:],0.0)
uw_min = min(u.xvalue[1:Nr,:],0.0)
uw_max = max(u.xvalue[1:Nr,:],0.0)
vn_min = min(u.yvalue[:,2:Ntheta+1],0.0)
vn_max = max(u.yvalue[:,2:Ntheta+1],0.0)
vs_min = min(u.yvalue[:,1:Ntheta],0.0)
vs_max = max(u.yvalue[:,1:Ntheta],0.0)

# calculate the coefficients for the internal cells, do not reshape yet
AE = re.*ue_min./(DRp.*rp)
AW = -rw.*uw_max./(DRp.*rp)
AN = vn_min./(DTHETAp.*rp)
AS = -vs_max./(DTHETAp.*rp)
APx = (re.*ue_max-rw.*uw_min)./(DRp.*rp)
APy = (vn_max-vs_min)./(DTHETAp.*rp)

# Also correct for the boundary cells (not the ghost cells)
# Left boundary:
APx[1,:] = APx[1,:]-rw[1,:].*uw_max[1,:]./(2.0*DRp[1]*rp[1,:])
AW[1,:] = AW[1,:]/2.0
# Right boundary:
AE[end,:] = AE[end,:]/2.0
APx[end,:] = APx[end,:]+re[end,:].*ue_min[end,:]./(2.0*DRp[end]*rp[end,:])
# Bottom boundary:
APy[:,1] = APy[:,1]-vs_max[:,1]./(2.0*DTHETAp[1]*rp[:,1])
AS[:,1] = AS[:,1]/2.0
# Top boundary:
AN[:,end] = AN[:,end]/2.0
APy[:,end] = APy[:,end]+vn_min[:,end]./(2.0*DTHETAp[end]*rp[:,end])

# now reshape
AE = reshape(AE,mnx)
AW = reshape(AW,mnx)
AN = reshape(AN,mny)
AS = reshape(AS,mny)
APx = reshape(APx,mnx)
APy = reshape(APy,mny)

# build the sparse matrix based on the numbering system
rowx_index = reshape(G[2:Nr+1,2:Ntheta+1],mnx) # main diagonal x
iix[1:3*mnx] = repmat(rowx_index,3)
rowy_index = reshape(G[2:Nr+1,2:Ntheta+1],mny) # main diagonal y
iiy[1:3*mny] = repmat(rowy_index,3)
jjx[1:3*mnx] = [reshape(G[1:Nr,2:Ntheta+1],mnx); reshape(G[2:Nr+1,2:Ntheta+1],mnx); reshape(G[3:Nr+2,2:Ntheta+1],mnx)]
jjy[1:3*mny] = [reshape(G[2:Nr+1,1:Ntheta],mny); reshape(G[2:Nr+1,2:Ntheta+1],mny); reshape(G[2:Nr+1,3:Ntheta+2],mny)]
sx[1:3*mnx] = [AW; APx; AE]
sy[1:3*mny] = [AS; APy; AN]

# calculate the TVD correction term
div_x = -(1./(DRp.*rp)).*(re.*(ue_max.*psiX_p[2:Nr+1,:]+ue_min.*psiX_m[2:Nr+1,:])-
              rw.*(uw_max.*psiX_p[1:Nr,:]+uw_min.*psiX_m[1:Nr,:]))
div_y = -(1./(DTHETAp.*rp)).*((vn_max.*psiY_p[:,2:Ntheta+1]+vn_min.*psiY_m[:,2:Ntheta+1])-
              (vs_max.*psiY_p[:,1:Ntheta]+vs_min.*psiY_m[:,1:Ntheta]))

# define the RHS Vector
RHS = zeros((Nr+2)*(Ntheta+2))
RHSx = zeros((Nr+2)*(Ntheta+2))
RHSy = zeros((Nr+2)*(Ntheta+2))

# assign the values of the RHS vector
RHS[rowx_index] = reshape(div_x+div_y,Nr*Ntheta)
RHSx[rowx_index] = reshape(div_x,Nr*Ntheta)
RHSy[rowy_index] = reshape(div_y,Nr*Ntheta)

# build the sparse matrix
kx = 3*mnx
ky = 3*mny
Mx = sparse(iix[1:kx], jjx[1:kx], sx[1:kx], (Nr+2)*(Ntheta+2), (Nr+2)*(Ntheta+2))
My = sparse(iiy[1:ky], jjy[1:ky], sy[1:ky], (Nr+2)*(Ntheta+2), (Nr+2)*(Ntheta+2))
M = Mx + My
(M, RHS, Mx, My, RHSx, RHSy)
end

# ============================= 3D Convection Term ===============================
function convectionTerm3D(u::FaceValue)
# u is a face variable
# extract data from the mesh structure
Nx = u.domain.dims[1]
Ny = u.domain.dims[2]
Nz = u.domain.dims[3]
G=reshape([1:(Nx+2)*(Ny+2)*(Nz+2);], Nx+2, Ny+2, Nz+2)
DXe = u.domain.cellsize.x[3:end]
DXw = u.domain.cellsize.x[1:end-2]
DXp = u.domain.cellsize.x[2:end-1]
DY = Array(Float64, 1, Ny+2)
DY[:] = u.domain.cellsize.y
DYn = DY[1,3:end]
DYs = DY[1,1:end-2]
DYp = DY[1,2:end-1]
DZ = Array(Float64, 1, 1, Nz+2)
DZ[:] = u.domain.cellsize.z
DZf = DZ[1,1,3:end]
DZb = DZ[1,1,1:end-2]
DZp = DZ[1,1,2:end-1]

# define the vectors to stores the sparse matrix data
iix = zeros(Int64, 3*(Nx+2)*(Ny+2)*(Nz+2))
jjx = zeros(Int64, 3*(Nx+2)*(Ny+2)*(Nz+2))
sx = zeros(Float64, 3*(Nx+2)*(Ny+2)*(Nz+2))
iiy = zeros(Int64, 3*(Nx+2)*(Ny+2)*(Nz+2))
jjy = zeros(Int64, 3*(Nx+2)*(Ny+2)*(Nz+2))
sy = zeros(Float64, 3*(Nx+2)*(Ny+2)*(Nz+2))
iiz = zeros(Int64, 3*(Nx+2)*(Ny+2)*(Nz+2))
jjz = zeros(Int64, 3*(Nx+2)*(Ny+2)*(Nz+2))
sz = zeros(Float64, 3*(Nx+2)*(Ny+2)*(Nz+2))
mnx = Nx*Ny*Nz
mny = Nx*Ny*Nz
mnz = Nx*Ny*Nz

# reassign the east, west, north, and south velocity vectors for the
# code readability
ue = u.xvalue[2:Nx+1,:,:]./(DXp+DXe)
uw = u.xvalue[1:Nx,:,:]./(DXp+DXw)
vn = u.yvalue[:,2:Ny+1,:]./(DYp+DYn)
vs = u.yvalue[:,1:Ny,:]./(DYp+DYs)
wf = u.zvalue[:,:,2:Nz+1]./(DZp+DZf)
wb = u.zvalue[:,:,1:Nz]./(DZp+DZb)

# calculate the coefficients for the internal cells
AE = reshape(ue,mnx)
AW = reshape(-uw,mnx)
AN = reshape(vn,mny)
AS = reshape(-vs,mny)
AF = reshape(wf,mnz)
AB = reshape(-wb,mnz)
APx = reshape((ue.*DXe-uw.*DXw)./DXp,mnx)
APy = reshape((vn.*DYn-vs.*DYs)./DYp,mny)
APz = reshape((wf.*DZf-wb.*DZb)./DZp,mnz)

# build the sparse matrix based on the numbering system
rowx_index = reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],mnx)  # main diagonal x
iix[1:3*mnx] = repmat(rowx_index,3)
rowy_index = reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],mny)  # main diagonal y
iiy[1:3*mny] = repmat(rowy_index,3)
rowz_index = reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],mnz)  # main diagonal z
iiz[1:3*mnz] = repmat(rowz_index,3)
jjx[1:3*mnx] = [reshape(G[1:Nx,2:Ny+1,2:Nz+1],mnx); reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],mnx); reshape(G[3:Nx+2,2:Ny+1,2:Nz+1],mnx)]
jjy[1:3*mny] = [reshape(G[2:Nx+1,1:Ny,2:Nz+1],mny); reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],mny); reshape(G[2:Nx+1,3:Ny+2,2:Nz+1],mny)]
jjz[1:3*mnz] = [reshape(G[2:Nx+1,2:Ny+1,1:Nz],mnz); reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],mnz); reshape(G[2:Nx+1,2:Ny+1,3:Nz+2],mnz)]
sx[1:3*mnx] = [AW; APx; AE]
sy[1:3*mny] = [AS; APy; AN]
sz[1:3*mnz] = [AB; APz; AF]

# build the sparse matrix
kx = 3*mnx
ky = 3*mny
kz = 3*mnz
Mx = sparse(iix[1:kx], jjx[1:kx], sx[1:kx], (Nx+2)*(Ny+2)*(Nz+2), (Nx+2)*(Ny+2)*(Nz+2))
My = sparse(iiy[1:ky], jjy[1:ky], sy[1:ky], (Nx+2)*(Ny+2)*(Nz+2), (Nx+2)*(Ny+2)*(Nz+2))
Mz = sparse(iiz[1:kz], jjz[1:kz], sz[1:kz], (Nx+2)*(Ny+2)*(Nz+2), (Nx+2)*(Ny+2)*(Nz+2))
M = Mx + My + Mz

(M, Mx, My, Mz)
end

# ============================= 3D Convection Upwind Term ===============================
function convectionUpwindTerm3D(u::FaceValue)
# u is a face variable
# extract data from the mesh structure
Nx = u.domain.dims[1]
Ny = u.domain.dims[2]
Nz = u.domain.dims[3]
G=reshape([1:(Nx+2)*(Ny+2)*(Nz+2);], Nx+2, Ny+2, Nz+2)
DXp = u.domain.cellsize.x[2:end-1]
DY = Array(Float64, 1, Ny+2)
DY[:] = u.domain.cellsize.y
DYp = DY[1,2:end-1]
DZ = Array(Float64, 1, 1, Nz+2)
DZ[:] = u.domain.cellsize.z
DZp = DZ[1,1,2:end-1]

# define the vectors to stores the sparse matrix data
iix = zeros(Int64, 3*(Nx+2)*(Ny+2)*(Nz+2))
jjx = zeros(Int64, 3*(Nx+2)*(Ny+2)*(Nz+2))
sx = zeros(Float64, 3*(Nx+2)*(Ny+2)*(Nz+2))
iiy = zeros(Int64, 3*(Nx+2)*(Ny+2)*(Nz+2))
jjy = zeros(Int64, 3*(Nx+2)*(Ny+2)*(Nz+2))
sy = zeros(Float64, 3*(Nx+2)*(Ny+2)*(Nz+2))
iiz = zeros(Int64, 3*(Nx+2)*(Ny+2)*(Nz+2))
jjz = zeros(Int64, 3*(Nx+2)*(Ny+2)*(Nz+2))
sz = zeros(Float64, 3*(Nx+2)*(Ny+2)*(Nz+2))
mnx = Nx*Ny*Nz
mny = Nx*Ny*Nz
mnz = Nx*Ny*Nz

# find the velocity direction for the upwind scheme
ue_min = min(u.xvalue[2:Nx+1,:,:],0)
ue_max = max(u.xvalue[2:Nx+1,:,:],0)
uw_min = min(u.xvalue[1:Nx,:,:],0)
uw_max = max(u.xvalue[1:Nx,:,:],0)
vn_min = min(u.yvalue[:,2:Ny+1,:],0)
vn_max = max(u.yvalue[:,2:Ny+1,:],0)
vs_min = min(u.yvalue[:,1:Ny,:],0)
vs_max = max(u.yvalue[:,1:Ny,:],0)
wf_min = min(u.zvalue[:,:,2:Nz+1],0)
wf_max = max(u.zvalue[:,:,2:Nz+1],0)
wb_min = min(u.zvalue[:,:,1:Nz],0)
wb_max = max(u.zvalue[:,:,1:Nz],0)

# calculate the coefficients for the internal cells
AE = ue_min./DXp
AW = -uw_max./DXp
AN = vn_min./DYp
AS = -vs_max./DYp
AF = wf_min./DZp
AB = -wb_max./DZp
APx = (ue_max-uw_min)./DXp
APy = (vn_max-vs_min)./DYp
APz = (wf_max-wb_min)./DZp

# Also correct for the boundary cells (not the ghost cells)
# Left boundary:
APx[1,:,:] = APx[1,:,:]-uw_max[1,:,:]/(2.0*DXp[1])
AW[1,:,:] = AW[1,:,:]/2.0
# Right boundary:
AE[end,:,:] = AE[end,:,:]/2.0
APx[end,:,:] = APx[end,:,:]+ue_min[end,:,:]/(2.0*DXp[end])
# Bottom boundary:
APy[:,1,:] = APy[:,1,:]-vs_max[:,1,:]/(2.0*DYp[1])
AS[:,1,:] = AS[:,1,:]/2.0
# Top boundary:
AN[:,end,:] = AN[:,end,:]/2.0
APy[:,end,:] = APy[:,end,:]+vn_min[:,end,:]/(2.0*DYp[end])
# Back boundary:
APz[:,:,1] = APz[:,:,1]-wb_max[:,:,1]/(2.0*DZp[1])
AB[:,:,1] = AB[:,:,1]/2.0
# Front boundary:
AF[:,:,end] = AF[:,:,end]/2.0
APz[:,:,end] = APz[:,:,end]+wf_min[:,:,end]/(2.0*DZp[end])

AE = reshape(AE,mnx)
AW = reshape(AW,mnx)
AN = reshape(AN,mny)
AS = reshape(AS,mny)
AF = reshape(AF,mnz)
AB = reshape(AB,mnz)
APx = reshape(APx,mnx)
APy = reshape(APy,mny)
APz = reshape(APz,mnz)

# build the sparse matrix based on the numbering system
rowx_index = reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],mnx)  # main diagonal x
iix[1:3*mnx] = repmat(rowx_index,3)
rowy_index = reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],mny)  # main diagonal y
iiy[1:3*mny] = repmat(rowy_index,3)
rowz_index = reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],mnz)  # main diagonal z
iiz[1:3*mnz] = repmat(rowz_index,3)
jjx[1:3*mnx] = [reshape(G[1:Nx,2:Ny+1,2:Nz+1],mnx); reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],mnx); reshape(G[3:Nx+2,2:Ny+1,2:Nz+1],mnx)]
jjy[1:3*mny] = [reshape(G[2:Nx+1,1:Ny,2:Nz+1],mny); reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],mny); reshape(G[2:Nx+1,3:Ny+2,2:Nz+1],mny)]
jjz[1:3*mnz] = [reshape(G[2:Nx+1,2:Ny+1,1:Nz],mnz); reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],mnz); reshape(G[2:Nx+1,2:Ny+1,3:Nz+2],mnz)]
sx[1:3*mnx] = [AW; APx; AE]
sy[1:3*mny] = [AS; APy; AN]
sz[1:3*mnz] = [AB; APz; AF]

# build the sparse matrix
kx = 3*mnx
ky = 3*mny
kz = 3*mnz
Mx = sparse(iix[1:kx], jjx[1:kx], sx[1:kx], (Nx+2)*(Ny+2)*(Nz+2), (Nx+2)*(Ny+2)*(Nz+2))
My = sparse(iiy[1:ky], jjy[1:ky], sy[1:ky], (Nx+2)*(Ny+2)*(Nz+2), (Nx+2)*(Ny+2)*(Nz+2))
Mz = sparse(iiz[1:kz], jjz[1:kz], sz[1:kz], (Nx+2)*(Ny+2)*(Nz+2), (Nx+2)*(Ny+2)*(Nz+2))
M = Mx + My + Mz;

(M, Mx, My, Mz)
end

# ============================= 3D Convection TVD Term ===============================
function convectionTvdTerm3D(u::FaceValue, phi::CellValue, FL::Function)
# u is a face variable
# a function to avoid division by zero
eps1 = 1.0e-20
fsign(phi_in) = (abs(phi_in).>=eps1).*phi_in+eps1*(phi_in.==0.0)+eps1*(abs(phi_in).<eps1).*sign(phi_in)
# extract data from the mesh structure
Nx = u.domain.dims[1]
Ny = u.domain.dims[2]
Nz = u.domain.dims[3]
G=reshape([1:(Nx+2)*(Ny+2)*(Nz+2);], Nx+2, Ny+2, Nz+2)
DXp = u.domain.cellsize.x[2:end-1]
DY = Array(Float64, 1, Ny+2)
DY[:] = u.domain.cellsize.y
DYp = DY[1,2:end-1]
DZ = Array(Float64, 1, 1, Nz+2)
DZ[:] = u.domain.cellsize.z
DZp = DZ[1,1,2:end-1]
dx=0.5*(u.domain.cellsize.x[1:end-1]+u.domain.cellsize.x[2:end])
dy=Array(Float64, 1, Ny+1)
dy[:]=0.5*(u.domain.cellsize.y[1:end-1]+u.domain.cellsize.y[2:end])
dz=Array(Float64, 1, 1, Nz+1)
dz[:]=0.5*(u.domain.cellsize.z[1:end-1]+u.domain.cellsize.z[2:end])
psiX_p = zeros(Nx+1,Ny,Nz)
psiX_m = zeros(Nx+1,Ny,Nz)
psiY_p = zeros(Nx,Ny+1,Nz)
psiY_m = zeros(Nx,Ny+1,Nz)
psiZ_p = zeros(Nx,Ny,Nz+1)
psiZ_m = zeros(Nx,Ny,Nz+1)

# define the vectors to stores the sparse matrix data
iix = zeros(Int64, 3*(Nx+2)*(Ny+2)*(Nz+2))
jjx = zeros(Int64, 3*(Nx+2)*(Ny+2)*(Nz+2))
sx = zeros(Float64, 3*(Nx+2)*(Ny+2)*(Nz+2))
iiy = zeros(Int64, 3*(Nx+2)*(Ny+2)*(Nz+2))
jjy = zeros(Int64, 3*(Nx+2)*(Ny+2)*(Nz+2))
sy = zeros(Float64, 3*(Nx+2)*(Ny+2)*(Nz+2))
iiz = zeros(Int64, 3*(Nx+2)*(Ny+2)*(Nz+2))
jjz = zeros(Int64, 3*(Nx+2)*(Ny+2)*(Nz+2))
sz = zeros(Float64, 3*(Nx+2)*(Ny+2)*(Nz+2))
mnx = Nx*Ny*Nz
mny = Nx*Ny*Nz
mnz = Nx*Ny*Nz

# calculate the upstream to downstream gradient ratios for u>0 (+ ratio)
# x direction
dphiX_p = (phi.value[2:Nx+2, 2:Ny+1, 2:Nz+1]-phi.value[1:Nx+1, 2:Ny+1, 2:Nz+1])./dx
rX_p = dphiX_p[1:end-1,:,:]./fsign(dphiX_p[2:end,:,:])
psiX_p[2:Nx+1,:,:] = 0.5*FL(rX_p).*(phi.value[3:Nx+2,2:Ny+1,2:Nz+1]-phi.value[2:Nx+1,2:Ny+1,2:Nz+1])
psiX_p[1,:,:] = 0.0  # left boundary
# y direction
dphiY_p = (phi.value[2:Nx+1, 2:Ny+2, 2:Nz+1]-phi.value[2:Nx+1, 1:Ny+1, 2:Nz+1])./dy
rY_p = dphiY_p[:,1:end-1,:]./fsign(dphiY_p[:,2:end,:])
psiY_p[:,2:Ny+1,:] = 0.5*FL(rY_p).*(phi.value[2:Nx+1,3:Ny+2,2:Nz+1]-phi.value[2:Nx+1, 2:Ny+1,2:Nz+1])
psiY_p[:,1,:] = 0.0  # Bottom boundary
# z direction
dphiZ_p = (phi.value[2:Nx+1, 2:Ny+1, 2:Nz+2]-phi.value[2:Nx+1, 2:Ny+1, 1:Nz+1])./dz
rZ_p = dphiZ_p[:,:,1:end-1]./fsign(dphiZ_p[:,:,2:end])
psiZ_p[:,:,2:Nz+1] = 0.5*FL(rZ_p).*(phi.value[2:Nx+1,2:Ny+1,3:Nz+2]-phi.value[2:Nx+1,2:Ny+1,2:Nz+1])
psiZ_p[:,:,1] = 0.0  # Back boundary

# calculate the upstream to downstream gradient ratios for u<0 (- ratio)
# x direction
rX_m = dphiX_p[2:end,:,:]./fsign(dphiX_p[1:end-1,:,:])
psiX_m[1:Nx,:,:] = 0.5*FL(rX_m).*(phi.value[1:Nx, 2:Ny+1, 2:Nz+1]-phi.value[2:Nx+1, 2:Ny+1, 2:Nz+1])
psiX_m[Nx+1,:,:] = 0.0  # right boundary
# y direction
rY_m = dphiY_p[:,2:end,:]./fsign(dphiY_p[:,1:end-1,:])
psiY_m[:,1:Ny,:] = 0.5*FL(rY_m).*(phi.value[2:Nx+1,1:Ny,2:Nz+1]-phi.value[2:Nx+1,2:Ny+1,2:Nz+1])
psiY_m[:,Ny+1,:] = 0.0  # top boundary
# z direction
rZ_m = dphiZ_p[:,:,2:end]./fsign(dphiZ_p[:,:,1:end-1])
psiZ_m[:,:,1:Nz] = 0.5*FL(rZ_m).*(phi.value[2:Nx+1,2:Ny+1,1:Nz]-phi.value[2:Nx+1,2:Ny+1,2:Nz+1])
psiZ_m[:,:,Nz+1] = 0.0  # front boundary

# find the velocity direction for the upwind scheme
ue_min = min(u.xvalue[2:Nx+1,:,:],0)
ue_max = max(u.xvalue[2:Nx+1,:,:],0)
uw_min = min(u.xvalue[1:Nx,:,:],0)
uw_max = max(u.xvalue[1:Nx,:,:],0)
vn_min = min(u.yvalue[:,2:Ny+1,:],0)
vn_max = max(u.yvalue[:,2:Ny+1,:],0)
vs_min = min(u.yvalue[:,1:Ny,:],0)
vs_max = max(u.yvalue[:,1:Ny,:],0)
wf_min = min(u.zvalue[:,:,2:Nz+1],0)
wf_max = max(u.zvalue[:,:,2:Nz+1],0)
wb_min = min(u.zvalue[:,:,1:Nz],0)
wb_max = max(u.zvalue[:,:,1:Nz],0)

# calculate the coefficients for the internal cells
AE = ue_min./DXp
AW = -uw_max./DXp
AN = vn_min./DYp
AS = -vs_max./DYp
AF = wf_min./DZp
AB = -wb_max./DZp
APx = (ue_max-uw_min)./DXp
APy = (vn_max-vs_min)./DYp
APz = (wf_max-wb_min)./DZp

# Also correct for the boundary cells (not the ghost cells)
# Left boundary:
APx[1,:,:] = APx[1,:,:]-uw_max[1,:,:]/(2.0*DXp[1])
AW[1,:,:] = AW[1,:,:]/2.0
# Right boundary:
AE[end,:,:] = AE[end,:,:]/2.0
APx[end,:,:] = APx[end,:,:]+ue_min[end,:,:]/(2.0*DXp[end])
# Bottom boundary:
APy[:,1,:] = APy[:,1,:]-vs_max[:,1,:]/(2.0*DYp[1])
AS[:,1,:] = AS[:,1,:]/2.0
# Top boundary:
AN[:,end,:] = AN[:,end,:]/2.0
APy[:,end,:] = APy[:,end,:]+vn_min[:,end,:]/(2.0*DYp[end])
# Back boundary:
APz[:,:,1] = APz[:,:,1]-wb_max[:,:,1]/(2.0*DZp[1])
AB[:,:,1] = AB[:,:,1]/2.0
# Front boundary:
AF[:,:,end] = AF[:,:,end]/2.0
APz[:,:,end] = APz[:,:,end]+wf_min[:,:,end]/(2.0*DZp[end])

AE = reshape(AE,mnx)
AW = reshape(AW,mnx)
AN = reshape(AN,mny)
AS = reshape(AS,mny)
AF = reshape(AF,mnz)
AB = reshape(AB,mnz)
APx = reshape(APx,mnx)
APy = reshape(APy,mny)
APz = reshape(APz,mnz)

# build the sparse matrix based on the numbering system
rowx_index = reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],mnx)  # main diagonal x
iix[1:3*mnx] = repmat(rowx_index,3)
rowy_index = reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],mny)  # main diagonal y
iiy[1:3*mny] = repmat(rowy_index,3)
rowz_index = reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],mnz)  # main diagonal z
iiz[1:3*mnz] = repmat(rowz_index,3)
jjx[1:3*mnx] = [reshape(G[1:Nx,2:Ny+1,2:Nz+1],mnx); reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],mnx); reshape(G[3:Nx+2,2:Ny+1,2:Nz+1],mnx)]
jjy[1:3*mny] = [reshape(G[2:Nx+1,1:Ny,2:Nz+1],mny); reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],mny); reshape(G[2:Nx+1,3:Ny+2,2:Nz+1],mny)]
jjz[1:3*mnz] = [reshape(G[2:Nx+1,2:Ny+1,1:Nz],mnz); reshape(G[2:Nx+1,2:Ny+1,2:Nz+1],mnz); reshape(G[2:Nx+1,2:Ny+1,3:Nz+2],mnz)]
sx[1:3*mnx] = [AW; APx; AE]
sy[1:3*mny] = [AS; APy; AN]
sz[1:3*mnz] = [AB; APz; AF]

# calculate the TVD correction term
div_x = -(1./DXp).*((ue_max.*psiX_p[2:Nx+1,:,:]+ue_min.*psiX_m[2:Nx+1,:,:])-
              (uw_max.*psiX_p[1:Nx,:,:]+uw_min.*psiX_m[1:Nx,:,:]))
div_y = -(1./DYp).*((vn_max.*psiY_p[:,2:Ny+1,:]+vn_min.*psiY_m[:,2:Ny+1,:])-
              (vs_max.*psiY_p[:,1:Ny,:]+vs_min.*psiY_m[:,1:Ny,:]))
div_z = -(1./DZp).*((wf_max.*psiZ_p[:,:,2:Nz+1]+wf_min.*psiZ_m[:,:,2:Nz+1])-
              (wb_max.*psiZ_p[:,:,1:Nz]+wb_min.*psiZ_m[:,:,1:Nz]))

# define the RHS Vector
RHS = zeros((Nx+2)*(Ny+2)*(Nz+2))
RHSx = zeros((Nx+2)*(Ny+2)*(Nz+2))
RHSy = zeros((Nx+2)*(Ny+2)*(Nz+2))
RHSz = zeros((Nx+2)*(Ny+2)*(Nz+2))

# assign the values of the RHS vector
row_index = rowx_index
RHS[row_index] = reshape(div_x+div_y+div_z,Nx*Ny*Nz)
RHSx[rowx_index] = reshape(div_x,Nx*Ny*Nz)
RHSy[rowy_index] = reshape(div_y,Nx*Ny*Nz)
RHSz[rowz_index] = reshape(div_z,Nx*Ny*Nz)

# build the sparse matrix
kx = 3*mnx
ky = 3*mny
kz = 3*mnz
Mx = sparse(iix[1:kx], jjx[1:kx], sx[1:kx], (Nx+2)*(Ny+2)*(Nz+2), (Nx+2)*(Ny+2)*(Nz+2))
My = sparse(iiy[1:ky], jjy[1:ky], sy[1:ky], (Nx+2)*(Ny+2)*(Nz+2), (Nx+2)*(Ny+2)*(Nz+2))
Mz = sparse(iiz[1:kz], jjz[1:kz], sz[1:kz], (Nx+2)*(Ny+2)*(Nz+2), (Nx+2)*(Ny+2)*(Nz+2))
M = Mx + My + Mz;

(M, RHS, Mx, My, Mz, RHSx, RHSy, RHSz)
end

# ============================= 3D Cylindrical Convection Term ===============================
function convectionTermCylindrical3D(u::FaceValue)
# u is a face variable
# extract data from the mesh structure
Nr = u.domain.dims[1]
Ntheta = u.domain.dims[2]
Nz = u.domain.dims[3]
G=reshape([1:(Nr+2)*(Ntheta+2)*(Nz+2);], Nr+2, Ntheta+2, Nz+2)
DRe = u.domain.cellsize.x[3:end]
DRw = u.domain.cellsize.x[1:end-2]
DRp = u.domain.cellsize.x[2:end-1]
DTHETA = Array(Float64, 1, Ntheta+2)
DTHETA[:] = u.domain.cellsize.y
DTHETAn = DTHETA[1,3:end]
DTHETAs = DTHETA[1,1:end-2]
DTHETAp = DTHETA[1,2:end-1]
DZ = Array(Float64, 1, 1, Nz+2)
DZ[:] = u.domain.cellsize.z
DZf = DZ[1,1,3:end]
DZb = DZ[1,1,1:end-2]
DZp = DZ[1,1,2:end-1]
rp = u.domain.cellcenters.x
rf = u.domain.facecenters.x

# define the vectors to stores the sparse matrix data
iix = zeros(Int64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
jjx = zeros(Int64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
sx = zeros(Float64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
iiy = zeros(Int64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
jjy = zeros(Int64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
sy = zeros(Float64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
iiz = zeros(Int64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
jjz = zeros(Int64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
sz = zeros(Float64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
mnx = Nr*Ntheta*Nz
mny = Nr*Ntheta*Nz
mnz = Nr*Ntheta*Nz

re = rf[2:Nr+1]
rw = rf[1:Nr]

# reassign the east, west, north, and south velocity vectors for the
# code readability
ue = rf[2:Nr+1].*u.xvalue[2:Nr+1,:,:]./(rp.*(DRp+DRe))
uw = rf[1:Nr].*u.xvalue[1:Nr,:,:]./(rp.*(DRp+DRw))
vn = u.yvalue[:,2:Ntheta+1,:]./(rp.*(DTHETAp+DTHETAn))
vs = u.yvalue[:,1:Ntheta,:]./(rp.*(DTHETAp+DTHETAs))
wf = u.zvalue[:,:,2:Nz+1]./(DZp+DZf)
wb = u.zvalue[:,:,1:Nz]./(DZp+DZb)

# calculate the coefficients for the internal cells
AE = reshape(ue,mnx)
AW = reshape(-uw,mnx)
AN = reshape(vn,mny)
AS = reshape(-vs,mny)
AF = reshape(wf,mnz)
AB = reshape(-wb,mnz)
APx = reshape((DRe.*ue-DRw.*uw)./DRp,mnx)
APy = reshape((DTHETAn.*vn-DTHETAs.*vs)./DTHETAp,mny)
APz = reshape((DZf.*wf-DZb.*wb)./DZp,mnz)

# build the sparse matrix based on the numbering system
rowx_index = reshape(G[2:Nr+1,2:Ntheta+1,2:Nz+1],mnx)  # main diagonal x
iix[1:3*mnx] = repmat(rowx_index,3)
rowy_index = reshape(G[2:Nr+1,2:Ntheta+1,2:Nz+1],mny)  # main diagonal y
iiy[1:3*mny] = repmat(rowy_index,3)
rowz_index = reshape(G[2:Nr+1,2:Ntheta+1,2:Nz+1],mnz)  # main diagonal z
iiz[1:3*mnz] = repmat(rowz_index,3)
jjx[1:3*mnx] = [reshape(G[1:Nr,2:Ntheta+1,2:Nz+1],mnx); reshape(G[2:Nr+1,2:Ntheta+1,2:Nz+1],mnx); reshape(G[3:Nr+2,2:Ntheta+1,2:Nz+1],mnx)]
jjy[1:3*mny] = [reshape(G[2:Nr+1,1:Ntheta,2:Nz+1],mny); reshape(G[2:Nr+1,2:Ntheta+1,2:Nz+1],mny); reshape(G[2:Nr+1,3:Ntheta+2,2:Nz+1],mny)]
jjz[1:3*mnz] = [reshape(G[2:Nr+1,2:Ntheta+1,1:Nz],mnz); reshape(G[2:Nr+1,2:Ntheta+1,2:Nz+1],mnz); reshape(G[2:Nr+1,2:Ntheta+1,3:Nz+2],mnz)]
sx[1:3*mnx] = [AW; APx; AE]
sy[1:3*mny] = [AS; APy; AN]
sz[1:3*mnz] = [AB; APz; AF]

# build the sparse matrix
kx = 3*mnx
ky = 3*mny
kz = 3*mnz
Mx = sparse(iix[1:kx], jjx[1:kx], sx[1:kx], (Nr+2)*(Ntheta+2)*(Nz+2), (Nr+2)*(Ntheta+2)*(Nz+2))
My = sparse(iiy[1:ky], jjy[1:ky], sy[1:ky], (Nr+2)*(Ntheta+2)*(Nz+2), (Nr+2)*(Ntheta+2)*(Nz+2))
Mz = sparse(iiz[1:kz], jjz[1:kz], sz[1:kz], (Nr+2)*(Ntheta+2)*(Nz+2), (Nr+2)*(Ntheta+2)*(Nz+2))
M = Mx + My + Mz

(M, Mx, My, Mz)
end

# ============================= 3D Cylindrical Convection Upwind Term ===============================
function convectionUpwindTermCylindrical3D(u::FaceValue)
# u is a face variable
# extract data from the mesh structure
Nr = u.domain.dims[1]
Ntheta = u.domain.dims[2]
Nz = u.domain.dims[3]
G=reshape([1:(Nr+2)*(Ntheta+2)*(Nz+2);], Nr+2, Ntheta+2, Nz+2)
DRe = u.domain.cellsize.x[3:end]
DRw = u.domain.cellsize.x[1:end-2]
DRp = u.domain.cellsize.x[2:end-1]
DTHETA = Array(Float64, 1, Ntheta+2)
DTHETA[:] = u.domain.cellsize.y
DTHETAn = DTHETA[1,3:end]
DTHETAs = DTHETA[1,1:end-2]
DTHETAp = DTHETA[1,2:end-1]
DZ = Array(Float64, 1, 1, Nz+2)
DZ[:] = u.domain.cellsize.z
DZf = DZ[1,1,3:end]
DZb = DZ[1,1,1:end-2]
DZp = DZ[1,1,2:end-1]
rp = u.domain.cellcenters.x
rf = u.domain.facecenters.x

# define the vectors to stores the sparse matrix data
iix = zeros(Int64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
jjx = zeros(Int64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
sx = zeros(Float64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
iiy = zeros(Int64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
jjy = zeros(Int64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
sy = zeros(Float64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
iiz = zeros(Int64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
jjz = zeros(Int64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
sz = zeros(Float64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
mnx = Nr*Ntheta*Nz
mny = Nr*Ntheta*Nz
mnz = Nr*Ntheta*Nz

re = rf[2:Nr+1,:]
rw = rf[1:Nr,:]

# find the velocity direction for the upwind scheme
ue_min = min(u.xvalue[2:Nr+1,:,:],0)
ue_max = max(u.xvalue[2:Nr+1,:,:],0)
uw_min = min(u.xvalue[1:Nr,:,:],0)
uw_max = max(u.xvalue[1:Nr,:,:],0)
vn_min = min(u.yvalue[:,2:Ntheta+1,:],0)
vn_max = max(u.yvalue[:,2:Ntheta+1,:],0)
vs_min = min(u.yvalue[:,1:Ntheta,:],0)
vs_max = max(u.yvalue[:,1:Ntheta,:],0)
wf_min = min(u.zvalue[:,:,2:Nz+1],0)
wf_max = max(u.zvalue[:,:,2:Nz+1],0)
wb_min = min(u.zvalue[:,:,1:Nz],0)
wb_max = max(u.zvalue[:,:,1:Nz],0)

# calculate the coefficients for the internal cells
AE = re.*ue_min./(DRp.*rp)
AW = -rw.*uw_max./(DRp.*rp)
AN = vn_min./(DTHETAp.*rp)
AS = -vs_max./(DTHETAp.*rp)
AF = wf_min./DZp
AB = -wb_max./DZp
APx = (re.*ue_max-rw.*uw_min)./(DRp.*rp)
APy = (vn_max-vs_min)./(DTHETAp.*rp)
APz = (wf_max-wb_min)./DZp

# Also correct for the boundary cells (not the ghost cells)
# Left boundary:
APx[1,:,:] = APx[1,:,:]-rw[1,:,:].*uw_max[1,:,:]./(2.0*DRp[1]*rp[1,:,:])
AW[1,:,:] = AW[1,:,:]/2.0
# Right boundary:
AE[end,:,:] = AE[end,:,:]/2.0
APx[end,:,:] = APx[end,:,:]+re[end,:,:].*ue_min[end,:,:]./(2.0*DRp[end]*rp[end,:,:])
# Bottom boundary:
APy[:,1,:] = APy[:,1,:]-vs_max[:,1,:]./(2.0*DTHETAp[1]*rp[:,1,:])
AS[:,1,:] = AS[:,1,:]/2.0
# Top boundary:
AN[:,end,:] = AN[:,end,:]/2.0
APy[:,end,:] = APy[:,end,:]+vn_min[:,end,:]./(2.0*DTHETAp[end]*rp[:,end,:])
# Back boundary:
APz[:,:,1] = APz[:,:,1]-wb_max[:,:,1]/(2.0*DZp[1])
AB[:,:,1] = AB[:,:,1]/2.0
# Front boundary:
AF[:,:,end] = AF[:,:,end]/2.0
APz[:,:,end] = APz[:,:,end]+wf_min[:,:,end]/(2.0*DZp[end])

AE = reshape(AE,mnx)
AW = reshape(AW,mnx)
AN = reshape(AN,mny)
AS = reshape(AS,mny)
AF = reshape(AF,mnz)
AB = reshape(AB,mnz)
APx = reshape(APx,mnx)
APy = reshape(APy,mny)
APz = reshape(APz,mnz)

# build the sparse matrix based on the numbering system
rowx_index = reshape(G[2:Nr+1,2:Ntheta+1,2:Nz+1],mnx)  # main diagonal x
iix[1:3*mnx] = repmat(rowx_index,3)
rowy_index = reshape(G[2:Nr+1,2:Ntheta+1,2:Nz+1],mny)  # main diagonal y
iiy[1:3*mny] = repmat(rowy_index,3)
rowz_index = reshape(G[2:Nr+1,2:Ntheta+1,2:Nz+1],mnz)  # main diagonal z
iiz[1:3*mnz] = repmat(rowz_index,3)
jjx[1:3*mnx] = [reshape(G[1:Nr,2:Ntheta+1,2:Nz+1],mnx); reshape(G[2:Nr+1,2:Ntheta+1,2:Nz+1],mnx); reshape(G[3:Nr+2,2:Ntheta+1,2:Nz+1],mnx)]
jjy[1:3*mny] = [reshape(G[2:Nr+1,1:Ntheta,2:Nz+1],mny); reshape(G[2:Nr+1,2:Ntheta+1,2:Nz+1],mny); reshape(G[2:Nr+1,3:Ntheta+2,2:Nz+1],mny)]
jjz[1:3*mnz] = [reshape(G[2:Nr+1,2:Ntheta+1,1:Nz],mnz); reshape(G[2:Nr+1,2:Ntheta+1,2:Nz+1],mnz); reshape(G[2:Nr+1,2:Ntheta+1,3:Nz+2],mnz)]
sx[1:3*mnx] = [AW; APx; AE]
sy[1:3*mny] = [AS; APy; AN]
sz[1:3*mnz] = [AB; APz; AF]

# build the sparse matrix
kx = 3*mnx
ky = 3*mny
kz = 3*mnz
Mx = sparse(iix[1:kx], jjx[1:kx], sx[1:kx], (Nr+2)*(Ntheta+2)*(Nz+2), (Nr+2)*(Ntheta+2)*(Nz+2))
My = sparse(iiy[1:ky], jjy[1:ky], sy[1:ky], (Nr+2)*(Ntheta+2)*(Nz+2), (Nr+2)*(Ntheta+2)*(Nz+2))
Mz = sparse(iiz[1:kz], jjz[1:kz], sz[1:kz], (Nr+2)*(Ntheta+2)*(Nz+2), (Nr+2)*(Ntheta+2)*(Nz+2))
M = Mx + My + Mz;

(M, Mx, My, Mz)
end

# ============================= 3D Cylindrical Convection TVD Term ===============================
function convectionTvdTermCylindrical3D(u::FaceValue, phi::CellValue, FL::Function)
# u is a face variable
# a function to avoid division by zero
eps1 = 1.0e-20
fsign(phi_in) = (abs(phi_in).>=eps1).*phi_in+eps1*(phi_in.==0.0)+eps1*(abs(phi_in).<eps1).*sign(phi_in)
# extract data from the mesh structure
Nr = u.domain.dims[1]
Ntheta = u.domain.dims[2]
Nz = u.domain.dims[3]
G=reshape([1:(Nr+2)*(Ntheta+2)*(Nz+2);], Nr+2, Ntheta+2, Nz+2)
DRe = u.domain.cellsize.x[3:end]
DRw = u.domain.cellsize.x[1:end-2]
DRp = u.domain.cellsize.x[2:end-1]
DTHETA = Array(Float64, 1, Ntheta+2)
DTHETA[:] = u.domain.cellsize.y
DTHETAn = DTHETA[1,3:end]
DTHETAs = DTHETA[1,1:end-2]
DTHETAp = DTHETA[1,2:end-1]
DZ = Array(Float64, 1, 1, Nz+2)
DZ[:] = u.domain.cellsize.z
DZf = DZ[1,1,3:end]
DZb = DZ[1,1,1:end-2]
DZp = DZ[1,1,2:end-1]
dr=0.5*(u.domain.cellsize.x[1:end-1]+u.domain.cellsize.x[2:end])
dtheta = Array(Float64, 1, Ntheta+1)
dtheta[:]=0.5*(u.domain.cellsize.y[1:end-1]+u.domain.cellsize.y[2:end])
dz = Array(Float64, 1, 1, Nz+1)
dz[:]=0.5*(u.domain.cellsize.z[1:end-1]+u.domain.cellsize.z[2:end])
psiX_p = zeros(Nr+1,Ntheta,Nz)
psiX_m = zeros(Nr+1,Ntheta,Nz)
psiY_p = zeros(Nr,Ntheta+1,Nz)
psiY_m = zeros(Nr,Ntheta+1,Nz)
psiZ_p = zeros(Nr,Ntheta,Nz+1)
psiZ_m = zeros(Nr,Ntheta,Nz+1)
rp = u.domain.cellcenters.x
rf = u.domain.facecenters.x

# define the vectors to stores the sparse matrix data
iix = zeros(Int64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
jjx = zeros(Int64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
sx = zeros(Float64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
iiy = zeros(Int64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
jjy = zeros(Int64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
sy = zeros(Float64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
iiz = zeros(Int64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
jjz = zeros(Int64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
sz = zeros(Float64, 3*(Nr+2)*(Ntheta+2)*(Nz+2))
mnx = Nr*Ntheta*Nz
mny = Nr*Ntheta*Nz
mnz = Nr*Ntheta*Nz

# calculate the upstream to downstream gradient ratios for u>0 (+ ratio)
# x direction
dphiX_p = (phi.value[2:Nr+2, 2:Ntheta+1, 2:Nz+1]-phi.value[1:Nr+1, 2:Ntheta+1, 2:Nz+1])./dr
rX_p = dphiX_p[1:end-1,:,:]./fsign(dphiX_p[2:end,:,:])
psiX_p[2:Nr+1,:,:] = 0.5*FL(rX_p).*(phi.value[3:Nr+2,2:Ntheta+1,2:Nz+1]-phi.value[2:Nr+1,2:Ntheta+1,2:Nz+1])
psiX_p[1,:,:] = 0  # left boundary
# y direction
dphiY_p = (phi.value[2:Nr+1, 2:Ntheta+2, 2:Nz+1]-phi.value[2:Nr+1, 1:Ntheta+1, 2:Nz+1])./dtheta
rY_p = dphiY_p[:,1:end-1,:]./fsign(dphiY_p[:,2:end,:])
psiY_p[:,2:Ntheta+1,:] = 0.5*FL(rY_p).*(phi.value[2:Nr+1,3:Ntheta+2,2:Nz+1]-phi.value[2:Nr+1, 2:Ntheta+1,2:Nz+1])
psiY_p[:,1,:] = 0.0  # Bottom boundary
# z direction
dphiZ_p = (phi.value[2:Nr+1, 2:Ntheta+1, 2:Nz+2]-phi.value[2:Nr+1, 2:Ntheta+1, 1:Nz+1])./dz
rZ_p = dphiZ_p[:,:,1:end-1]./fsign(dphiZ_p[:,:,2:end])
psiZ_p[:,:,2:Nz+1] = 0.5*FL(rZ_p).*(phi.value[2:Nr+1,2:Ntheta+1,3:Nz+2]-phi.value[2:Nr+1,2:Ntheta+1,2:Nz+1])
psiZ_p[:,:,1] = 0.0  # Back boundary

# calculate the upstream to downstream gradient ratios for u<0 (- ratio)
# x direction
rX_m = dphiX_p[2:end,:,:]./fsign(dphiX_p[1:end-1,:,:])
psiX_m[1:Nr,:,:] = 0.5*FL(rX_m).*(phi.value[1:Nr, 2:Ntheta+1, 2:Nz+1]-phi.value[2:Nr+1, 2:Ntheta+1, 2:Nz+1])
psiX_m[Nr+1,:,:] = 0.0  # right boundary
# y direction
rY_m = dphiY_p[:,2:end,:]./fsign(dphiY_p[:,1:end-1,:])
psiY_m[:,1:Ntheta,:] = 0.5*FL(rY_m).*(phi.value[2:Nr+1,1:Ntheta,2:Nz+1]-phi.value[2:Nr+1,2:Ntheta+1,2:Nz+1])
psiY_m[:,Ntheta+1,:] = 0.0  # top boundary
# z direction
rZ_m = dphiZ_p[:,:,2:end]./fsign(dphiZ_p[:,:,1:end-1])
psiZ_m[:,:,1:Nz] = 0.5*FL(rZ_m).*(phi.value[2:Nr+1,2:Ntheta+1,1:Nz]-phi.value[2:Nr+1,2:Ntheta+1,2:Nz+1])
psiZ_m[:,:,Nz+1] = 0.0  # front boundary

re = rf[2:Nr+1]
rw = rf[1:Nr]

# find the velocity direction for the upwind scheme
ue_min = min(u.xvalue[2:Nr+1,:,:],0)
ue_max = max(u.xvalue[2:Nr+1,:,:],0)
uw_min = min(u.xvalue[1:Nr,:,:],0)
uw_max = max(u.xvalue[1:Nr,:,:],0)
vn_min = min(u.yvalue[:,2:Ntheta+1,:],0)
vn_max = max(u.yvalue[:,2:Ntheta+1,:],0)
vs_min = min(u.yvalue[:,1:Ntheta,:],0)
vs_max = max(u.yvalue[:,1:Ntheta,:],0)
wf_min = min(u.zvalue[:,:,2:Nz+1],0)
wf_max = max(u.zvalue[:,:,2:Nz+1],0)
wb_min = min(u.zvalue[:,:,1:Nz],0)
wb_max = max(u.zvalue[:,:,1:Nz],0)

# calculate the coefficients for the internal cells
AE = re.*ue_min./(DRp.*rp)
AW = -rw.*uw_max./(DRp.*rp)
AN = vn_min./(DTHETAp.*rp)
AS = -vs_max./(DTHETAp.*rp)
AF = wf_min./DZp
AB = -wb_max./DZp
APx = (re.*ue_max-rw.*uw_min)./(DRp.*rp)
APy = (vn_max-vs_min)./(DTHETAp.*rp)
APz = (wf_max-wb_min)./DZp

# Also correct for the boundary cells (not the ghost cells)
# Left boundary:
APx[1,:,:] = APx[1,:,:]-rw[1,:,:].*uw_max[1,:,:]./(2.0*DRp[1]*rp[1,:,:])
AW[1,:,:] = AW[1,:,:]/2.0
# Right boundary:
AE[end,:,:] = AE[end,:,:]/2.0
APx[end,:,:] = APx[end,:,:]+re[end,:,:].*ue_min[end,:,:]./(2.0*DRp[end]*rp[end,:,:])
# Bottom boundary:
APy[:,1,:] = APy[:,1,:]-vs_max[:,1,:]./(2.0*DTHETAp[1]*rp[:,1,:])
AS[:,1,:] = AS[:,1,:]/2.0
# Top boundary:
AN[:,end,:] = AN[:,end,:]/2.0
APy[:,end,:] = APy[:,end,:]+vn_min[:,end,:]./(2.0*DTHETAp[end]*rp[:,end,:])
# Back boundary:
APz[:,:,1] = APz[:,:,1]-wb_max[:,:,1]/(2.0*DZp[1])
AB[:,:,1] = AB[:,:,1]/2.0
# Front boundary:
AF[:,:,end] = AF[:,:,end]/2.0
APz[:,:,end] = APz[:,:,end]+wf_min[:,:,end]/(2.0*DZp[end])

AE = reshape(AE,mnx)
AW = reshape(AW,mnx)
AN = reshape(AN,mny)
AS = reshape(AS,mny)
AF = reshape(AF,mnz)
AB = reshape(AB,mnz)
APx = reshape(APx,mnx)
APy = reshape(APy,mny)
APz = reshape(APz,mnz)

# build the sparse matrix based on the numbering system
rowx_index = reshape(G[2:Nr+1,2:Ntheta+1,2:Nz+1],mnx)  # main diagonal x
iix[1:3*mnx] = repmat(rowx_index,3)
rowy_index = reshape(G[2:Nr+1,2:Ntheta+1,2:Nz+1],mny)  # main diagonal y
iiy[1:3*mny] = repmat(rowy_index,3)
rowz_index = reshape(G[2:Nr+1,2:Ntheta+1,2:Nz+1],mnz)  # main diagonal z
iiz[1:3*mnz] = repmat(rowz_index,3)
jjx[1:3*mnx] = [reshape(G[1:Nr,2:Ntheta+1,2:Nz+1],mnx); reshape(G[2:Nr+1,2:Ntheta+1,2:Nz+1],mnx); reshape(G[3:Nr+2,2:Ntheta+1,2:Nz+1],mnx)]
jjy[1:3*mny] = [reshape(G[2:Nr+1,1:Ntheta,2:Nz+1],mny); reshape(G[2:Nr+1,2:Ntheta+1,2:Nz+1],mny); reshape(G[2:Nr+1,3:Ntheta+2,2:Nz+1],mny)]
jjz[1:3*mnz] = [reshape(G[2:Nr+1,2:Ntheta+1,1:Nz],mnz); reshape(G[2:Nr+1,2:Ntheta+1,2:Nz+1],mnz); reshape(G[2:Nr+1,2:Ntheta+1,3:Nz+2],mnz)]
sx[1:3*mnx] = [AW; APx; AE]
sy[1:3*mny] = [AS; APy; AN]
sz[1:3*mnz] = [AB; APz; AF]

# calculate the TVD correction term
div_x = -(1./(DRp.*rp)).*(re.*(ue_max.*psiX_p[2:Nr+1,:,:]+ue_min.*psiX_m[2:Nr+1,:,:])-
              rw.*(uw_max.*psiX_p[1:Nr,:,:]+uw_min.*psiX_m[1:Nr,:,:]))
div_y = -(1./(DTHETAp.*rp)).*((vn_max.*psiY_p[:,2:Ntheta+1,:]+vn_min.*psiY_m[:,2:Ntheta+1,:])-
              (vs_max.*psiY_p[:,1:Ntheta,:]+vs_min.*psiY_m[:,1:Ntheta,:]))
div_z = -(1./DZp).*((wf_max.*psiZ_p[:,:,2:Nz+1]+wf_min.*psiZ_m[:,:,2:Nz+1])-
              (wb_max.*psiZ_p[:,:,1:Nz]+wb_min.*psiZ_m[:,:,1:Nz]))

# define the RHS Vector
RHS = zeros((Nr+2)*(Ntheta+2)*(Nz+2))
RHSx = zeros((Nr+2)*(Ntheta+2)*(Nz+2))
RHSy = zeros((Nr+2)*(Ntheta+2)*(Nz+2))
RHSz = zeros((Nr+2)*(Ntheta+2)*(Nz+2))

# assign the values of the RHS vector
row_index = rowx_index
RHS[row_index] = reshape(div_x+div_y+div_z,Nr*Ntheta*Nz)
RHSx[rowx_index] = reshape(div_x,Nr*Ntheta*Nz)
RHSy[rowy_index] = reshape(div_y,Nr*Ntheta*Nz)
RHSz[rowz_index] = reshape(div_z,Nr*Ntheta*Nz)

# build the sparse matrix
kx = 3*mnx
ky = 3*mny
kz = 3*mnz
Mx = sparse(iix[1:kx], jjx[1:kx], sx[1:kx], (Nr+2)*(Ntheta+2)*(Nz+2), (Nr+2)*(Ntheta+2)*(Nz+2))
My = sparse(iiy[1:ky], jjy[1:ky], sy[1:ky], (Nr+2)*(Ntheta+2)*(Nz+2), (Nr+2)*(Ntheta+2)*(Nz+2))
Mz = sparse(iiz[1:kz], jjz[1:kz], sz[1:kz], (Nr+2)*(Ntheta+2)*(Nz+2), (Nr+2)*(Ntheta+2)*(Nz+2))
M = Mx + My + Mz;

(M, RHS, Mx, My, Mz, RHSx, RHSy, RHSz)
end