Skip to content

Commit c981b3b

Browse files
xtalaxxtalax
committed
Moved around some comments
1 parent 67474ef commit c981b3b

File tree

2 files changed

+166
-9
lines changed

2 files changed

+166
-9
lines changed

‎src/discretization/MOL_discretization.jl‎

Lines changed: 1 addition & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -87,16 +87,8 @@ function SciMLBase.symbolic_discretize(pdesys::PDESystem, discretization::Method
8787
end
8888
return [left, right]
8989
end
90-
#--------------------------------------------------
91-
# * The stencil is the tappoints of the finite difference operator, relative to the current index
92-
stencil(j, order) = CartesianIndices(Tuple(map(x -> -x:x, (1:nspace.==j) * (order÷2))))
90+
9391

94-
# TODO: Generalize central difference handling to allow higher even order derivatives
95-
# The central neighbour indices should add the stencil to II, unless II is too close
96-
# to an edge in which case we need to shift away from the edge
97-
98-
99-
# Calculate buffers
10092
#TODO: Update this when allowing higher order derivatives to correctly deal with boundary upwinding
10193
interior = s.Igrid[[let bcs = get_bc_counts(i)
10294
(1 + first(bcs)):length(g)-last(bcs)

‎src/discretization/generate_finite_difference_rules_new.jl‎

Lines changed: 165 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,165 @@
1+
function calculate_weights_cartesian(order::Int, x0::T, xs::AbstractVector, idxs::AbstractVector) where T<:Real
2+
# Cartesian domain: use Fornberg
3+
calculate_weights(order, x0, vec(xs[idxs]))
4+
end
5+
function calculate_weights_spherical(order::Int, x0::T, x::AbstractVector, idxs::AbstractVector) where T<:Real
6+
# Spherical domain: see #367
7+
# https://web.mit.edu/braatzgroup/analysis_of_finite_difference_discretization_schemes_for_diffusion_in_spheres_with_variable_diffusivity.pdf
8+
# Only order 2 is implemented
9+
@assert order == 2
10+
# Only 2nd order discretization is implemented
11+
# We can't activate this assertion for now because the rules try to create the spherical Laplacian
12+
# before checking whether there is a spherical Laplacian
13+
# this could be fixed by dispatching on domain type when we have different domain types
14+
# but for now everything is an Interval
15+
# @assert length(x) == 3
16+
# TODO: nonlinear diffusion in a spherical domain
17+
i = idxs[2]
18+
dx1 = x[i] - x[i-1]
19+
dx2 = x[i+1] - x[i]
20+
i0 = i - 1 # indexing starts at 0 in the paper and starts at 1 in julia
21+
1 / (i0 * dx1 * dx2) * [i0-1, -2i0, i0+1]
22+
end
23+
24+
function generate_finite_difference_rules(II, s, pde, discretization)
25+
approx_order = discretization.centered_order
26+
I1 = oneunit(first(s.Igrid))
27+
28+
Imin(I1, order) = first(s.Igrid) + I1 * (order ÷ 2)
29+
Imax(I1, order) = last(s.Igrid) - I1 * (order ÷ 2)
30+
31+
#--------------------------------------------------
32+
# * The stencil is the tappoints of the finite difference operator, relative to the current index
33+
stencil(j, order) = CartesianIndices(Tuple(map(x -> -x:x, (1:length(s.nottime) .== j) * (order ÷ 2))))
34+
35+
# TODO: Generalize central difference handling to allow higher even order derivatives
36+
# The central neighbour indices should add the stencil to II, unless II is too close
37+
# to an edge in which case we need to shift away from the edge
38+
39+
central_neighbor_idxs(II, j, order) = stencil(j, order) .+ max(Imin(I1, order), min(II, Imax(I1, order)))
40+
central_weights_cartesian(d_order, II, j) = calculate_weights_cartesian(d_order, s.grid[j][II[j]], s.grid[j], vec(map(i -> i[j],
41+
central_neighbor_idxs(II, j, approx_order))))
42+
central_deriv(d_order, II, j, k) = dot(central_weights(d_order, II, j), s.discvars[k][central_neighbor_idxs(II, j, approx_order)])
43+
44+
central_deriv_cartesian(d_order, II, j, k) = dot(central_weights_cartesian(d_order, II, j), s.discvars[k][central_neighbor_idxs(II, j, approx_order)])
45+
46+
# spherical Laplacian has a hardcoded order of 2 (only 2nd order is implemented)
47+
# both for derivative order and discretization order
48+
central_weights_spherical(II, j) = calculate_weights_spherical(2, s.grid[j][II[j]], s.grid[j], vec(map(i -> i[j], central_neighbor_idxs(II, j, 2))))
49+
central_deriv_spherical(II, j, k) = dot(central_weights_spherical(II, j), s.discvars[k][central_neighbor_idxs(II, j, 2)])
50+
51+
# get a sorted list derivative order such that highest order is first. This is useful when substituting rules
52+
# starting from highest to lowest order.
53+
d_orders(x) = reverse(sort(collect(union(differential_order(pde.rhs, x), differential_order(pde.lhs, x)))))
54+
55+
# orders = [d_orders(s) for s in s.nottime]
56+
# D = [CenteredDifference{1} for i in 1:length(s.nottime), d in dorders()]
57+
58+
# 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)]
59+
central_deriv_rules_cartesian = Array{Pair{Num,Num},1}()
60+
for (j, x) in enumerate(s.nottime)
61+
rs = [(Differential(x)^d)(u) => central_deriv_cartesian(d, II, j, k) for d in d_orders(x), (k, u) in enumerate(s.vars)]
62+
for rule in rs
63+
push!(central_deriv_rules_cartesian, rule)
64+
end
65+
end
66+
67+
central_deriv_rules_spherical = [Differential(r)(r^2 * Differential(r)(u)) / r^2 => central_deriv_spherical(II, j, k)
68+
for (j, r) in enumerate(s.nottime), (k, u) in enumerate(s.vars)]
69+
70+
valrules = map_symbolic_to_discrete(II, s)
71+
72+
# ! Use DerivativeOperator to get the coefficients, write a new constructor
73+
# TODO: upwind rules needs interpolation into `@rule`
74+
# upwinding_rules = [@rule(*(~~a,$(Differential(s.nottime[j]))(u),~~b) => IfElse.ifelse(*(~~a..., ~~b...,)>0,
75+
# *(~~a..., ~~b..., dot(reverse_weights(II,j),s.vars[k][central_neighbor_idxs(II,j)[1:2]])),
76+
# *(~~a..., ~~b..., dot(forward_weights(II,j),s.vars[k][central_neighbor_idxs(II,j)[2:3]]))))
77+
# for j in 1:length(s.nottime), k in 1:length(pdesys.s.vars)]
78+
79+
80+
return vcat(nonlinear_laplacian_rules(II, I1, pde, s, discretization, central_deriv_spherical),
81+
vec(central_deriv_rules_cartesian),
82+
vec(central_deriv_rules_spherical),
83+
valrules)
84+
end
85+
86+
function nonlinear_laplacian_rules(II, I1, pde, s, discretization, central_deriv_spherical)
87+
approx_order = discretization.centered_order
88+
# # Discretization of non-linear laplacian.
89+
# d/dx( a du/dx ) ~ (a(x+1/2) * (u[i+1] - u[i]) - a(x-1/2) * (u[i] - u[i-1]) / dx^2
90+
stencil(j, order) = CartesianIndices(Tuple(map(x -> -x:x, (1:length(s.nottime) .== j) * (order ÷ 2))))
91+
Imin(I1, order) = first(s.Igrid) + I1 * (order ÷ 2)
92+
Imax(I1, order) = last(s.Igrid) - I1 * (order ÷ 2)
93+
94+
central_neighbor_idxs(II, j, order) = stencil(j, order) .+ max(Imin(I1, order), min(II, Imax(I1, order)))
95+
# TODO: upwind rules needs interpolation into `@rule`
96+
97+
98+
forward_weights(II, j) = calculate_weights(discretization.upwind_order, 0.0, s.grid[j][[II[j], II[j] + 1]])
99+
reverse_weights(II, j) = -calculate_weights(discretization.upwind_order, 0.0, s.grid[j][[II[j] - 1, II[j]]])
100+
101+
102+
reverse_finite_difference(II, j, k) = -dot(reverse_weights(II, j), s.discvars[k][central_neighbor_idxs(II, j, approx_order)[1:2]]) / s.dxs[j]
103+
forward_finite_difference(II, j, k) = dot(forward_weights(II, j), s.discvars[k][central_neighbor_idxs(II, j, approx_order)[2:3]]) / s.dxs[j]
104+
105+
# TODO: improve interpolation of g(x) = u(x) for calculating u(x+-dx/2)
106+
interpolate_discrete_depvar(II, j, k, l) = sum([s.discvars[k][central_neighbor_idxs(II, j, approx_order)][i] for i in (l == 1 ? [2, 3] : [1, 2])]) / 2.0
107+
# interpolate_discrete_indvar returns middle space values. E.g. x(i-1/2) or y(i+1/2).
108+
interpolate_discrete_indvar(II, j, l) = (s.grid[j][II[j]] + s.grid[j][II[j]+l]) / 2.0
109+
# Dependent variable rules
110+
map_depvars_to_discrete(II, j, l) = [s.vars[k] => interpolate_discrete_depvar(II, j, k, l) for k = 1:length(s.vars)]
111+
# Independent variable rules
112+
map_indvars_to_discrete(II, l) = [s.nottime[j] => interpolate_discrete_indvar(II, j, l) for j = 1:length(s.nottime)]
113+
114+
build_discrete_symbolic_expression(II, expr, j, l) = Num(substitute(substitute(expr, map_depvars_to_discrete(II, j, l)), map_indvars_to_discrete(II, l)))
115+
116+
function cartesian_deriv(expr, II, j, k)
117+
discrete_expression_upwind = build_discrete_symbolic_expression(II, expr, j, 1)
118+
discrete_expression_downwind = build_discrete_symbolic_expression(II, expr, j, -1)
119+
# We need to do
120+
discrete_expression_combined = [discrete_expression_downwind, discrete_expression_upwind]
121+
combined_finite_difference = [reverse_finite_difference(II, j, k), forward_finite_difference(II, j, k)]
122+
return dot(discrete_expression_combined, combined_finite_difference)
123+
end
124+
125+
function generate_cartesian_deriv_rules()
126+
rules = [@rule ($(Differential(iv))(*(~~a, $(Differential(iv))(dv), ~~b))) =>
127+
cartesian_deriv(*(~~a..., ~~b...), II, j, k) for (j, iv) in enumerate(s.nottime) for (k, dv) in enumerate(s.vars)]
128+
129+
rules = vcat(vec(rules), vec(
130+
[@rule ($(Differential(iv))($(Differential(iv))(dv) / ~a)) =>
131+
cartesian_deriv(1 / ~a, II, j, k) for (j, iv) in enumerate(s.nottime) for (k, dv) in enumerate(s.vars)]))
132+
133+
return rules
134+
end
135+
136+
function generate_spherical_deriv_rules()
137+
rules = [@rule *(~~a, ($(Differential(iv))((iv^2) * $(Differential(iv))(dv))), ~~b) / (iv^2) =>
138+
*(~a..., central_deriv_spherical(II, j, k), ~~b...)
139+
for (j, iv) in enumerate(s.nottime) for (k, dv) in enumerate(s.vars)]# r^-2 needs to be handled separately
140+
rules = vcat(vec(rules), vec(
141+
[@rule *(~~a, (iv^-2) * ($(Differential(iv))((iv^2) * $(Differential(iv))(dv))), ~~b) =>
142+
*(~a..., central_deriv_spherical(II, j, k), ~~b...) for (j, iv) in enumerate(s.nottime) for (k, dv) in enumerate(s.vars)]))
143+
return rules
144+
end
145+
146+
cartesian_deriv_rules = vec(generate_cartesian_deriv_rules())
147+
spherical_deriv_rules = vec(generate_spherical_deriv_rules())
148+
149+
# If lhs or rhs is a linear combination, apply rules term by term
150+
rhs_arg = istree(pde.rhs) && (SymbolicUtils.operation(pde.rhs) == +) ? SymbolicUtils.arguments(pde.rhs) : [pde.rhs]
151+
lhs_arg = istree(pde.lhs) && (SymbolicUtils.operation(pde.lhs) == +) ? SymbolicUtils.arguments(pde.lhs) : [pde.lhs]
152+
153+
nonlinlap_rules = []
154+
155+
for term in vcat(lhs_arg, rhs_arg)
156+
for rule in vcat(cartesian_deriv_rules, spherical_deriv_rules)
157+
if rule(term) !== nothing
158+
push!(nonlinlap_rules, term => rule(term))
159+
end
160+
end
161+
end
162+
163+
164+
return vec(nonlinlap_rules)
165+
end

0 commit comments

Comments
 (0)