Skip to content

Commit 8e93a88

Browse files
xtalaxxtalax
committed
Added boundary matching, began BoundaryHandler
1 parent 66ddbad commit 8e93a88

File tree

3 files changed

+106
-14
lines changed

3 files changed

+106
-14
lines changed

‎src/bcs/generate_bc_eqs.jl‎

Lines changed: 82 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,29 @@
11
### INITIAL AND BOUNDARY CONDITIONS ###
2+
3+
abstract type AbstractBoundary end
4+
5+
struct LowerBoundary <: AbstractBoundary
6+
end
7+
8+
struct UpperBoundary<: AbstractBoundary
9+
end
10+
11+
struct CompleteBoundary <: AbstractBoundary
12+
end#
13+
14+
struct PeriodicBoundary <: AbstractBoundary
15+
end
16+
17+
struct BoundaryHandler{hasperiodic}
18+
boundaries::Dict{Num, AbstractBoundary}
19+
end
20+
221
"""
3-
Mutates bceqs and u0 by finding relevant equations and discretizing them
22+
Mutates bceqs and u0 by finding relevant equations and discretizing them.
23+
TODO: return a handler for use with generate_finite_difference_rules
424
"""
5-
function generate_u0_and_bceqs!!(u0, bceqs, bcs, s, depvar_ops, tspan, derivweights)
25+
function BoundaryHandler!!(u0, bceqs, bcs, s, depvar_ops, tspan, derivweights)
26+
627
t=s.time
728

829
if t === nothing
@@ -11,29 +32,81 @@ function generate_u0_and_bceqs!!(u0, bceqs, bcs, s, depvar_ops, tspan, derivweig
1132
initmaps = substitute.(s.vars,[t=>tspan[1]])
1233
end
1334

35+
# Create some rules to match which bundary/variable a bc concerns
36+
# ? Is it nessecary to check whether all other args are present?
37+
lower_boundary_rules = vec([@rule operation(u)(~~a, lowerboundary(x), ~~b) => IfElse.ifelse(all(y-> y in vcat(~~a, ~~b), setdiff(x, arguments(u))), x, nothing) for x in setdiff(arguments(u), t), u in s.vars])
38+
39+
upper_boundary_rules = vec([@rule operation(u)(~~a, upperboundary(x), ~~b) => IfElse.ifelse(all(y-> y in vcat(~~a, ~~b), setdiff(x, arguments(u))), x, nothing) for x in setdiff(arguments(u), t), u in s.vars])
40+
1441
# Generate initial conditions and bc equations
1542
for bc in bcs
1643
bcdepvar = first(get_depvars(bc.lhs, depvar_ops))
17-
if any(u->isequal(operation(u),operation(bcdepvar)), s.vars)
44+
if any(u -> isequal(operation(u), operation(bcdepvar)), s.vars)
1845
if t !== nothing && operation(bc.lhs) isa Sym && !any(x -> isequal(x, t.val), arguments(bc.lhs))
1946
# initial condition
20-
# Assume in the form `u(...) ~ ...` for now
47+
# * Assume in the form `u(...) ~ ...` for now
48+
# * Assume that the initial condition is not in terms of the initial derivative
2149
initindex = findfirst(isequal(bc.lhs), initmaps)
2250
if initindex !== nothing
2351
push!(u0,vec(s.discvars[s.vars[initindex]] .=> substitute.((bc.rhs,),gridvals(s))))
2452
end
2553
else
26-
# boundary condition
27-
# 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
28-
push!(bceqs, vec(map(s.Iedge) do II
29-
rules = generate_finite_difference_rules(II, s, bc, derivweights)
30-
rules = vcat(rules, generate_bc_rules(II, s, bc , derivweights))
54+
# Split out additive terms
55+
rhs_arg = istree(pde.rhs) && (SymbolicUtils.operation(pde.rhs) == +) ? SymbolicUtils.arguments(pde.rhs) : [pde.rhs]
56+
lhs_arg = istree(pde.lhs) && (SymbolicUtils.operation(pde.lhs) == +) ? SymbolicUtils.arguments(pde.lhs) : [pde.lhs]
57+
58+
u_, x_ = (nothing, nothing)
59+
terms = vcat(lhs_arg,rhs_arg)
60+
# * Assume that the term of the condition is applied additively and has no multiplier/divisor/power etc.
61+
boundary = nothing
62+
# Check whether the bc is on the lower boundary, or periodic
63+
for term in terms, r in lower_boundary_rules
64+
if r(term) !== nothing
65+
u_, x_ = (term, r(term))
66+
boundary = :lower
67+
for term_ in setdiff(terms, term)
68+
for r in upper_boundary_rules
69+
if r(term_) !== nothing
70+
# boundary = :periodic
71+
#TODO: Add handling for perioodic boundary conditions here
72+
end
73+
end
74+
break
75+
end
76+
end
77+
for term in terms, r in upper_boundary_rules
78+
if r(term) !== nothing
79+
u_, x_ = (term, r(term))
80+
boundary = :upper
81+
break
82+
end
83+
end
84+
85+
@assert boundary !== nothing "Boundary condition ${bc} is not on a boundary of the domain, or is not a valid boundary condition"
86+
87+
push!(bceqs, vec(map(s.Iedge[x_][boundary]) do II
88+
rules = generate_bc_rules(II, s, bc, u_, boundary, derivweights)
89+
rules = vcat(rules, generate_finite_difference_rules(II, s, bc, derivweights))
3190

3291
substitute(bc.lhs, rules) ~ substitute(bc.rhs, rules)
3392
end))
3493
end
94+
else
95+
throw(ArgumentError("No active variables in boundary condition $bc lhs, please ensure that bcs are "))
96+
end
97+
end
98+
99+
function get_active_variable(bc, s, depvar_ops)
100+
bcdepvar = first(get_depvars(bc.lhs, depvar_ops))
101+
out = Array{typeof(operation(s.vars[1]))}()
102+
for u in s.vars
103+
if u isa Sym && isequal(operation(u), operation(bcdepvar))
104+
push!(out, u)
35105
end
36106
end
107+
return out
37108
end
109+
110+
38111

39112
#function generate_u0_and_bceqs_with_rules!!(u0, bceqs, bcs, t, s, depvar_ops)

‎src/discretization/differential_discretizer.jl‎

Lines changed: 6 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,8 @@
22
struct DifferentialDiscretizer{T, D1}
33
approx_order::Int
44
map::D1
5-
nonlinlapmap::Dict{Num, NTuple{2, DiffEqOperators.DerivativeOperator}}
5+
halfoffsetmap::Dict{Num, DiffEqOperators.DerivativeOperator}
6+
interpmap::Dict{Num, DiffEqOperators.DerivativeOperator}
67
orders::Dict{Num, Vector{Int}}
78
end
89

@@ -13,6 +14,7 @@ function DifferentialDiscretizer(pde, s, discretization)
1314
# 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)]
1415
differentialmap = Array{Pair{Num,DiffEqOperators.DerivativeOperator},1}()
1516
nonlinlap = []
17+
interp = []
1618
orders = []
1719
# Hardcoded to centered difference, generate weights for each differential
1820
# TODO: Add handling for upwinding
@@ -22,10 +24,11 @@ function DifferentialDiscretizer(pde, s, discretization)
2224
rs = [(Differential(x)^d) => CompleteCenteredDifference(d, approx_order, s.dxs[x] ) for d in last(orders).second]
2325

2426
differentialmap = vcat(differentialmap, rs)
25-
push!(nonlinlap, x => (CompleteHalfCenteredDifference(0, approx_order, s.dxs[x]), CompleteHalfCenteredDifference(1, approx_order, s.dxs[x])))
27+
push!(nonlinlap, x => CompleteHalfCenteredDifference(0, approx_order, s.dxs[x])
28+
push!(interp, x => CompleteHalfCenteredDifference(1, approx_order, s.dxs[x]))
2629
end
2730

28-
return DifferentialDiscretizer{eltype(orders), typeof(Dict(differentialmap))}(approx_order, Dict(differentialmap), Dict(nonlinlap), Dict(orders))
31+
return DifferentialDiscretizer{eltype(orders), typeof(Dict(differentialmap))}(approx_order, Dict(differentialmap), Dict(nonlinlap), Dict(interp), Dict(orders))
2932
end
3033

3134

‎src/discretization/generate_rules.jl‎

Lines changed: 18 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -47,7 +47,8 @@ function cartesian_nonlinear_laplacian(expr, II, derivweights, s, x, u)
4747
# See scheme 1, namely the term without the 1/r dependence. See also #354 and #371 in DiffEqOperators, the previous home of this package.
4848

4949
jx = j, x = (s.x2i(x), x)
50-
inner_interpolater, D_inner = derivweights.nonlinlapmap[x]
50+
D_inner = derivweights.halfoffsetmap[x]
51+
inner_interpolater = derivweights.interpmap[x]
5152

5253
# Get the outer weights and stencil to generate the required
5354
outerweights, outerstencil = get_half_offset_weights_and_stencil(D_inner, II, s, 0, jx)
@@ -149,7 +150,6 @@ There are of course more specific schemes that are used to improve stability/spe
149150
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.
150151
"""
151152
function generate_finite_difference_rules(II, s, pde, derivweights)
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)]
@@ -205,3 +205,19 @@ function generate_finite_difference_rules(II, s, pde, derivweights)
205205
return rules
206206
end
207207

208+
function generate_bc_rules(II, dim, s, bc, edgemaps)
209+
# ! Recognise which dim a BC is on, and use that to get the axiesvals at the boundary
210+
# ! loop through the Iedge at that boundary
211+
# ! replace the symbolic variables at either end of the boundary with the appropriate discvars, interpolated if nessecary
212+
# ! eventually move to multi dim interpolations to improve validity
213+
214+
# depvarbcmaps will dictate what to replace the variable terms with in the bcs
215+
# replace u(t,0) with u₁, etc
216+
if grid_align == center_align
217+
depvarbcmaps = reduce(vcat,[substitute(depvar, edgevals(s)) .=> edgevar for (depvar, edgevar) in zip(s.vars, edgevars(s, II))])
218+
elseif grid_align == edge_align
219+
220+
end
221+
222+
varrules = edgemaps(s)
223+
rules = vcat(depvarbcmaps, edgemaps)

0 commit comments

Comments
 (0)