From e45418916aed2b9c5e24d20af4625c04b3cbb4a1 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Wed, 9 Dec 2020 21:27:59 -0500 Subject: [PATCH] #300 use Method of Manufactured Solutions --- Project.toml | 1 + src/MOL_discretization.jl | 36 +++++++------ test/MOL_1D_Linear_Diffusion.jl | 96 +++++++++++++++++++-------------- 3 files changed, 76 insertions(+), 57 deletions(-) diff --git a/Project.toml b/Project.toml index b765faf98..de5b02fbf 100644 --- a/Project.toml +++ b/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" diff --git a/src/MOL_discretization.jl b/src/MOL_discretization.jl index 180155111..25c508697 100644 --- a/src/MOL_discretization.jl +++ b/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,11 +244,13 @@ 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 ######################### @@ -259,13 +258,8 @@ function DiffEqBase.discretize(pdesys::PDESystem,discretization::MOLFiniteDiffer 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 diff --git a/test/MOL_1D_Linear_Diffusion.jl b/test/MOL_1D_Linear_Diffusion.jl index a100efca8..4ff6047be 100644 --- a/test/MOL_1D_Linear_Diffusion.jl +++ b/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