Adding NonLinearDiffusion operator to derivate_operator.jl returning discretized function

src/DiffEqOperators.jl

@@ -59,7 +59,7 @@ end
 export MatrixFreeOperator
 export AnalyticalJacVecOperator, JacVecOperator, getops
 export AbstractDerivativeOperator, DerivativeOperator,
-       CenteredDifference, UpwindDifference
+       CenteredDifference, UpwindDifference, nonlinear_diffusion, nonlinear_diffusion!
 export DirichletBC, Dirichlet0BC, NeumannBC, Neumann0BC, RobinBC, GeneralBC, MultiDimBC, PeriodicBC,
        MultiDimDirectionalBC, ComposedMultiDimBC
 export compose, decompose, perpsize

src/derivative_operators/derivative_operator.jl

@@ -26,6 +26,36 @@ struct DerivativeOperator{T<:Real,N,Wind,T2,S1,S2<:SArray,T3,F} <: AbstractDeriv
     coeff_func              :: F
 end
 
+function nonlinear_diffusion!(du::AbstractVector{T}, second_differential_order::Int, first_differential_order::Int, approx_order::Int,
+    p::AbstractVector{T}, q::AbstractVector{T}, dx::Union{T , AbstractVector{T} , Real},
+    nknots::Int) where {T<:Real, N}
+    #q is given by bc*u , u being the unknown function
+    #p is given as some function of q , p being the diffusion coefficient
+
+    @assert approx_order>1 "approximation_order must be greater than 1."
+    if first_differential_order > 0 
+    du .= (CenteredDifference(first_differential_order,approx_order,dx,nknots)*q).*(CenteredDifference(second_differential_order,approx_order,dx,nknots)*p)
+    else 
+    du .= q[2:(nknots + 1)].*(CenteredDifference(second_differential_order,approx_order,dx,nknots)*p)
+    end
+
+    for l = 1:(second_differential_order - 1)
+    du .= du .+ binomial(second_differential_order,l)*(CenteredDifference(l + first_differential_order,approx_order,dx,nknots)*q).*(CenteredDifference(second_differential_order - l,approx_order,dx,nknots)*p)
+    end
+
+    du .= du .+ (CenteredDifference(first_differential_order + second_differential_order,approx_order,dx,nknots)*q).*p[2:(nknots + 1)]
+
+end
+
+# An out of place workaround for the mutating version
+function nonlinear_diffusion(second_differential_order::Int, first_differential_order::Int, approx_order::Int,
+    p::AbstractVector{T}, q::AbstractVector{T}, dx::Union{T , AbstractVector{T} , Real},
+    nknots::Int) where {T<:Real, N}
+
+    du = similar(q,length(q) - 2)
+    return nonlinear_diffusion!(du,second_differential_order,first_differential_order,approx_order,p,q,dx,nknots)
+end
+
 struct CenteredDifference{N} end
 
 function CenteredDifference{N}(derivative_order::Int,

test/Fast_Diffusion.jl

@@ -0,0 +1,50 @@
+# # Fast-Diffusion Problem
+#
+# This example demonstrates the use of 'NonLinearDiffusion' operator to solve time-dependent non-linear diffusion PDEs with coefficient having dependence on the unknown function. 
+# Here we consider a fast diffusion problem with Dirichlet BCs on unit interval:
+# ∂ₜu = ∂ₓ(k*∂ₓu) 
+# k = 1/u²   
+# u(x=0,t) = exp(-t)
+# u(x=1,t) = 1/(1.0 + exp(2t))
+# u(x, t=0) = u₀(x)
+using Test
+using DiffEqOperators, OrdinaryDiffEq
+
+@testset "1D Nonlinear fast diffusion equation" begin
+    # The analytical solution for this is given by :
+
+    u_analytic(x, t) = 1 / sqrt(x^2 + exp(2*t))
+
+    #
+    # Reproducing it numerically 
+    #
+
+
+    nknots = 100
+    h = 1.0/(nknots+1)
+    knots = range(h, step=h, length=nknots)
+    n = 1                                   # Outer differential order
+    m = 1                                   # Inner differential order
+    approx_ord = 2                               
+
+    u0 = u_analytic.(knots,0.0)
+    du = similar(u0)
+    t0 = 0.0
+    t1 = 1.0
+
+    function f(du,u,p,t)
+        bc = DirichletBC(exp(-t),(1.0 + exp(2*t))^(-0.5))
+        l = bc*u
+        k = l.^(-2)                        # Diffusion Coefficient
+        nonlinear_diffusion!(du,n,m,approx_ord,k,l,h,nknots)
+    end
+
+    prob = ODEProblem(f, u0, (t0, t1))
+    alg = KenCarp4()
+    sol = solve(prob,alg)
+
+    for t in t0:0.1:t1
+        @test u_analytic.(knots, t) ≈ sol(t) rtol=1e-3
+    end
+
+end

test/runtests.jl

@@ -22,6 +22,7 @@ if GROUP == "All" || GROUP == "Interface"
     @time @safetestset "Validate Boundary Padded Array Concretization" begin include("boundary_padded_array.jl") end
     @time @safetestset "Validate Higher Dimensional Boundary Extension" begin include("multi_dim_bc_test.jl") end
     @time @safetestset "2nd order check" begin include("2nd_order_check.jl") end
+    @time @safetestset "Non-linear Diffusion" begin include("Fast_Diffusion.jl") end
     #@time @safetestset "KdV" begin include("KdV.jl") end # KdV times out and all fails
     #@time @safetestset "Heat Equation" begin include("heat_eqn.jl") end
     @time @safetestset "Matrix-Free Operators" begin include("matrixfree.jl") end