.

Author: xtalax

src/bcs/generate_bc_eqs.jl

@@ -4,9 +4,9 @@ function _boundary_rules(s, orders, x, val)
 
     args = substitute.(args, (x=>val,))
 
-    rules = [operation(u)(args...) => (operation(u)(args...), x) for u in s.ū]
+    rules = [operation(u)(args...) => [operation(u)(args...), x] for u in s.ū]
 
-    return vcat(rules, vec([(Differential(x)^d)(operation(u)(args...)) => (operation(u)(args...), x) for d in orders[x], u in s.ū]))
+    return vcat(rules, vec([(Differential(x)^d)(operation(u)(args...)) => [operation(u)(args...), x] for d in orders[x], u in s.ū]))
 end
 
 function generate_boundary_matching_rules(s, orders)

src/discretization/generate_finite_difference_rules.jl

@@ -142,7 +142,7 @@ end
 end
 
 @inline function generate_nonlinlap_rules(II, s, derivweights, terms)
-    rules = [@rule *(~~c, $(Differential(x))(*(~~a, $(Differential(x))(u), ~~b)), ~~d) => cartesian_nonlinear_laplacian(*(a..., b...), II, derivweights, s, x, u) for x in s.x̄, u in s.ū]
+    rules = [@rule *(~~c, $(Differential(x))(*(~~a, $(Differential(x))(u), ~~b)), ~~d) => *(~~c,cartesian_nonlinear_laplacian(*(a..., b...), II, derivweights, s, x, u), ~~d) for x in s.x̄, u in s.ū]
 
     rules = [@rule $(Differential(x))(*(~~a, $(Differential(x))(u), ~~b)) => cartesian_nonlinear_laplacian(*(a..., b...), II, derivweights, s, x, u) for x in s.x̄, u in s.ū]
 
@@ -162,6 +162,10 @@ end
 @inline function generate_spherical_diffusion_rules(II, s, derivweights, terms)
     rules = vec([@rule *(~~a, 1/(r^2), ($(Differential(r))(*(~~c, (r^2), ~~d, $(Differential(r))(u), ~~e))), ~~b) => *(~a..., spherical_diffusion(*(~c..., ~d..., ~e...), II, derivweights, s, r, u), ~b...)
             for r in s.x̄, u in s.ū])
+
+    rules = vcat(rules, vec([@rule *(~~a, (r^2)^-2, ($(Differential(r))(*(~~c, (r^2), ~~d, $(Differential(r))(u), ~~e))), ~~b) => *(~a..., spherical_diffusion(*(~c..., ~d..., ~e...), II, derivweights, s, r, u), ~b...)
+            for r in s.x̄, u in s.ū]))
+
     rules = vcat(rules, vec([@rule /(($(Differential(r))(*(~~c, (r^2), ~~d, $(Differential(r))(u), ~~e))), (r^2)) => spherical_diffusion(*(~c..., ~d..., ~e...), II, derivweights, s, r, u)
     for r in s.x̄, u in s.ū]))
 

test/MOL_original.jl

@@ -1,11 +1,12 @@
 using ModelingToolkit: operation, istree, arguments, variable
 import DomainSets
+using DiffEqOperators, DiffEqBase, SciMLBase
 
 # Method of lines discretization scheme
 
 @enum GridAlign center_align edge_align
 
-struct MOLFiniteDifference_origial{T,T2} <: DiffEqBase.AbstractDiscretization
+struct MOLFiniteDifference_original{T,T2} <: DiffEqBase.AbstractDiscretization
     dxs::T
     time::T2
     upwind_order::Int
@@ -41,7 +42,7 @@ function calculate_weights_spherical_original(order::Int, x0::T, x::AbstractVect
 end
 
 
-function symbolic_discretize_original(pdesys::ModelingToolkit.PDESystem,discretization::DiffEqOperators.MOLFiniteDifference)
+function symbolic_discretize_original(pdesys::ModelingToolkit.PDESystem,discretization::MOLFiniteDifference_original)
     grid_align = discretization.grid_align
     pdeeqs = pdesys.eqs isa Vector ? pdesys.eqs : [pdesys.eqs]
     t = discretization.time
@@ -398,11 +399,13 @@ function symbolic_discretize_original(pdesys::ModelingToolkit.PDESystem,discreti
         return sys, nothing
     else
         sys = ODESystem(vcat(alleqs,unique(bceqs)),t,vec(reduce(vcat,vec(alldepvarsdisc))),ps,defaults=Dict(defaults),name=pdesys.name)
+        @show alleqs
+        @show bceqs
         return sys, tspan
     end
 end
 
-function discretize_original(pdesys::ModelingToolkit.PDESystem,discretization::DiffEqOperators.MOLFiniteDifference)
+function discretize_original(pdesys::ModelingToolkit.PDESystem,discretization::MOLFiniteDifference_original)
     sys, tspan = SciMLBase.symbolic_discretize(pdesys,discretization)
     if tspan == nothing
         return prob = NonlinearProblem(sys, ones(length(sys.states)))

test/pde_systems/MOL_1D_Linear_Diffusion.jl

@@ -231,7 +231,7 @@ end
     discretization_edge = MOLFiniteDifference([x=>dx_],t;grid_align=edge_align, approx_order=6)
 
     # Convert the PDE problem into an ODE problem
-    for disc ∈ [discretization, discretization_edge]
+    for (j,disc) ∈ enumerate([discretization, discretization_edge])
         prob = discretize(pdesys,disc)
 
         # Solve ODE problem
@@ -243,7 +243,15 @@ end
             x = ((0.0-dx_/2): dx_ : (Float64(π)+dx_/2))[2:end-1]
         end
         t = sol.t
-
+        # Plots
+        using Plots
+        anim = @animate for (i,T) in enumerate(t) 
+            exact = u_exact(x, T)
+            plot(x, exact, seriestype = :scatter,label="Analytic solution")
+            plot!(x, sol.u[i], label="Numeric solution")
+            plot!(x, log.(abs.(exact-sol.u[i])), label="Error at t = $(t[i])")
+        end
+        gif(anim, "plots/MOL_Linear_Diffusion_1D_Test03a_$disc.gif", fps = 5)
         # Test against exact solution
         # exact integral based on Neumann BCs
         integral_u_exact = t -> sum(sol.u[1] * dx[2]) + 2 * (exp(-t) - 1)
@@ -328,7 +336,8 @@ end
     # PDE system
     @named pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t,x)])
 
-    # Method of lines discretization
+    # Method of lines discretizationinclusions
+
     dx = 0.01
     order = 4
     discretization = MOLFiniteDifference([x=>dx],t; approx_order=order)