Skip to content

Commit 66ddbad

Browse files
xtalaxxtalax
committed
Fix indexing
1 parent 3739604 commit 66ddbad

File tree

7 files changed

+40
-34
lines changed

7 files changed

+40
-34
lines changed

‎src/MethodOfLines.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,7 @@ module MethodOfLines
1313
include("discretization/fornberg.jl")
1414
include("discretization/discretize_vars.jl")
1515
include("discretization/differential_discretizer.jl")
16-
include("discretization/generate_finite_difference_rules.jl")
16+
include("discretization/generate_rules.jl")
1717
include("bcs/generate_bc_eqs.jl")
1818

1919
include("discretization/MOL_discretization.jl")

‎src/bcs/generate_bc_eqs.jl‎

Lines changed: 6 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -2,9 +2,8 @@
22
"""
33
Mutates bceqs and u0 by finding relevant equations and discretizing them
44
"""
5-
function generate_u0_and_bceqs!!(u0, bceqs, bcs, s, depvar_ops, tspan, discretization)
6-
grid_align = discretization.grid_align
7-
t=discretization.time
5+
function generate_u0_and_bceqs!!(u0, bceqs, bcs, s, depvar_ops, tspan, derivweights)
6+
t=s.time
87

98
if t === nothing
109
initmaps = s.vars
@@ -21,13 +20,15 @@ function generate_u0_and_bceqs!!(u0, bceqs, bcs, s, depvar_ops, tspan, discretiz
2120
# Assume in the form `u(...) ~ ...` for now
2221
initindex = findfirst(isequal(bc.lhs), initmaps)
2322
if initindex !== nothing
24-
push!(u0,vec(s.discvars[initindex] .=> substitute.((bc.rhs,),gridvals(s))))
23+
push!(u0,vec(s.discvars[s.vars[initindex]] .=> substitute.((bc.rhs,),gridvals(s))))
2524
end
2625
else
2726
# boundary condition
2827
# TODO: Seperate out Iedge in to individual boundaries and map each seperately to save time and catch the periodic case, have a look at the old maps for how to match the bcs
2928
push!(bceqs, vec(map(s.Iedge) do II
30-
rules = generate_finite_difference_rules(II, s, bc, discretization)
29+
rules = generate_finite_difference_rules(II, s, bc, derivweights)
30+
rules = vcat(rules, generate_bc_rules(II, s, bc , derivweights))
31+
3132
substitute(bc.lhs, rules) ~ substitute(bc.rhs, rules)
3233
end))
3334
end

‎src/discretization/MOL_discretization.jl‎

Lines changed: 8 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -70,20 +70,21 @@ function SciMLBase.symbolic_discretize(pdesys::PDESystem, discretization::Method
7070
# Get the grid
7171
s = DiscreteSpace(domain, depvars, indvars, nottime, discretization)
7272

73+
# Get the finite difference weights
7374
derivweights = DifferentialDiscretizer(pde, s, discretization)
7475

75-
interior = s.Igrid[[2:(length(s, x)-1) for x in s.nottime]]
76+
# Get the boundary conditions
77+
generate_u0_and_bceqs!!(u0, bceqs, pdesys.bcs, s, depvar_ops, tspan, derivweights)
7678

77-
### PDE EQUATIONS ###
78-
# Create a stencil in the required dimension centered around 0
79-
# e.g. (-1,0,1) for 2nd order, (-2,-1,0,1,2) for 4th order, etc
80-
# For all cartesian indices in the interior, generate finite difference rules
79+
# Find the indexes on the interior
80+
interior = Iinterior(s)
81+
82+
# Discretize the equation on the interior
8183
pdeeqs = vec(map(interior) do II
8284
rules = generate_finite_difference_rules(II, s, pde, derivweights)
8385
substitute(pde.lhs,rules) ~ substitute(pde.rhs,rules)
8486
end)
8587

86-
generate_u0_and_bceqs!!(u0, bceqs, pdesys.bcs, s, depvar_ops, tspan, discretization)
8788
push!(alleqs,pdeeqs)
8889
push!(alldepvarsdisc, reduce(vcat, s.discvars))
8990
end
@@ -107,6 +108,7 @@ function SciMLBase.symbolic_discretize(pdesys::PDESystem, discretization::Method
107108
return sys, nothing
108109
else
109110
# * In the end we have reduced the problem to a system of equations in terms of Dt that can be solved by the `solve` method.
111+
println(alleqs, bceqs)
110112
sys = ODESystem(vcat(alleqs, unique(bceqs)), t, vec(reduce(vcat, vec(alldepvarsdisc))), ps, defaults=Dict(defaults), name=pdesys.name)
111113
return sys, tspan
112114
end

‎src/discretization/differential_discretizer.jl‎

Lines changed: 10 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -28,20 +28,22 @@ function DifferentialDiscretizer(pde, s, discretization)
2828
return DifferentialDiscretizer{eltype(orders), typeof(Dict(differentialmap))}(approx_order, Dict(differentialmap), Dict(nonlinlap), Dict(orders))
2929
end
3030

31+
3132
# ufunc is a function that returns the correct discretization indexed at Itap, it is designed this way to allow for central differences of arbitrary expressions which may be needed in some schemes
3233
function central_difference(D, II, s, jx, u, ufunc)
3334
j, x = jx
3435
# unit index in direction of the derivative
3536
I1 = unitindices(nparams(s))[j]
3637
# offset is important due to boundary proximity
38+
3739
if II[j] <= D.boundary_point_count
3840
weights = D.low_boundary_coefs[II[j]]
39-
offset = D.boundary_point_count - II[j] + 1
40-
Itap = [II + (i+offset)*I1 for i in half_range(D.boundary_stencil_length)]
41+
offset = 1 - II[j]
42+
Itap = [II + (i+offset)*I1 for i in 0:(D.boundary_stencil_length-1)]
4143
elseif II[j] > (length(s, x) - D.boundary_point_count)
4244
weights = D.high_boundary_coefs[length(s, x)-II[j]+1]
43-
offset = length(s, x) - II[j] - D.boundary_point_count
44-
Itap = [II + (i+offset)*I1 for i in half_range(D.boundary_stencil_length)]
45+
offset = length(s, x) - II[j]
46+
Itap = [II + (i+offset)*I1 for i in (-D.boundary_stencil_length+1):1:0]
4547
else
4648
weights = D.stencil_coefs
4749
Itap = [II + i*I1 for i in half_range(D.stencil_length)]
@@ -60,13 +62,12 @@ function _get_weights_and_stencil(D, II, I1, s, k, j, x)
6062
# The low boundary coeffs has a heirarchy of coefficients following: number of indices from boundary -> which half offset point does it correspond to -> weights
6163
if II[j] <= (D.boundary_point_count-1)
6264
weights = D.low_boundary_coefs[II[j]][k]
63-
offset = D.boundary_point_count - II[j] + 1
64-
# ? Is this offset correct?
65-
Itap = [II + (i+offset)*I1 for i in half_range(D.boundary_stencil_length)]
65+
offset = 1 - II[j]
66+
Itap = [II + (i+offset)*I1 for i in 0:(D.boundary_stencil_length-1)]
6667
elseif II[j] > (length(s, x) - D.boundary_point_count - 1)
6768
weights = D.high_boundary_coefs[length(s,x)-II[j]+1][k]
68-
offset = length(s, x) - II[j] - D.boundary_point_count
69-
Itap = [II + (i+offset)*I1 for i in half_range(D.boundary_stencil_length)]
69+
offset = length(s, x) - II[j]
70+
Itap = [II + (i+offset)*I1 for i in (-D.boundary_stencil_length+1):1:0]
7071
else
7172
weights = D.stencil_coefs
7273
Itap = [II + (i+1)*I1 for i in half_range(D.stencil_length)]

‎src/discretization/discretize_vars.jl‎

Lines changed: 10 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -10,7 +10,7 @@ end
1010
Return a map of of variables to the gridpoints at the edge of the domain
1111
"""
1212
@inline function get_edgevals(s, i)
13-
return [s.nottime[i] => first(s.axies[i]), s.nottime[i] => last(s.axies[i])]
13+
return [s.nottime[i] => first(s.axies[nottime[i]]), s.nottime[i] => last(s.axies[nottime[i]])]
1414
end
1515

1616
@inline function subvar(depvar, edge_vals)
@@ -19,6 +19,7 @@ end
1919
struct DiscreteSpace{N,M}
2020
vars
2121
discvars
22+
time
2223
nottime # Note that these aren't necessarily @parameters
2324
params
2425
axies
@@ -41,8 +42,8 @@ params(s::DiscreteSpace{N,M}) where {N,M}= s.params
4142
grid_idxs(s::DiscreteSpace) = CartesianIndices(((axes(g)[1] for g in s.grid)...,))
4243
edge_idxs(s::DiscreteSpace{N}) where {N} = reduce(vcat, [[vcat([Colon() for j = 1:i-1], 1, [Colon() for j = i+1:N]), vcat([Colon() for j = 1:i-1], length(s.axies[i]), [Colon() for j = i+1:N])] for i = 1:N])
4344

44-
axiesvals(s::DiscreteSpace{N}) where {N} = map(y -> [s.nottime[i] => s.axies[i][y.I[i]] for i = 1:N], s.Iaxies)
45-
gridvals(s::DiscreteSpace{N}) where {N} = map(y -> [s.nottime[i] => s.grid[i][y.I[i]] for i = 1:N], s.Igrid)
45+
axiesvals(s::DiscreteSpace{N}) where {N} = map(y -> [s.nottime[i] => s.axies[s.nottime[i]][y.I[i]] for i = 1:N], s.Iaxies)
46+
gridvals(s::DiscreteSpace{N}) where {N} = map(y -> [s.nottime[i] => s.grid[s.nottime[i]][y.I[i]] for i = 1:N], s.Igrid)
4647

4748
## Boundary methods ##
4849
edgevals(s::DiscreteSpace{N}) where {N} = reduce(vcat, [get_edgevals(s.nottime, s.axies, i) for i = 1:N])
@@ -53,8 +54,11 @@ edgevars(s::DiscreteSpace) = [[d[e...] for e in s.Iedge] for d in s.discvars]
5354
return Dict(bclocs(s) .=> [axiesvals(s)[e...] for e in s.Iedge])
5455
end
5556

57+
Iinterior(s::DiscreteSpace) = s.Igrid[[2:(length(s, x)-1) for x in s.nottime]...]
58+
5659
map_symbolic_to_discrete(II::CartesianIndex, s::DiscreteSpace{N,M}) where {N,M} = vcat([s.vars[k] => s.discvars[k][II] for k = 1:M], [s.nottime[j] => s.grid[j][II[j]] for j = 1:N])
5760

61+
# ? How rude is this? Makes Iedge work
5862

5963
# TODO: Allow other grids
6064

@@ -91,7 +95,7 @@ function DiscreteSpace(domain, depvars, indvars, nottime, discretization)
9195
# Build symbolic variables
9296
Iaxies = CartesianIndices(((axes(s.second)[1] for s in axies)...,))
9397
Igrid = CartesianIndices(((axes(g.second)[1] for g in grid)...,))
94-
@show depvars
98+
9599
depvarsdisc = map(depvars) do u
96100
if t === nothing
97101
sym = nameof(SymbolicUtils.operation(u))
@@ -104,14 +108,12 @@ function DiscreteSpace(domain, depvars, indvars, nottime, discretization)
104108
end
105109
end
106110

107-
108111
# Build symbolic maps for boundaries
109-
Iedge = reduce(vcat, [[Igrid[vcat([Colon() for j = 1:i-1], 1, [Colon() for j = i+1:nspace])],
110-
Igrid[vcat([Colon() for j = 1:i-1], length(axies[i].second), [Colon() for j = i+1:nspace])]] for i = 1:nspace])
112+
Iedge = vcat((vcat(vec(selectdim(Igrid, dim, 1)), vec(selectdim(Igrid, dim, length(grid[dim].second)))) for dim in 1:nspace)...)
111113

112114
nottime2dim = [nottime[i] => i for i in 1:nspace]
113115
dim2nottime = [i => nottime[i] for i in 1:nspace]
114-
return DiscreteSpace{nspace,length(depvars)}(depvars, Dict(depvarsdisc), nottime, indvars, Dict(axies), Dict(grid), Dict(dxs), Iaxies, Igrid, Iedge, Dict(nottime2dim))
116+
return DiscreteSpace{nspace,length(depvars)}(depvars, Dict(depvarsdisc), discretization.time, nottime, indvars, Dict(axies), Dict(grid), Dict(dxs), Iaxies, Igrid, Iedge, Dict(nottime2dim))
115117
end
116118

117119

‎src/discretization/generate_finite_difference_rules.jl‎ renamed to ‎src/discretization/generate_rules.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -149,7 +149,7 @@ There are of course more specific schemes that are used to improve stability/spe
149149
Please submit an issue if you know of any special cases that are not implemented, with links to papers and/or code that demonstrates the special case.
150150
"""
151151
function generate_finite_difference_rules(II, s, pde, derivweights)
152-
152+
@show II
153153
valrules = vcat([u => s.discvars[u][II] for u in s.vars],
154154
[x => s.grid[x][II[j]] for (j,x) in enumerate(s.nottime)])
155155
# central_deriv_rules = [(Differential(s)^2)(u) => central_deriv(2,II,j,k) for (j,s) in enumerate(s.nottime), (k,u) in enumerate(s.vars)]

‎test/runtests.jl‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -8,10 +8,10 @@ const is_TRAVIS = haskey(ENV,"TRAVIS")
88

99
@time begin
1010

11-
# if GROUP == "All" || GROUP == "Component Tests"
12-
# @time @safetestset "Test for regression against original code" begin include("regression_test.jl") end
13-
# @time @safetestset "MOLFiniteDifference Utils" begin include("utils_test.jl") end
14-
# end
11+
if GROUP == "All" || GROUP == "Component Tests"
12+
# @time @safetestset "Test for regression against original code" begin include("regression_test.jl") end
13+
#@time @safetestset "MOLFiniteDifference Utils" begin include("utils_test.jl") end
14+
end
1515

1616
if GROUP == "All" || GROUP == "Integration Tests"
1717
@time @safetestset "MOLFiniteDifference Interface" begin include("pde_systems/MOLtest.jl") end

0 commit comments

Comments
 (0)