#300 use Method of Manufactured Solutions

Project.toml

@@ -13,6 +13,7 @@ LazyBandedMatrices = "d7e5e226-e90b-4449-9968-0f923699bf6f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 NNlib = "872c559c-99b0-510c-b3b7-b6c96a88d5cd"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"

src/MOL_discretization.jl

@@ -56,15 +56,13 @@ function get_bcs(bcs,tdomain,domain)
             ))
         end
         # Create value
-        if !(bcs[i].rhs isa ModelingToolkit.Symbolic)
-            γ = :(var=$(bcs[i].rhs))
-        else
-            γ = toexpr(bcs[i].rhs)
-        end
+        γ = toexpr(bcs[i].rhs)
         # Assign
         if key == "ic"
+            # Initial conditions always take the form u(0, x) ~ γ
             lhs_deriv_depvars_bcs[var][key] = γ
         else
+            # Boundary conditions can be more general
             lhs_deriv_depvars_bcs[var][key] = (α, β, γ)
         end
     end
@@ -232,8 +230,7 @@ function DiffEqBase.discretize(pdesys::PDESystem,discretization::MOLFiniteDiffer
     # TODO: extend to Neumann BCs and Robin BCs
     lhs_deriv_depvars_bcs = get_bcs(pdesys.bcs,tdomain,domain[2])
     u_ic = Array{Float64}(undef,len_of_indep_vars[2]-2,num_dep_vars)
-    u_left_bc = Array{Any}(undef,num_dep_vars)
-    u_right_bc = Array{Any}(undef,num_dep_vars)
+    robin_bc_func = Array{Any}(undef,num_dep_vars)
     Q = Array{RobinBC}(undef,num_dep_vars)
 
     for var in keys(dep_var_idx)
@@ -247,25 +244,22 @@ function DiffEqBase.discretize(pdesys::PDESystem,discretization::MOLFiniteDiffer
         # Boundary conditions depend on time and so will be evaluated at t within the
         # ODE function for the discretized PDE (below)
         # We use @RuntimeGeneratedFunction to do this efficiently
-        # For now we restrict to the case where only γ depends on time
         αl, βl, γl = bcs["left bc"]
         αr, βr, γr = bcs["right bc"]
-        u_left_bc[j] = (αl, βl, @RuntimeGeneratedFunction(:(t -> $(γl))))
-        u_right_bc[j] =  (αl, βl, @RuntimeGeneratedFunction(:(t -> $(γr))))
+        # Right now there is only one independent variable, so dx is always dx[2]
+        # i.e. the spatial variable
+        robin_bc_func[j] = @RuntimeGeneratedFunction(:(t -> begin
+            RobinBC(($(αl), $(βl), $(γl)), ($(αr), $(βr), $(γr)), $(dx[2]))
+        end))
     end
 
     ### Define the discretized PDE as an ODE function #########################
 
     function f(du,u,p,t)
 
         # Boundary conditions can vary with respect to time (but not space)
-        # For now assume that only γ depends on t
         for j in 1:num_dep_vars
-            αl, βl, γl = u_left_bc[j]
-            αr, βr, γr = u_right_bc[j]
-            # Right now there is only one independent variable, so dx is always dx[2]
-            # i.e. the spatial variable
-            Q[j] = RobinBC((αl, βl, γl(t)), (αr, βr, γr(t)), dx[2])
+            Q[j] = robin_bc_func[j](t)
         end
 
         for (var,disc) in dep_var_disc
@@ -281,5 +275,13 @@ function DiffEqBase.discretize(pdesys::PDESystem,discretization::MOLFiniteDiffer
     end
 
     # Return problem ##########################################################
-    return PDEProblem(ODEProblem(f,u_ic,(tdomain.lower,tdomain.upper),nothing),Q,X)
+    # The second entry, robin_bc_func, is stored as the "extrapolation" in the
+    # PDEProblem. Since this can in general be time-dependent, it must be evaluated
+    # at the correct time when used, e.g.
+    #   prob.extrapolation[1](t[i])*sol[:,1,i]
+    return PDEProblem(
+        ODEProblem(f,u_ic,(tdomain.lower,tdomain.upper),nothing),
+        robin_bc_func,
+        X
+    )
 end

test/MOL_1D_Linear_Diffusion.jl

@@ -7,6 +7,9 @@ using ModelingToolkit,DiffEqOperators,DiffEqBase,LinearAlgebra,Test,OrdinaryDiff
 
 # Tests
 @testset "Test 00: Dt(u(t,x)) ~ Dxx(u(t,x))" begin
+    # Method of Manufactured Solutions
+    u_exact = (x,t) -> exp.(-t) * cos.(x)
+
     # Parameters, variables, and derivatives
     @parameters t x
     @variables u(..)
@@ -15,13 +18,13 @@ using ModelingToolkit,DiffEqOperators,DiffEqBase,LinearAlgebra,Test,OrdinaryDiff
 
     # 1D PDE and boundary conditions
     eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
-    bcs = [u(0,x) ~ -x*(x-1)*sin(x),
-           u(t,0) ~ 0.0,
-           u(t,1) ~ 0.0]
+    bcs = [u(0,x) ~ cos(x),
+           u(t,0) ~ exp(-t),
+           u(t,Float64(pi)) ~ -exp(-t)]
 
     # Space and time domains
     domains = [t ∈ IntervalDomain(0.0,1.0),
-               x ∈ IntervalDomain(0.0,1.0)]
+               x ∈ IntervalDomain(0.0,Float64(pi))]
 
     # PDE system
     pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
@@ -36,20 +39,23 @@ using ModelingToolkit,DiffEqOperators,DiffEqBase,LinearAlgebra,Test,OrdinaryDiff
 
     # Solve ODE problem
     sol = solve(prob,Tsit5(),saveat=0.1)
+    x = prob.space[2]
+    t = sol.t
 
     # Plot and save results
     # using Plots
-    # plot(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,1]))
-    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,2]))
-    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,3]))
-    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,4]))
+    # plot()
+    # for i in 1:4
+    #     plot!(x,Array(prob.extrapolation[1](t[i])*sol.u[i]))
+    #     scatter!(x, u_exact(x, t[i]))
+    # end
     # savefig("MOL_1D_Linear_Diffusion_Test00.png")
 
-    # Test
-    n = size(sol,1)
-    t_f = size(sol,3)
-
-    @test sol[:,1,t_f] ≈ zeros(n) atol = 0.001;
+    # Test against exact solution
+    for i in 1:size(t,1)
+        u_approx = Array(prob.extrapolation[1](t[i])*sol.u[i])
+        @test u_approx ≈ u_exact(x, t[i]) atol=0.01
+    end
 end
 
 @testset "Test 01: Dt(u(t,x)) ~ D*Dxx(u(t,x))" begin
@@ -86,11 +92,14 @@ end
     sol = solve(prob,Tsit5(),saveat=0.1)
 
     # Plot and save results
+    # x = prob.space[2]
+    # t = sol.t
+
     # using Plots
-    # plot(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,1]))
-    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,2]))
-    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,3]))
-    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,4]))
+    # plot()
+    # for i in 1:4
+    #     plot!(x,Array(prob.extrapolation[1](t[i])*sol.u[i]))
+    # end
     # savefig("MOL_1D_Linear_Diffusion_Test01.png")
 
     # Test
@@ -138,11 +147,14 @@ end
     sol = solve(prob,Tsit5(),saveat=0.1)
 
     # Plot and save results
+    # x = prob.space[2]
+    # t = sol.t
+
     # using Plots
-    # plot(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,1]))
-    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,2]))
-    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,3]))
-    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,4]))
+    # plot()
+    # for i in 1:4
+    #     plot!(x,Array(prob.extrapolation[1](t[i])*sol.u[i]))
+    # end
     # savefig("MOL_1D_Linear_Diffusion_Test02.png")
 
     # Test
@@ -151,7 +163,10 @@ end
     @test sol[:,1,t_f] ≈ zeros(n) atol = 0.1;
 end
 
-@testset "Test 04: Dt(u(t,x)) ~ Dxx(u(t,x)), Dx(u(t,0)) ~ 0, Dx(u(t,1)) ~ 0" begin
+@testset "Test 03: Dt(u(t,x)) ~ Dxx(u(t,x)), Dx(u(t,0)) ~ 0, Dx(u(t,1)) ~ 0" begin
+    # Method of Manufactured Solutions
+    u_exact = (x,t) -> exp.(-t) * sin.(x)
+
     # Parameters, variables, and derivatives
     @parameters t x
     @variables u(..)
@@ -161,13 +176,13 @@ end
 
     # 1D PDE and boundary conditions
     eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
-    bcs = [u(0,x) ~ 0.5 + sin(2pi*x),
-           Dx(u(t,0)) ~ 0.0,
-           Dx(u(t,1)) ~ 0.0]
+    bcs = [u(0,x) ~ sin(x),
+           Dx(u(t,0)) ~ exp(-t),
+           Dx(u(t,Float64(pi))) ~ -exp(-t)]
 
     # Space and time domains
     domains = [t ∈ IntervalDomain(0.0,1.0),
-               x ∈ IntervalDomain(0.0,1.0)]
+               x ∈ IntervalDomain(0.0,Float64(pi))]
 
     # PDE system
     pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
@@ -182,22 +197,23 @@ end
 
     # Solve ODE problem
     sol = solve(prob,Tsit5(),saveat=0.1)
+    x = prob.space[2]
+    t = sol.t
 
     # Plot and save results
-    # using Plots
-    # plot(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,1]))
-    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,2]))
-    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,3]))
-    # plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,4]))
-    # savefig("MOL_1D_Linear_Diffusion_Test04.png")
-
-    # Test
-    # With zero flux at the boundaries, the solution should converge to the average of
-    # its initial condition, 0.5
-    n = size(sol,1)
-    t_f = size(sol,3)
-
-    @test sol[:,1,t_f] ≈ 0.5ones(n) atol = 0.001;
+    using Plots
+    plot()
+    for i in 1:4
+        plot!(x,Array(prob.extrapolation[1](t[i])*sol.u[i]))
+        scatter!(x, u_exact(x, t[i]))
+    end
+    savefig("MOL_1D_Linear_Diffusion_Test03.png")
+
+    # Test against exact solution
+    for i in 1:size(t,1)
+        u_approx = Array(prob.extrapolation[1](t[i])*sol.u[i])
+        @test u_approx ≈ u_exact(x, t[i]) atol=0.01
+    end
 end
 
 @testset "Test errors" begin