Skip to content

Commit 47f0659

Browse files
xtalaxxtalax
committed
stash
1 parent d364719 commit 47f0659

File tree

6 files changed

+333
-327
lines changed

6 files changed

+333
-327
lines changed

‎src/discretization/MOL_discretization.jl‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -68,11 +68,11 @@ function SciMLBase.symbolic_discretize(pdesys::PDESystem, discretization::Method
6868
indvars = first(allindvars)
6969
nspace = length(nottime)
7070

71-
s = DiscreteSpace(pdesys.domain, depvars, indvars, nottime, grid_align, discretization)
71+
s = DiscreteSpace(pdesys.domain, depvars, indvars, nottime, discretization)
7272

7373
derivweights = DifferentialDiscretizer(pde, s, discretization)
7474

75-
interior = s.Igrid[[2:(length(grid[x])-1) for x in s.nottime]
75+
interior = s.Igrid[[2:(length(grid[x])-1) for x in s.nottime]]
7676

7777
### PDE EQUATIONS ###
7878
# Create a stencil in the required dimension centered around 0

‎src/discretization/generate_finite_difference_rules.jl‎

Lines changed: 166 additions & 85 deletions
Original file line numberDiff line numberDiff line change
@@ -1,108 +1,188 @@
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]))
1+
"""
2+
`interpolate_discrete_param`
3+
4+
Interpolate gridpoints by taking the average of the values of the discrete points, or if the offset is outside the grid, extrapolate the value with dx.
5+
"""
6+
function interpolate_discrete_param(II, s, itap, j, x)
7+
# * This will need to be updated to dispatch on grid type when grids become more general
8+
offset = itap+1/2
9+
if (II[j]+itap) < 1
10+
return s.grid[x][1]+s.dxs[x]*offset
11+
elseif (II[j]+itap) > (length(x) - 1)
12+
return s.grid[x][length(x)]+s.dxs[x]*offset
13+
else
14+
return (s.grid[x][II[j]+itap]+s.grid[x][II[j]+itap+1])/2
15+
end
416
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]
17+
18+
"""
19+
`cartesian_nonlinear_laplacian`
20+
21+
Differential(x)(expr(x)*Differential(x)(u(x)))
22+
23+
Given an internal multiplying expression `expr`, return the correct finite difference equation for the nonlinear laplacian at the location in the grid given by `II`.
24+
25+
The inner derivative is discretized with the half offset centered scheme, giving the derivative at interpolated grid points offset by dx/2 from the regular grid.
26+
27+
The outer derivative is discretized with the centered scheme, giving the nonlinear laplacian at the grid point `II`.
28+
For first order returns something like this:
29+
`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`
30+
31+
For 4th order, returns something like this:
32+
```
33+
first_finite_diffs = [a(x-3/2)*finitediff(u, i-3/2),
34+
a(x-1/2)*finitediff(u, i-1/2),
35+
a(x+1/2)*finitediff(u, i+1/2),
36+
a(x+3/2)*finitediff(u, i+3/2)]
37+
38+
dot(central_finite_diff_weights, first_finite_diffs)
39+
```
40+
41+
where `finitediff(u, i)` is the finite difference at the interpolated point `i` in the grid.
42+
43+
And so on.
44+
"""
45+
function cartesian_nonlinear_laplacian(expr, II, derivweights, s, x, u)
46+
# Based on the paper https://web.mit.edu/braatzgroup/analysis_of_finite_difference_discretization_schemes_for_diffusion_in_spheres_with_variable_diffusivity.pdf
47+
# See scheme 1, namely the term without the 1/r dependence. See also #354 and #371 in DiffEqOperators, the previous home of this package.
48+
49+
jx = j, x = (s.x2i(x), x)
50+
inner_interpolater, D_inner = derivweights.nonlinlapmap[x]
51+
52+
# Get the outer weights and stencil to generate the required
53+
outerweights, outerstencil = get_half_offset_weights_and_stencil(D_inner, II, s, 0, jx)
54+
# Index offsets of each stencil in the inner finite difference to get the correct stencil for each needed half grid point, 0 corresopnds to x+1/2
55+
itaps = getindex.(outerstencil, (j,))
56+
57+
# Get the correct weights and stencils for this II
58+
inner_deriv_weights_and_stencil = [get_half_offset_weights_and_stencil(D_inner, II, s, itap, jx) for itap in itaps]
59+
interp_weights_and_stencil = [get_half_offset_weights_and_stencil(inner_interpolater, II, s, itap, jx) for itap in itaps]
60+
61+
# map variables to symbolically inerpolated/extrapolated expressions
62+
map_vars_to_interpolated(stencil, weights) = [v => dot(weights, s.discvars[v][stencil]) for v in s.vars]
63+
64+
# Map parameters to interpolated values. Using simplistic extrapolation/interpolation for now as grids are uniform
65+
map_params_to_interpolated(itap) = [z => interpolate_discrete_param(II, s, itap, i, z) for (i,z) in s.nottime]
66+
67+
# Take the inner finite difference
68+
inner_difference = [dot(inner_weights, s.discvars[u][inner_stencil]) for (inner_weights, inner_stencil) in inner_deriv_weights_and_stencil]
69+
70+
# Symbolically interpolate the multiplying expression
71+
interpolated_expr = [Num(substitute(substitute(expr, map_vars_to_interpolated(interpstencil, interpweights)), map_params_to_interpolated(itap)))
72+
for (itap,(interpweights, interpstencil)) in zip(itaps, interp_weights_and_stencil)]
73+
74+
# multiply the inner finite difference by the interpolated expression, and finally take the outer finite difference
75+
return dot(weights, inner_difference .* interpolated_expr)
2276
end
2377

24-
function generate_finite_difference_rules(II, s, pde, discretization)
25-
approx_order = discretization.centered_order
26-
I1 = oneunit(first(s.Igrid))
27-
Imin(order) = first(s.Igrid) + I1 * (order ÷ 2)
28-
Imax(order) = last(s.Igrid) - I1 * (order ÷ 2)
29-
stencil(j, order) = CartesianIndices(Tuple(map(x -> -x:x, (1:length(s.nottime) .== j) * (order ÷ 2))))
30-
# Use max and min to apply buffers
31-
central_neighbor_idxs(II,j,order) = stencil(j,order) .+ max(Imin(order),min(II,Imax(order)))
32-
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],central_neighbor_idxs(II,j,approx_order))))
33-
central_deriv(d_order, II,j,k) = dot(central_weights(d_order, II,j),s.discvars[k][central_neighbor_idxs(II,j,approx_order)])
78+
"""
79+
`cartesian_nonlinear_laplacian`
80+
81+
Differential(x)(expr(x)*Differential(x)(u(x)))
82+
83+
Given an internal multiplying expression `expr`, return the correct finite difference equation for the nonlinear laplacian at the location in the grid given by `II`.
84+
85+
The inner derivative is discretized with the half offset centered scheme, giving the derivative at interpolated grid points offset by dx/2 from the regular grid.
86+
87+
The outer derivative is discretized with the centered scheme, giving the nonlinear laplacian at the grid point `II`.
88+
For first order returns something like this:
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+
91+
For 4th order, returns something like this:
92+
```
93+
first_finite_diffs = [a(x-3/2)*finitediff(u, i-3/2),
94+
a(x-1/2)*finitediff(u, i-1/2),
95+
a(x+1/2)*finitediff(u, i+1/2),
96+
a(x+3/2)*finitediff(u, i+3/2)]
97+
98+
dot(central_finite_diff_weights, first_finite_diffs)
99+
```
34100
35-
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)])
101+
where `finitediff(u, i)` is the finite difference at the interpolated point `i` in the grid.
36102
37-
# spherical Laplacian has a hardcoded order of 2 (only 2nd order is implemented)
38-
# both for derivative order and discretization order
39-
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))))
40-
central_deriv_spherical(II,j,k) = dot(central_weights_spherical(II,j),s.discvars[k][central_neighbor_idxs(II,j,2)])
103+
And so on.
104+
"""
105+
function spherical_diffusion(innerexpr, II, derivweights, s, r, u)
106+
# Based on the paper https://web.mit.edu/braatzgroup/analysis_of_finite_difference_discretization_schemes_for_diffusion_in_spheres_with_variable_diffusivity.pdf
107+
D_1 = derivweights.map[Differential(r)]
108+
D_2 = derivweights.map[Differential(r)^2]
41109

42-
# get a sorted list derivative order such that highest order is first. This is useful when substituting rules
43-
# starting from highest to lowest order.
44-
d_orders(order) = reverse(sort(collect(union(differential_order(pde.rhs, order), differential_order(pde.lhs, order)))))
110+
# What to replace parameter x with given I
111+
_rsubs(x, I) = x => s.grid[x][I[s.x2i(x)]]
112+
# Full rules for substituting parameters in the inner expression
113+
rsubs(I) = vcat([v => s.discvars[v][I] for v in s.vars], [_rsubs(x, I) for x in params(s)])
114+
# Discretization func for u
115+
ufunc_u(v, I, x) = s.discvars[v][I]
45116

117+
# 2nd order finite difference in u
118+
exprhere = substitute(innerexpr, rsubs(II))
119+
# Catch the r ≈ 0 case
120+
if substitute(r, _rsubs(r, II)) ≈ 0
121+
D_2_u = central_difference(D_2, II, s, (s.x2i(r), r), u, ufunc_u)
122+
return 3exprhere*D_2_u # See appendix B of the paper
123+
124+
D_1_u = central_difference(D_1, II, s, (s.x2i[r], r), u, ufunc_u)
125+
# See scheme 1 in appendix A of the paper
126+
127+
return exprhere*(D_1_u/substitute(r, _rsubs(r, II)) + cartesian_nonlinear_laplacian(innerexpr, II, derivweights, s, r, u))
128+
end
129+
130+
"""
131+
`generate_finite_difference_rules`
132+
133+
Generate a vector of finite difference rules to dictate what to replace variables in the `pde` with at the gridpoint `II`.
134+
135+
Care is taken to make sure that the rules only use points that are actually in the discretized grid by progressively up/downwinding the stencils when the gridpoint `II` is close to the boundary.
136+
137+
There is a genral catch all ruleset that uses the cartesian centered difference scheme for derivatives, and simply the discretized variable at the given gridpoint for particular variables.
138+
139+
There are of course more specific schemes that are used to improve stability/speed/accuracy when particular forms are encountered in the PDE. These rules are applied first to override the general ruleset.
140+
141+
##Currently implemented special cases are as follows:
142+
- Spherical derivatives
143+
- Nonlinear laplacian uses a half offset centered scheme for the inner derivative to improve stability
144+
- Spherical nonlinear laplacian.
145+
146+
##Planned special cases include:
147+
- Up/Downwind schemes to be used for odd ordered derivatives multiplied by a coefficient, downwinding when the coefficient is positive, and upwinding when the coefficient is negative.
148+
149+
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.
150+
"""
151+
function generate_finite_difference_rules(II, s, pde, derivweights)
152+
153+
valrules = vcat([u => s.discvars[u][II] for u in s.vars],
154+
[x => s.grid[x][II[s.x2i[x]]] for x in params(s)])
46155
# 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)]
47-
central_deriv_rules_cartesian = Array{Pair{Num,Num},1}()
48-
for (j,x) in enumerate(s.nottime)
49-
rs = [(Differential(x)^d)(u) => central_deriv_cartesian(d,II,j,k) for d in d_orders(x), (k,u) in enumerate(s.vars)]
50-
for r in rs
51-
push!(central_deriv_rules_cartesian, r)
52-
end
53-
end
54156

55-
central_deriv_rules_spherical = [Differential(x)(x^2*Differential(x)(u))/x^2 => central_deriv_spherical(II,j,k)
56-
for (j,x) in enumerate(s.nottime), (k,u) in enumerate(s.vars)]
157+
central_ufunc(u, I, x) = s.discvars[u][I]
158+
central_deriv_rules_cartesian = Array{Pair{Num,Num},1}()
159+
for x in s.nottime
160+
j = s.x2i[x]
161+
rs = [(Differential(x)^d)(u) => central_difference(derivweights.map[Differential(x)^d], II, s, (j,x), u, central_ufunc) for d in derivweights.orders[x], u in s.vars]
57162

58-
valrules = vcat([s.vars[k] => s.discvars[k][II] for k in 1:length(s.vars)],
59-
[s.nottime[j] => s.grid[j][II[j]] for j in 1:nparams(s)])
163+
central_deriv_rules_cartesian = vcat(central_deriv_rules_cartesian, rs)
164+
end
60165

61166
# TODO: upwind rules needs interpolation into `@rule`
62-
forward_weights(II,j) = calculate_weights(discretization.upwind_order, 0.0, s.grid[j][[II[j],II[j]+1]])
63-
reverse_weights(II,j) = calculate_weights(discretization.upwind_order, 0.0, s.grid[j][[II[j]-1,II[j]]])
64-
# upwinding_rules = [@rule(*(~~a,$(Differential(s.nottime[j]))(u),~~b) => IfElse.ifelse(*(~~a..., ~~b...,)>0,
167+
#forward_weights(II,j) = calculate_weights(discretization.upwind_order, 0.0, s.grid[j][[II[j],II[j]+1]])
168+
#reverse_weights(II,j) = calculate_weights(discretization.upwind_order, 0.0, s.grid[j][[II[j]-1,II[j]]])
169+
# upwinding_rules = [@rule(*(~~a, $(Differential(s.nottime[j]))(u),~~b) => IfElse.ifelse(*(~~a..., ~~b...,)>0,
65170
# *(~~a..., ~~b..., dot(reverse_weights(II,j),s.vars[k][central_neighbor_idxs(II,j)[1:2]])),
66171
# *(~~a..., ~~b..., dot(forward_weights(II,j),s.vars[k][central_neighbor_idxs(II,j)[2:3]]))))
67172
# for j in 1:nparams(s), k in 1:length(pdesys.s.vars)]
68173

69174
## Discretization of non-linear laplacian.
70-
# 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
71-
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]
72-
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]
73-
# TODO: improve interpolation of g(x) = u(x) for calculating u(x+-dx/2)
74-
interpolate_discrete_depvar(II, s, 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.
75-
# iv_mid returns middle space values. E.g. x(i-1/2) or y(i+1/2).
76-
interpolate_discrete_indvar(II, j, l) = (s.grid[j][II[j]] + s.grid[j][II[j]+l]) / 2.0
77-
# Dependent variable rules
78-
map_vars_to_discrete(II, j, k, l) = [s.vars[k] => interpolate_discrete_depvar(II, s, j, k, l) for k in 1:length(s.vars)]
79-
# Independent variable rules
80-
map_params_to_discrete(II, j, l) = [s.nottime[j] => interpolate_discrete_indvar(II, j, l) for j in 1:length(s.nottime)]
81-
# Replacement rules: new approach
82-
# Calc
83-
function discrete_cartesian(expr, j, k)
84-
disc_downwind = Num(substitute(substitute(expr, map_vars_to_discrete(II, j, k, -1)), map_params_to_discrete(II, j, -1)))
85-
disc_upwind = Num(substitute(substitute(expr, map_vars_to_discrete(II, j, k, 1)), map_params_to_discrete(II, j, 1)))
86-
return dot([disc_downwind, disc_upwind],
87-
[reverse_finite_difference(II, j, k), forward_finite_difference(II, j, k)])
88-
end
89175

90-
cartesian_deriv_rules = [@rule ($(Differential(iv))(*(~~a, $(Differential(iv))(dv), ~~b))) => discrete_cartesian(*(a..., b...),j,k) for (j, iv) in enumerate(s.nottime) for (k, dv) in enumerate(s.vars)]
176+
91177

92-
cartesian_deriv_rules = vcat(vec(cartesian_deriv_rules),vec(
93-
[@rule ($(Differential(iv))($(Differential(iv))(dv)/~a)) =>
94-
discrete_cartesian(1/~a,j,k)
95-
for (j, iv) in enumerate(s.nottime) for (k, dv) in enumerate(s.vars)]))
178+
cartesian_deriv_rules = [@rule ($(Differential(x))(*(~~a, $(Differential(x))(u), ~~b))) => cartesian_nonlinear_laplacian(*(a..., b...), II, derivweights, s, x, u) for x in s.nottime, u in s.vars]
96179

97-
spherical_deriv_rules = [@rule *(~~a, ($(Differential(iv))((iv^2)*$(Differential(iv))(dv))), ~~b) / (iv^2) =>
98-
*(~a..., central_deriv_spherical(II, j, k), ~b...)
99-
for (j, iv) in enumerate(s.nottime) for (k, dv) in enumerate(s.vars)]
180+
cartesian_deriv_rules = vcat(vec(cartesian_deriv_rules),vec(
181+
[@rule ($(Differential(x))($(Differential(x))(u)/~a)) => cartesian_nonlinear_laplacian(1/~a, II, derivweights, s, x, u) for x in s.nottime, u in s.vars]))
100182

101-
# r^-2 needs to be handled separately
102-
spherical_deriv_rules = vcat(vec(spherical_deriv_rules),vec(
103-
[@rule *(~~a, (iv^-2) * ($(Differential(iv))((iv^2)*$(Differential(iv))(dv))), ~~b) =>
104-
*(~a..., central_deriv_spherical(II, j, k), ~b...)
105-
for (j, iv) in enumerate(s.nottime) for (k, dv) in enumerate(s.vars)]))
183+
spherical_deriv_rules = [@rule *(~~a, (r^-2), ($(Differential(r))(*(~~c, (r^2), ~~d, $(Differential(r))(u), ~~e))), ~~b) =>
184+
*(~a..., spherical_diffusion(*(~c..., ~d..., ~e...), II, derivweights, s, r, u), ~b...)
185+
for r in s.nottime, u in s.vars]
106186

107187
rhs_arg = istree(pde.rhs) && (SymbolicUtils.operation(pde.rhs) == +) ? SymbolicUtils.arguments(pde.rhs) : [pde.rhs]
108188
lhs_arg = istree(pde.lhs) && (SymbolicUtils.operation(pde.lhs) == +) ? SymbolicUtils.arguments(pde.lhs) : [pde.lhs]
@@ -125,4 +205,5 @@ function generate_finite_difference_rules(II, s, pde, discretization)
125205
vec(central_deriv_rules_spherical),
126206
valrules)
127207
return rules
128-
end
208+
end
209+

0 commit comments

Comments
 (0)