Add advection equation (#46)

* add advection equation

* format

* don't save gif

* format

* adjust tolerance for CI

* add test with vector advection velocity

docs/src/interpolation.md

@@ -68,7 +68,7 @@ depends only on one variable. A *radial basis function* kernel is a translation-
 
 Many radial symmetric kernels come with a parameter, the so-called *shape parameter* ``\varepsilon``, which can be used to control the "flatness"
 of the kernel. The shape parameter simply acts as a multiplicative factor to the norm, i.e. for a general radial-symmetric kernel we take
-``K(x, y) = \phi(\varepsilon\|x - y\|)``.
+``K(x, y) = \phi(\varepsilon\|x - y\|_2)``.
 
 The completion of the linear space of functions that is spanned by the basis given a specific kernel and a domain ``\Omega``,
 ``\mathcal{H}_{K, \Omega} = \overline{\text{span}\{K(\cdot, x), x\in\Omega\}}``, is called *native space* and is a (reproducing kernel) Hilbert space (RKHS),

examples/PDEs/advection_1d_basic.jl

@@ -0,0 +1,36 @@
+using KernelInterpolation
+using OrdinaryDiffEq
+using Plots
+
+# source term of advection equation
+f(t, x, equations) = 0.0
+pde = AdvectionEquation((0.5,), f)
+
+# initial condition
+u(t, x, equations) = exp(-100.0 * (x[1] - equations.advection_velocity[1] * t - 0.3)^2)
+
+n = 20
+nodeset_inner = homogeneous_hypercube(n, 0.01, 1.0)
+# only provide boundary condition at left boundary
+nodeset_boundary = NodeSet([0.0])
+g(t, x) = u(t, x, pde)
+
+kernel = WendlandKernel{1}(3, shape_parameter = 1.0)
+sd = Semidiscretization(pde, nodeset_inner, g, nodeset_boundary, u, kernel)
+tspan = (0.0, 1.0)
+ode = semidiscretize(sd, tspan)
+
+sol = solve(ode, Rosenbrock23(), saveat = 0.01)
+titp = TemporalInterpolation(sol)
+
+many_nodes = homogeneous_hypercube(200; dim = 1)
+
+plot()
+for t in 0.0:0.2:1.0
+    plot!(many_nodes, titp(t), training_nodes = false, linewidth = 2, color = :blue,
+          label = t == 0.0 ? "numerical" : "")
+    scatter!(many_nodes, u.(Ref(t), many_nodes, Ref(pde)), linewidth = 2,
+             linestyle = :dot, color = :red, markersize = 2,
+             label = t == 0.0 ? "analytical" : "")
+end
+plot!()

examples/PDEs/advection_3d_basic.jl

@@ -0,0 +1,40 @@
+using KernelInterpolation
+using OrdinaryDiffEq
+using LinearAlgebra: norm
+using WriteVTK: WriteVTK, paraview_collection
+
+# source term of advection equation
+f(t, x, equations) = 0.0
+pde = AdvectionEquation((0.5, 0.5, 0.5), f)
+
+# initial condition
+function u(t, x, equations)
+    exp(-20.0 * norm(x .- equations.advection_velocity .* t .- [0.3, 0.3, 0.3])^2)
+end
+
+n = 10
+nodeset_inner = homogeneous_hypercube(n, 0.01, 1.0; dim = 3)
+# don't provide any boundary condition
+nodeset_boundary = empty_nodeset(3, Float64)
+g(t, x) = 0.0
+
+kernel = WendlandKernel{3}(3, shape_parameter = 1.0)
+sd = Semidiscretization(pde, nodeset_inner, g, nodeset_boundary, u, kernel)
+tspan = (0.0, 1.0)
+ode = semidiscretize(sd, tspan)
+
+sol = solve(ode, Rosenbrock23(), saveat = 0.01)
+titp = TemporalInterpolation(sol)
+
+many_nodes = homogeneous_hypercube(15; dim = 3)
+OUT = "out"
+ispath(OUT) || mkpath(OUT)
+pvd = paraview_collection(joinpath(OUT, "solution"))
+for t in sol.t
+    KernelInterpolation.add_to_pvd(joinpath(OUT,
+                                            "advection_3d_basic_$(lpad(round(Int, t * 100), 4, '0'))"),
+                                   pvd, t, many_nodes,
+                                   titp(t).(many_nodes), u.(Ref(t), many_nodes, Ref(pde)),
+                                   keys = ["numerical", "analytical"])
+end
+WriteVTK.vtk_save(pvd)

src/KernelInterpolation.jl

@@ -2,7 +2,7 @@ module KernelInterpolation
 
 using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect
 using ForwardDiff: ForwardDiff
-using LinearAlgebra: Symmetric, norm, tr, muladd
+using LinearAlgebra: Symmetric, norm, tr, muladd, dot
 using Printf: @sprintf
 using ReadVTK: VTKFile, get_points, get_point_data, get_data
 using RecipesBase: RecipesBase, @recipe, @series
@@ -33,11 +33,11 @@ export GaussKernel, MultiquadricKernel, InverseMultiquadricKernel,
        Matern52Kernel, Matern72Kernel, RieszKernel,
        TransformationKernel, ProductKernel, SumKernel
 export phi, Phi, order
-export Laplacian
-export PoissonEquation, HeatEquation
+export Gradient, Laplacian
+export PoissonEquation, AdvectionEquation, HeatEquation
 export SpatialDiscretization, Semidiscretization, semidiscretize
-export NodeSet, separation_distance, dim, eachdim, values_along_dim, distance_matrix,
-       random_hypercube, random_hypercube_boundary, homogeneous_hypercube,
+export NodeSet, empty_nodeset, separation_distance, dim, eachdim, values_along_dim,
+       distance_matrix, random_hypercube, random_hypercube_boundary, homogeneous_hypercube,
        homogeneous_hypercube_boundary, random_hypersphere, random_hypersphere_boundary
 export interpolation_kernel, nodeset, coefficients, kernel_coefficients,
        polynomial_coefficients, polynomial_basis, polyvars, system_matrix,

src/differential_operators.jl

@@ -1,7 +1,7 @@
 abstract type AbstractDifferentialOperator end
 
 function (D::AbstractDifferentialOperator)(kernel::RadialSymmetricKernel, x, y)
-    @assert length(x) == length(y)
+    @assert length(x) == length(y) == dim(kernel)
     return save_call(D, kernel, x .- y)
 end
 
@@ -16,10 +16,27 @@ function save_call(D::AbstractDifferentialOperator, kernel::RadialSymmetricKerne
     return D(kernel, x)
 end
 
+"""
+    Gradient()
+
+The gradient operator. It can be called with a `RadialSymmetricKernel` and points
+`x` and `y` to evaluate the gradient of the `kernel` at `x - y`.
+"""
+struct Gradient <: AbstractDifferentialOperator
+end
+
+function Base.show(io::IO, ::Gradient)
+    print(io, "∇")
+end
+
+function (::Gradient)(kernel::RadialSymmetricKernel, x)
+    return ForwardDiff.gradient(x -> Phi(kernel, x), x)
+end
+
 """
     Laplacian()
 
-The Laplacian operator. It can be called with an `RadialSymmetricKernel` and points
+The Laplacian operator. It can be called with a `RadialSymmetricKernel` and points
 `x` and `y` to evaluate the Laplacian of the `kernel` at `x - y`.
 """
 struct Laplacian <: AbstractDifferentialOperator

src/discretization.jl

@@ -144,6 +144,8 @@ function semidiscretize(semi::Semidiscretization, tspan)
                                  Ref(semi.spatial_discretization.equations))
     c0 = semi.cache.kernel_matrix \ u0
     iip = true # is-inplace, i.e., we modify a vector when calling rhs!
+    # TODO: This defines an ODEProblem with a mass matrix, which is singular, i.e. the problem is a DAE.
+    # Many ODE solvers do not support DAEs.
     f = ODEFunction{iip}(rhs!, mass_matrix = semi.cache.mass_matrix)
     return ODEProblem(f, c0, tspan, semi)
 end

src/equations.jl

@@ -46,6 +46,36 @@ function rhs(t, nodeset::NodeSet, equations::AbstractTimeDependentEquation)
     end
 end
 
+@doc raw"""
+    AdvectionEquation(advection_velocity)
+
+Advection equation with advection velocity `advection_velocity`. The advection equation is defined as
+```math
+    \partial_t u + \mathbf{a}\cdot\nabla u = f,
+```
+where ``\mathbf{a}`` is the advection velocity and ``f`` a source term.
+"""
+struct AdvectionEquation{RealT, F} <: AbstractTimeDependentEquation where {RealT, F}
+    advection_velocity::Vector{RealT}
+    f::F
+
+    function AdvectionEquation(advection_velocity::Vector{RealT}, f) where {RealT}
+        return new{RealT, typeof(f)}(advection_velocity, f)
+    end
+
+    function AdvectionEquation(advection_velocity::NTuple, f)
+        return new{eltype(advection_velocity), typeof(f)}(collect(advection_velocity), f)
+    end
+end
+
+function Base.show(io::IO, ::AdvectionEquation)
+    print(io, "∂_t u + a⋅∇u = f")
+end
+
+function (equations::AdvectionEquation)(kernel::RadialSymmetricKernel, x, y)
+    return dot(equations.advection_velocity, Gradient()(kernel, x, y))
+end
+
 @doc raw"""
     HeatEquation(diffusivity, f)
 

src/nodes.jl

@@ -54,6 +54,15 @@ function NodeSet(nodes::AbstractVector{RealT}) where {RealT}
     return NodeSet([[node] for node in nodes])
 end
 
+"""
+    empty_nodeset(Dim, RealT)
+
+Create an empty [`NodeSet`](@ref).
+"""
+function empty_nodeset(Dim, RealT)
+    NodeSet{Dim, RealT}(Vector{MVector{Dim, RealT}}[], Inf)
+end
+
 function separation_distance(nodes::Vector{MVector{Dim, RealT}}) where {Dim, RealT}
     r = Inf
     for (i, x) in enumerate(nodes)

test/Project.toml

@@ -7,6 +7,7 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
 [compat]
 Aqua = "0.8.3"

test/test_examples_pde.jl

@@ -28,6 +28,19 @@ EXAMPLES_DIR = joinpath(examples_dir(), "PDEs")
                               l2=0.051463428227268404, linf=0.005234201701416197,
                               pde_test=true)
     end
+
+    @ki_testset "advection_1d_basic.jl" begin
+        @test_include_example(joinpath(EXAMPLES_DIR, "advection_1d_basic.jl"),
+                              l2=0.029243246170576723, linf=0.005575735203507293,
+                              pde_test=true, atol=2e-4) # stability issues
+    end
+
+    @ki_testset "advection_3d_basic.jl" begin
+        @test_include_example(joinpath(EXAMPLES_DIR, "advection_3d_basic.jl"),
+                              l2=0.05532511824096249, linf=0.004383701006257845,
+                              pde_test=true, tspan=(0.0, 0.1))
+    end
 end
 
+rm("out"; force = true, recursive = true)
 end # module

test/test_unit.jl

@@ -187,6 +187,15 @@ using Plots
         end
         @test length(point_data) == 0
         @test_nowarn rm("nodeset1.vtu", force = true)
+        nodeset1_1 = @test_nowarn empty_nodeset(2, Float64)
+        @test length(nodeset1_1) == 0
+        @test dim(nodeset1_1) == 2
+        @test separation_distance(nodeset1_1) == Inf
+        @test_nowarn push!(nodeset1_1, [0.0, 0.0])
+        @test length(nodeset1_1) == 1
+        @test separation_distance(nodeset1_1) == Inf
+        @test_nowarn push!(nodeset1_1, [1.0, 0.0])
+        @test separation_distance(nodeset1_1) == 0.5
 
         nodeset2 = @test_nowarn NodeSet([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]])
         @test dim(nodeset2) == 2
@@ -630,6 +639,9 @@ using Plots
         l = @test_nowarn Laplacian()
         @test_nowarn println(l)
         @test_nowarn display(l)
+        g = @test_nowarn Gradient()
+        @test_nowarn println(g)
+        @test_nowarn display(g)
         # Test if automatic differentiation gives same results as analytical derivatives
         # Define derivatives of Gauss kernel analytically and use them in the solver
         # instead of automatic differentiation
@@ -664,13 +676,16 @@ using Plots
         end
         kernel = GaussKernel{2}(shape_parameter = 0.5)
         x1 = [0.4, 0.6]
-        @test isapprox(AnalyticalLaplacian()(kernel, x1), Laplacian()(kernel, x1))
+        @test isapprox(l(kernel, x1), AnalyticalLaplacian()(kernel, x1))
+        @test isapprox(g(kernel, x1), [-0.17561908618411226, -0.2634286292761684])
         kernel = GaussKernel{3}(shape_parameter = 0.5)
         x2 = [0.1, 0.2, 0.3]
-        @test isapprox(AnalyticalLaplacian()(kernel, x2), Laplacian()(kernel, x2))
+        @test isapprox(l(kernel, x2), AnalyticalLaplacian()(kernel, x2))
+        @test isapprox(g(kernel, x2),
+                       [-0.04828027081287833, -0.09656054162575665, -0.14484081243863497])
         kernel = GaussKernel{4}(shape_parameter = 0.5)
         x3 = rand(4)
-        @test isapprox(AnalyticalLaplacian()(kernel, x3, x3), Laplacian()(kernel, x3, x3))
+        @test isapprox(l(kernel, x3, x3), AnalyticalLaplacian()(kernel, x3, x3))
     end
 
     @testset "PDEs" begin
@@ -692,6 +707,15 @@ using Plots
         # time-dependent PDEs
         # Passing a function
         f(t, x, equations) = x[1] + x[2] + t
+        advection = @test_nowarn AdvectionEquation((2.0, 0.5), f)
+        advection = @test_nowarn AdvectionEquation([2.0, 0.5], f)
+        @test_nowarn println(advection)
+        @test_nowarn display(advection)
+        @test KernelInterpolation.rhs(1.0, nodeset, advection) == [1.0, 2.0, 2.0, 3.0]
+        # Passing a vector
+        @test_nowarn advection = HeatEquation((2.0, 0.5), [1.0, 2.0, 2.0, 4.0])
+        @test KernelInterpolation.rhs(1.0, nodeset, advection) == [1.0, 2.0, 2.0, 4.0]
+
         heat = @test_nowarn HeatEquation(2.0, f)
         @test_nowarn println(heat)
         @test_nowarn display(heat)