linear_scalar_advection.jl

@doc raw"""
    LinearScalarAdvectionEquation2D

The linear scalar advection equation
```math
\partial_t u + a_1 \partial_1 u + a_2 \partial_2 u = 0
```
in two space dimensions with constant velocity `a`.
"""
struct LinearScalarAdvectionEquation2D{RealT<:Real} <: AbstractLinearScalarAdvectionEquation{2, 1}
  advectionvelocity::SVector{2, RealT}
end

function LinearScalarAdvectionEquation2D(a::NTuple{2,<:Real})
  LinearScalarAdvectionEquation2D(SVector(a))
end

function LinearScalarAdvectionEquation2D(a1::Real, a2::Real)
  LinearScalarAdvectionEquation2D(SVector(a1, a2))
end

# TODO Taal refactor, remove old constructors and replace them with default values
function LinearScalarAdvectionEquation2D()
  a = convert(SVector{2,Float64}, parameter("advectionvelocity"))
  LinearScalarAdvectionEquation2D(a)
end


get_name(::LinearScalarAdvectionEquation2D) = "LinearScalarAdvectionEquation2D"
varnames_cons(::LinearScalarAdvectionEquation2D) = SVector("scalar")
varnames_prim(::LinearScalarAdvectionEquation2D) = SVector("scalar")

# Calculates translated coordinates `x` for a periodic domain
function x_trans_periodic_2d(x, domain_length = SVector(2, 2), center = SVector(0, 0))
  x_normalized = x .- center
  x_shifted = x_normalized .% domain_length
  x_offset = ((x_shifted .< -0.5*domain_length) - (x_shifted .> 0.5*domain_length)) .* domain_length
  return center + x_shifted + x_offset
end

# Set initial conditions at physical location `x` for time `t`
function initial_condition_gauss(x, t, equation::LinearScalarAdvectionEquation2D)
  # Store translated coordinate for easy use of exact solution
  x_trans = x_trans_periodic_2d(x - equation.advectionvelocity * t)

  return @SVector [exp(-(x_trans[1]^2 + x_trans[2]^2))]
end

function initial_condition_convergence_test(x, t, equation::LinearScalarAdvectionEquation2D)
  # Store translated coordinate for easy use of exact solution
  x_trans = x - equation.advectionvelocity * t

  c = 1.0
  A = 0.5
  L = 2
  f = 1/L
  omega = 2 * pi * f
  scalar = c + A * sin(omega * sum(x_trans))
  return @SVector [scalar]
end

function initial_condition_sin_sin(x, t, equation::LinearScalarAdvectionEquation2D)
  # Store translated coordinate for easy use of exact solution
  x_trans = x - equation.advectionvelocity * t

  scalar = sin(2 * pi * x_trans[1]) * sin(2 * pi * x_trans[2])
  return @SVector [scalar]
end

function initial_condition_constant(x, t, equation::LinearScalarAdvectionEquation2D)
  return @SVector [2.0]
end

function initial_condition_linear_x_y(x, t, equation::LinearScalarAdvectionEquation2D)
  # Store translated coordinate for easy use of exact solution
  x_trans = x - equation.advectionvelocity * t

  return @SVector [sum(x_trans)]
end

function initial_condition_linear_x(x, t, equation::LinearScalarAdvectionEquation2D)
  # Store translated coordinate for easy use of exact solution
  x_trans = x - equation.advectionvelocity * t

  return @SVector [x_trans[1]]
end

function initial_condition_linear_y(x, t, equation::LinearScalarAdvectionEquation2D)
  # Store translated coordinate for easy use of exact solution
  x_trans = x - equation.advectionvelocity * t

  return @SVector [x_trans[2]]
end


# Apply boundary conditions
function boundary_condition_gauss(u_inner, orientation, direction, x, t, surface_flux_function,
                                   equation::LinearScalarAdvectionEquation2D)
  u_boundary = initial_condition_gauss(x, t, equation)

  # Calculate boundary flux
  if direction in (2, 4) # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, orientation, equation)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_boundary, u_inner, orientation, equation)
  end

  return flux
end


function boundary_condition_linear_x_y(u_inner, orientation, direction, x, t,
                                              surface_flux_function,
                                              equation::LinearScalarAdvectionEquation2D)
  u_boundary = initial_condition_linear_x_y(x, t, equation)

  # Calculate boundary flux
  if direction in (2, 4) # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, orientation, equation)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_boundary, u_inner, orientation, equation)
  end

  return flux
end


function boundary_condition_linear_x(u_inner, orientation, direction, x, t,
                                              surface_flux_function,
                                              equation::LinearScalarAdvectionEquation2D)
  u_boundary = initial_condition_linear_x(x, t, equation)

  # Calculate boundary flux
  if direction in (2, 4) # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, orientation, equation)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_boundary, u_inner, orientation, equation)
  end

  return flux
end


function boundary_condition_linear_y(u_inner, orientation, direction, x, t,
                                              surface_flux_function,
                                              equation::LinearScalarAdvectionEquation2D)
  u_boundary = initial_condition_linear_y(x, t, equation)

  # Calculate boundary flux
  if direction in (2, 4) # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, orientation, equation)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_boundary, u_inner, orientation, equation)
  end

  return flux
end


# Pre-defined source terms should be implemented as
# function source_terms_WHATEVER(ut, u, x, element_id, t, n_nodes, equation::LinearScalarAdvectionEquation2D)


# Calculate 1D flux in for a single point
@inline function calcflux(u, orientation, equation::LinearScalarAdvectionEquation2D)
  a = equation.advectionvelocity[orientation]
  return a * u
end


function flux_lax_friedrichs(u_ll, u_rr, orientation, equation::LinearScalarAdvectionEquation2D)
  a = equation.advectionvelocity[orientation]
  return 0.5 * ( a * (u_ll + u_rr) - abs(a) * (u_rr - u_ll) )
end


# Determine maximum stable time step based on polynomial degree and CFL number
function calc_max_dt(u, element_id, invjacobian, cfl,
                     equation::LinearScalarAdvectionEquation2D, dg)
  λ_max = maximum(abs, equation.advectionvelocity)
  return cfl * 2 / (nnodes(dg) * invjacobian * λ_max)
end

@inline have_constant_speed(::LinearScalarAdvectionEquation2D) = Val(true)

@inline function max_abs_speeds(eq::LinearScalarAdvectionEquation2D)
  return abs.(eq.advectionvelocity)
end


# Convert conservative variables to primitive
@inline cons2prim(u, equation::LinearScalarAdvectionEquation2D) = u

# Convert conservative variables to entropy variables
@inline cons2entropy(u, equation::LinearScalarAdvectionEquation2D) = u


# Calculate entropy for a conservative state `cons`
@inline entropy(u::Real, ::LinearScalarAdvectionEquation2D) = 0.5 * u^2
@inline entropy(u, equation::LinearScalarAdvectionEquation2D) = entropy(u[1], equation)


# Calculate total energy for a conservative state `cons`
@inline energy_total(u::Real, ::LinearScalarAdvectionEquation2D) = 0.5 * u^2
@inline energy_total(u, equation::LinearScalarAdvectionEquation2D) = energy_total(u[1], equation)