Add an analytic solution field for PDESystem (#2082)

Author: xtalax

src/systems/pde/pdesystem.jl

@@ -61,6 +61,19 @@ struct PDESystem <: ModelingToolkit.AbstractMultivariateSystem
     """
     systems::Vector
     """
+    analytic: A vector of explicit symbolic expressions for the analytic solutions of each
+    dependent variable. e.g. `analytic = [u(t, x) ~ a*sin(c*t) * cos(k*x)]`
+    """
+    analytic::Any
+    """
+    analytic_func: A vector of functions for the analytic solutions of each dependent
+    variable. Will be generated from `analytic` if not provided. Should have the same
+    argument signature as the variable, and a `ps` argument as the last argument,
+    which takes an indexable of parameter values in the order you specified them in `ps`.
+    e.g. `analytic_func = [u(t, x) => (ps, t, x) -> ps[1]*sin(ps[2]*t) * cos(ps[3]*x)]`
+    """
+    analytic_func::Any
+    """
     name: the name of the system
     """
     name::Symbol
@@ -78,15 +91,40 @@ struct PDESystem <: ModelingToolkit.AbstractMultivariateSystem
                                    systems = [],
                                    connector_type = nothing,
                                    metadata = nothing,
+                                   analytic = nothing,
+                                   analytic_func = nothing,
                                    gui_metadata = nothing,
                                    checks::Union{Bool, Int} = true,
                                    name)
         if checks == true || (checks & CheckUnits) > 0
             all_dimensionless([dvs; ivs; ps]) || check_units(eqs)
         end
+
         eqs = eqs isa Vector ? eqs : [eqs]
-        new(eqs, bcs, domain, ivs, dvs, ps, defaults, connector_type, systems, name,
-            metadata, gui_metadata)
+
+        if !isnothing(analytic)
+            analytic = analytic isa Vector ? analytic : [analytic]
+            if length(analytic) != length(dvs)
+                throw(ArgumentError("The number of analytic solutions must match the number of dependent variables"))
+            end
+
+            if isnothing(analytic_func)
+                analytic_func = map(analytic) do eq
+                    args = arguments(eq.lhs)
+                    p = ps isa SciMLBase.NullParameters ? [] : map(a -> a.first, ps)
+                    args = vcat(DestructuredArgs(p), args)
+                    ex = Func(args, [], eq.rhs) |> toexpr
+                    eq.lhs => @RuntimeGeneratedFunction(ex)
+                end
+            end
+        end
+
+        if !isnothing(analytic_func)
+            analytic_func = analytic_func isa Dict ? analytic_func : analytic_func |> Dict
+        end
+
+        new(eqs, bcs, domain, ivs, dvs, ps, defaults, connector_type, systems, analytic,
+            analytic_func, name, metadata, gui_metadata)
     end
 end
 

test/pde.jl

@@ -1,4 +1,4 @@
-using ModelingToolkit, DiffEqBase, LinearAlgebra
+using ModelingToolkit, DiffEqBase, LinearAlgebra, Test
 
 # Define some variables
 @parameters t x
@@ -13,7 +13,17 @@ bcs = [u(0, x) ~ -h * x * (x - 1) * sin(x),
 domains = [t ∈ (0.0, 1.0),
     x ∈ (0.0, 1.0)]
 
-@named pdesys = PDESystem(eq, bcs, domains, [t, x], [u])
+analytic = [u(t, x) ~ -h * x * (x - 1) * sin(x) * exp(-2 * h * t)]
+analytic_function = (ps, t, x) -> -ps[1] * x * (x - 1) * sin(x) * exp(-2 * ps[1] * t)
+
+@named pdesys = PDESystem(eq, bcs, domains, [t, x], [u], [h => 1], analytic = analytic)
 @show pdesys
 
 @test all(isequal.(independent_variables(pdesys), [t, x]))
+
+dx = 0:0.1:1
+dt = 0:0.1:1
+
+# Test generated analytic_func
+@test all(pdesys.analytic_func[u(t, x)]([2], disct, discx) ≈
+          analytic_function([2], disct, discx) for disct in dt, discx in dx)

test/runtests.jl

@@ -20,6 +20,7 @@ using SafeTestsets, Test
 @safetestset "SteadyStateSystem Test" begin include("steadystatesystems.jl") end
 @safetestset "SDESystem Test" begin include("sdesystem.jl") end
 @safetestset "NonlinearSystem Test" begin include("nonlinearsystem.jl") end
+@safetestset "PDE Construction Test" begin include("pde.jl") end
 @safetestset "JumpSystem Test" begin include("jumpsystem.jl") end
 @safetestset "Constraints Test" begin include("constraints.jl") end
 @safetestset "Reduction Test" begin include("reduction.jl") end
@@ -31,7 +32,6 @@ using SafeTestsets, Test
 @safetestset "State Selection Test" begin include("state_selection.jl") end
 @safetestset "Symbolic Event Test" begin include("symbolic_events.jl") end
 @safetestset "Stream Connnect Test" begin include("stream_connectors.jl") end
-@safetestset "PDE Construction Test" begin include("pde.jl") end
 @safetestset "Lowering Integration Test" begin include("lowering_solving.jl") end
 @safetestset "Test Big System Usage" begin include("bigsystem.jl") end
 @safetestset "Depdendency Graph Test" begin include("dep_graphs.jl") end