Skip to content
This repository was archived by the owner on Jul 19, 2023. It is now read-only.

[WIP] Progress on upwindig rules (issue 355) #390

Closed
wants to merge 17 commits into from
Closed

[WIP] Progress on upwindig rules (issue 355) #390

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
f24dbbe
start fixing issue 355
emmanuellujan May 11, 2021
c8dac10
Enable linear convection tests.
emmanuellujan May 11, 2021
062a12e
small change
emmanuellujan May 11, 2021
5391309
small change
emmanuellujan May 12, 2021
a3d9bdf
merge
emmanuellujan Jun 6, 2021
0e26e4f
Add comment in MOLtest.jl execution line. This test is not ready to b…
emmanuellujan Jun 6, 2021
2984230
IntervalDomain -> Interval
emmanuellujan Jun 6, 2021
bd64ad6
issue with @assert tdomain.domain isa IntervalDomain in MOL_discretiz…
emmanuellujan Jun 7, 2021
afc6410
.
emmanuellujan Jun 7, 2021
71880e5
removing comments
emmanuellujan Jun 8, 2021
7a8bcc7
addressing Interval issues
emmanuellujan Jun 8, 2021
f4d79fb
testing if checks work in github removing Interval from MOL_discretiz…
emmanuellujan Jun 9, 2021
8fbfaa3
Adding 2D linear convection-diffusion test case.
emmanuellujan Jun 10, 2021
5f58b05
2D linear convection-diffusion test requires discretize first derivat…
emmanuellujan Jun 10, 2021
10142b1
Progress in upwind autodiscretizatio. Issues with : -Dx(v(t,x)) * u(t,x)
emmanuellujan Jun 19, 2021
45485da
solving conflicts
emmanuellujan Jul 8, 2021
0dcc323
bound MTK to <5.20
emmanuellujan Jul 8, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ DomainSets = "0.5"
ForwardDiff = "0.10"
LazyArrays = "0.17, 0.18, 0.19, 0.20, 0.21"
LazyBandedMatrices = "0.3, 0.4, 0.5"
ModelingToolkit = "4, 5"
ModelingToolkit = "4 - 5.20"
NNlib = "0.6, 0.7"
NonlinearSolve = "0.3.7"
RuntimeGeneratedFunctions = "0.4, 0.5"
Expand Down
77 changes: 55 additions & 22 deletions src/MOLFiniteDifference/MOL_discretization.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
#using ModelingToolkit: operation, istree, arguments, Interval, infimum, supremum
using ModelingToolkit: operation, istree, arguments
import DomainSets

Expand Down Expand Up @@ -91,7 +92,9 @@ function SciMLBase.symbolic_discretize(pdesys::ModelingToolkit.PDESystem,discret
space = map(nottime) do x
xdomain = pdesys.domain[findfirst(d->isequal(x, d.variables),pdesys.domain)]
dx = discretization.dxs[findfirst(dxs->isequal(x, dxs[1].val),discretization.dxs)][2]

dx isa Number ? (DomainSets.infimum(xdomain.domain):dx:DomainSets.supremum(xdomain.domain)) : dx

end
dxs = map(nottime) do x
dx = discretization.dxs[findfirst(dxs->isequal(x, dxs[1].val),discretization.dxs)][2]
Expand Down Expand Up @@ -188,7 +191,9 @@ function SciMLBase.symbolic_discretize(pdesys::ModelingToolkit.PDESystem,discret
for bc in pdesys.bcs
bcdepvar = first(get_depvars(bc.lhs, depvar_ops))
if any(u->isequal(operation(u),operation(bcdepvar)),depvars)

if t != nothing && operation(bc.lhs) isa Sym && !any(x -> isequal(x, t.val), arguments(bc.lhs))

# initial condition
# Assume in the form `u(...) ~ ...` for now
i = findfirst(isequal(bc.lhs),initmaps)
Expand Down Expand Up @@ -258,7 +263,7 @@ function SciMLBase.symbolic_discretize(pdesys::ModelingToolkit.PDESystem,discret
Imax(order) = last(grid_indices) - I1 * (order÷2)

interior = grid_indices[[let bcs = get_bc_counts(i)
(1 + first(bcs)):length(g)-last(bcs)
(1 + first(bcs)):length(g)-last(bcs)#-1
end
for (i,g) in enumerate(grid)]...]
eqs = vec(map(interior) do II
Expand Down Expand Up @@ -291,16 +296,19 @@ function SciMLBase.symbolic_discretize(pdesys::ModelingToolkit.PDESystem,discret
central_deriv_rules_spherical = [Differential(s)(s^2*Differential(s)(u))/s^2 => central_deriv_spherical(II,j,k)
for (j,s) in enumerate(nottime), (k,u) in enumerate(depvars)]

# Val rules ############################################################
valrules = vcat([depvars[k] => depvarsdisc[k][II] for k in 1:length(depvars)],
[nottime[j] => grid[j][II[j]] for j in 1:nspace])

# TODO: upwind rules needs interpolation into `@rule`
forward_weights(II,j) = DiffEqOperators.calculate_weights(discretization.upwind_order, 0.0, grid[j][[II[j],II[j]+1]])

# Upwind rules #########################################################
reverse_weights(II,j) = DiffEqOperators.calculate_weights(discretization.upwind_order, 0.0, grid[j][[II[j]-1,II[j]]])
# upwinding_rules = [@rule(*(~~a,$(Differential(nottime[j]))(u),~~b) => IfElse.ifelse(*(~~a..., ~~b...,)>0,
# *(~~a..., ~~b..., dot(reverse_weights(II,j),depvars[k][central_neighbor_idxs(II,j)[1:2]])),
# *(~~a..., ~~b..., dot(forward_weights(II,j),depvars[k][central_neighbor_idxs(II,j)[2:3]]))))
# for j in 1:nspace, k in 1:length(pdesys.depvars)]
forward_weights(II,j) = DiffEqOperators.calculate_weights(discretization.upwind_order, 0.0, grid[j][[II[j],II[j]+1]])
#forward_weights(II,j) = -1.0 * reverse_weights(II,j) # TODO: check this function
side = 1.0
upwinding_rules_tmp = [@rule(*(~~a,$(Differential(iv))(dv),~~b) => Base.ifelse(*(side, ~~a..., ~~b...,)>0,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

to flip, I think you need @acrule @shashi ?

*(~~a..., ~~b..., dot(reverse_weights(II,j),depvarsdisc[k][central_neighbor_idxs(II,j,approx_order)[1:2]])),
*(~~a..., ~~b..., dot(forward_weights(II,j),depvarsdisc[k][central_neighbor_idxs(II,j,approx_order)[2:3]]))))
for (j, iv) in enumerate(nottime) for (k, dv) in enumerate(depvars)]

## Discretization of non-linear laplacian.
# 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
Expand All @@ -315,28 +323,52 @@ function SciMLBase.symbolic_discretize(pdesys::ModelingToolkit.PDESystem,discret
# Independent variable rules
r_mid_indep(II, j, l) = [nottime[j] => iv_mid(II, j, l) for j in 1:length(nottime)]
# Replacement rules: new approach
rules = [@rule ($(Differential(iv))(*(~~a, $(Differential(iv))(dv), ~~b))) =>
dot([Num(substitute(substitute(*(~~a..., ~~b...), r_mid_dep(II, j, k, -1)), r_mid_indep(II, j, -1))),
Num(substitute(substitute(*(~~a..., ~~b...), r_mid_dep(II, j, k, 1)), r_mid_indep(II, j, 1)))],
[-b1(II, j, k), b2(II, j, k)])
for (j, iv) in enumerate(nottime) for (k, dv) in enumerate(depvars)]
rhs_arg = istree(eq.rhs) && (SymbolicUtils.operation(eq.rhs) == +) ? SymbolicUtils.arguments(eq.rhs) : [eq.rhs]
lhs_arg = istree(eq.lhs) && (SymbolicUtils.operation(eq.lhs) == +) ? SymbolicUtils.arguments(eq.lhs) : [eq.lhs]
nonlinlap_rules_tmp = [@rule ($(Differential(iv))(*(~~a, $(Differential(iv))(dv), ~~b))) =>
dot([Num(substitute(substitute(*(~~a..., ~~b...), r_mid_dep(II, j, k, -1)), r_mid_indep(II, j, -1))),
Num(substitute(substitute(*(~~a..., ~~b...), r_mid_dep(II, j, k, 1)), r_mid_indep(II, j, 1)))],
[-b1(II, j, k), b2(II, j, k)])
for (j, iv) in enumerate(nottime) for (k, dv) in enumerate(depvars)]

# Post-processing @rules for applying `substitute` (see below) #########
lhs_arg = (SymbolicUtils.istree(eq.lhs) && SymbolicUtils.operation(eq.lhs) == +) ?
SymbolicUtils.arguments(eq.lhs) : [eq.lhs]
rhs_arg = (SymbolicUtils.istree(eq.rhs) && SymbolicUtils.operation(eq.rhs) == +) ?
SymbolicUtils.arguments(eq.rhs) : [eq.rhs]
nonlinlap_rules = []
lhs_upwinding_rules = []
rhs_upwinding_rules = []
for t in vcat(lhs_arg)
side = 1.0
for r in upwinding_rules_tmp
if r(t) != nothing
push!(lhs_upwinding_rules, t => r(t))
end
end
end
for t in vcat(rhs_arg)
side = -1.0
for r in upwinding_rules_tmp
if r(t) != nothing
push!(rhs_upwinding_rules, t => r(t))
end
end
end
for t in vcat(lhs_arg,rhs_arg)
for r in rules
if r(t) !== nothing
for r in nonlinlap_rules_tmp
if r(t) != nothing
push!(nonlinlap_rules, t => r(t))
end
end
end

# Applying rules to the equation #######################################
rules = vcat(vec(nonlinlap_rules),
vec(central_deriv_rules_cartesian),
vec(central_deriv_rules_spherical),
valrules)

substitute(eq.lhs,rules) ~ substitute(eq.rhs,rules)
vec(central_deriv_rules_cartesian),
vec(central_deriv_rules_spherical),
valrules)
lhs_tmp = substitute(eq.lhs,lhs_upwinding_rules)
rhs_tmp = substitute(eq.rhs,rhs_upwinding_rules)
substitute(lhs_tmp,rules) ~ substitute(rhs_tmp,rules)

end)
push!(alleqs,eqs)
Expand Down Expand Up @@ -374,3 +406,4 @@ function SciMLBase.discretize(pdesys::ModelingToolkit.PDESystem,discretization::
return prob = ODEProblem(simpsys,Pair[],tspan)
end
end

91 changes: 46 additions & 45 deletions test/MOL/MOL_1D_HigherOrder.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -74,49 +74,50 @@ end
prob = discretize(pdesys,discretization)
end

@testset "Kuramoto–Sivashinsky equation" begin
@parameters x, t
@variables u(..)
Dt = Differential(t)
Dx = Differential(x)
Dx2 = Differential(x)^2
Dx3 = Differential(x)^3
Dx4 = Differential(x)^4

α = 1
β = 4
γ = 1
eq = Dt(u(x,t)) ~ -u(x,t)*Dx(u(x,t)) - α*Dx2(u(x,t)) - β*Dx3(u(x,t)) - γ*Dx4(u(x,t))

u_analytic(x,t;z = -x/2+t) = 11 + 15*tanh(z) -15*tanh(z)^2 - 15*tanh(z)^3
du(x,t;z = -x/2+t) = 15/2*(tanh(z) + 1)*(3*tanh(z) - 1)*sech(z)^2

bcs = [u(x,0) ~ u_analytic(x,0),
u(-10,t) ~ u_analytic(-10,t),
u(10,t) ~ u_analytic(10,t),
Dx(u(-10,t)) ~ du(-10,t),
Dx(u(10,t)) ~ du(10,t)]
#@testset "Kuramoto–Sivashinsky equation" begin
# @parameters x, t
# @variables u(..)
# Dt = Differential(t)
# Dx = Differential(x)
# Dx2 = Differential(x)^2
# Dx3 = Differential(x)^3
# Dx4 = Differential(x)^4

# α = 1
# β = 4
# γ = 1
# eq = Dt(u(x,t)) ~ -u(x,t)*Dx(u(x,t)) - α*Dx2(u(x,t)) - β*Dx3(u(x,t)) - γ*Dx4(u(x,t))

# u_analytic(x,t;z = -x/2+t) = 11 + 15*tanh(z) -15*tanh(z)^2 - 15*tanh(z)^3
# du(x,t;z = -x/2+t) = 15/2*(tanh(z) + 1)*(3*tanh(z) - 1)*sech(z)^2

# bcs = [u(x,0) ~ u_analytic(x,0),
# u(-10,t) ~ u_analytic(-10,t),
# u(10,t) ~ u_analytic(10,t),
# Dx(u(-10,t)) ~ du(-10,t),
# Dx(u(10,t)) ~ du(10,t)]

# # Space and time domains
# domains = [x ∈ Interval(-10.0,10.0),
# t ∈ Interval(0.0,1.0)]
# # Discretization
# dx = 0.4; dt = 0.2

# discretization = MOLFiniteDifference([x=>dx],t;centered_order=4,grid_align=center_align)
# pdesys = PDESystem(eq,bcs,domains,[x,t],[u(x,t)])
# prob = discretize(pdesys,discretization)

# sol = solve(prob,Tsit5(),saveat=0.1,dt=dt)

# @test sol.retcode == :Success

# xs = domains[1].domain.lower+dx+dx:dx:domains[1].domain.upper-dx-dx
# ts = sol.t

# u_predict = sol.u
# u_real = [[u_analytic(x, t) for x in xs] for t in ts]
# u_diff = u_real - u_predict
# @test_broken u_diff[:] ≈ zeros(length(u_diff)) atol=0.01;
# #plot(xs, u_diff)
#end

# Space and time domains
domains = [x ∈ Interval(-10.0,10.0),
t ∈ Interval(0.0,1.0)]
# Discretization
dx = 0.4; dt = 0.2

discretization = MOLFiniteDifference([x=>dx],t;centered_order=4,grid_align=center_align)
pdesys = PDESystem(eq,bcs,domains,[x,t],[u(x,t)])
prob = discretize(pdesys,discretization)

sol = solve(prob,Tsit5(),saveat=0.1,dt=dt)

@test sol.retcode == :Success

xs = domains[1].domain.lower+dx+dx:dx:domains[1].domain.upper-dx-dx
ts = sol.t

u_predict = sol.u
u_real = [[u_analytic(x, t) for x in xs] for t in ts]
u_diff = u_real - u_predict
@test_broken u_diff[:] ≈ zeros(length(u_diff)) atol=0.01;
#plot(xs, u_diff)
end
Loading