#378 use enum

src/DiffEqOperators.jl

@@ -67,6 +67,6 @@ export compose, decompose, perpsize
 export discretize, symbolic_discretize
 
 export GhostDerivativeOperator
-export MOLFiniteDifference
+export MOLFiniteDifference, center_align, edge_align
 export BoundaryConditionError
 end # module

src/MOLFiniteDifference/MOL_discretization.jl

@@ -1,15 +1,18 @@
 using ModelingToolkit: operation, istree, arguments
 # Method of lines discretization scheme
+
+@enum GridAlign center_align edge_align
+center_align
 struct MOLFiniteDifference{T,T2} <: DiffEqBase.AbstractDiscretization
     dxs::T
     time::T2
     upwind_order::Int
     centered_order::Int
-    grid_align::String
+    grid_align::GridAlign
 end
 
 # Constructors. If no order is specified, both upwind and centered differences will be 2nd order
-MOLFiniteDifference(dxs, time; upwind_order = 1, centered_order = 2, grid_align="center") =
+MOLFiniteDifference(dxs, time; upwind_order = 1, centered_order = 2, grid_align=center_align) =
     MOLFiniteDifference(dxs, time, upwind_order, centered_order, grid_align)
 
 function SciMLBase.symbolic_discretize(pdesys::ModelingToolkit.PDESystem,discretization::DiffEqOperators.MOLFiniteDifference)
@@ -34,11 +37,11 @@ function SciMLBase.symbolic_discretize(pdesys::ModelingToolkit.PDESystem,discret
     tspan = (tdomain.domain.lower,tdomain.domain.upper)
 
     # Define the grid on which the dependent variables will be evaluated (see #378)
-    # "center" is recommended for Dirichlet BCs
-    # "edge" is recommended for Neumann BCs (spatial discretization is conservative)
-    if grid_align == "center"
+    # center_align is recommended for Dirichlet BCs
+    # edge_align is recommended for Neumann BCs (spatial discretization is conservative)
+    if grid_align == center_align
         grid = space
-    elseif grid_align == "edge"
+    elseif grid_align == edge_align
         # boundary conditions implementation assumes centered_order=2
         @assert discretization.centered_order==2
         # construct grid including ghost nodes beyond outer edges
@@ -72,7 +75,7 @@ function SciMLBase.symbolic_discretize(pdesys::ModelingToolkit.PDESystem,discret
 
     # Generate map from variable (e.g. u(t,0)) to discretized variable (e.g. u₁(t))
     subvar(depvar) = substitute.((depvar,),edgevals)
-    if grid_align == "center"
+    if grid_align == center_align
         # depvarbcmaps will dictate what to replace the variable terms with in the bcs
         # replace u(t,0) with u₁, etc
         depvarbcmaps = reduce(vcat,[subvar(depvar) .=> edgevar for (depvar, edgevar) in zip(depvars, edgevars)])
@@ -94,7 +97,7 @@ function SciMLBase.symbolic_discretize(pdesys::ModelingToolkit.PDESystem,discret
         depvarderivbcmaps = reduce(vcat,[subderivar(depvar, s) .=> derivars[i]
                                          for (i, depvar) in enumerate(depvars) for s in nottime])
 
-        if grid_align == "edge"
+        if grid_align == edge_align
             # Constructs symbolic spatially discretized terms of the form e.g. (u₁ + u₂) / 2 
             bcvars = [[dot(ones(2)/2,depvar[left_idxs]), dot(ones(2)/2,depvar[right_idxs(1)])]
                       for depvar in depvarsdisc]

test/MOL/MOL_1D_Diffusion.jl

@@ -32,7 +32,7 @@ using ModelingToolkit: Differential
     dx = range(0.0,Float64(π),length=30)
     order = 2
     discretization = MOLFiniteDifference([x=>dx],t)
-    discretization_edge = MOLFiniteDifference([x=>dx],t;grid_align="edge")
+    discretization_edge = MOLFiniteDifference([x=>dx],t;grid_align=edge_align)
     # Explicitly specify order of centered difference
     discretization_centered = MOLFiniteDifference([x=>dx],t;centered_order=order)
     # Higher order centered difference
@@ -45,7 +45,7 @@ using ModelingToolkit: Differential
         # Solve ODE problem
         sol = solve(prob,Tsit5(),saveat=0.1)
 
-        if disc.grid_align == "center"
+        if disc.grid_align == center_align
             x = dx[2:end-1]
         else
             x = (dx[1:end-1]+dx[2:end])/2
@@ -171,7 +171,7 @@ end
     dx = range(0.0,Float64(π),length=300)
     order = 2
     discretization = MOLFiniteDifference([x=>dx],t)
-    discretization_edge = MOLFiniteDifference([x=>dx],t;grid_align="center")
+    discretization_edge = MOLFiniteDifference([x=>dx],t;grid_align=center_align)
 
     # Convert the PDE problem into an ODE problem
     for disc in [discretization, discretization_edge]
@@ -180,7 +180,7 @@ end
         # Solve ODE problem
         sol = solve(prob,Tsit5(),saveat=0.1)
 
-        if disc.grid_align == "center"
+        if disc.grid_align == center_align
             x_sol = dx[2:end-1]
         else
             x_sol = (dx[1:end-1]+dx[2:end])/2
@@ -225,7 +225,7 @@ end
     dx = range(0.0,Float64(π),length=30)
     order = 2
     discretization = MOLFiniteDifference([x=>dx],t)
-    discretization_edge = MOLFiniteDifference([x=>dx],t;grid_align="edge")
+    discretization_edge = MOLFiniteDifference([x=>dx],t;grid_align=edge_align)
 
     # Convert the PDE problem into an ODE problem
     for disc ∈ [discretization, discretization_edge]
@@ -234,7 +234,7 @@ end
         # Solve ODE problem
         sol = solve(prob,Tsit5(),saveat=0.1)
 
-        if disc.grid_align == "center"
+        if disc.grid_align == center_align
             x = dx[2:end-1]
         else
             x = (dx[1:end-1]+dx[2:end])/2
@@ -247,7 +247,6 @@ end
         for i in 1:length(sol)
             exact = u_exact(x, t[i])
             u_approx = sol.u[i]
-            @show maximum(abs.(u_approx - exact))
             @test u_approx ≈ exact atol=0.01
             # test mass conservation
             integral_u_approx = sum(u_approx * dx[2])
@@ -329,7 +328,7 @@ end
     dx = 0.01
     order = 2
     discretization = MOLFiniteDifference([x=>dx],t)
-    discretization_edge = MOLFiniteDifference([x=>dx],t;grid_align="edge")
+    discretization_edge = MOLFiniteDifference([x=>dx],t;grid_align=edge_align)
 
     for disc ∈ [discretization, discretization_edge]
         # Convert the PDE problem into an ODE problem
@@ -338,7 +337,7 @@ end
         # Solve ODE problem
         sol = solve(prob,Tsit5(),saveat=0.1)
         x = (-1:dx:1)
-        if disc.grid_align == "center"
+        if disc.grid_align == center_align
             x = x[2:end-1]
         else
             x = (x[1:end-1].+x[2:end])/2

test/MOL/MOLtest.jl

@@ -16,7 +16,7 @@ domains = [t ∈ IntervalDomain(0.0,1.0),
            x ∈ IntervalDomain(0.0,1.0)]
 
 pdesys = PDESystem(eqs,bcs,domains,[t,x],[u(t,x),v(t,x)])
-discretization = MOLFiniteDifference([x=>0.1],t;grid_align="edge")
+discretization = MOLFiniteDifference([x=>0.1],t;grid_align=edge_align)
 prob = discretize(pdesys,discretization) # This gives an ODEProblem since it's time-dependent
 sol = solve(prob,Tsit5())