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

Commit a08a6cd

Browse files
ChrisRackauckasChrisRackauckas
authored
Merge pull request #250 from emmanuellujan/wip_pde_autodiscretization
[WIP] First steps in PDE autodiscretization.
2 parents fb9e717 + dc7b601 commit a08a6cd

6 files changed

Lines changed: 711 additions & 18 deletions

File tree

‎src/MOL_discretization.jl‎

Lines changed: 231 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -1,22 +1,236 @@
1+
# Method of lines discretization scheme
12
struct MOLFiniteDifference{T} <: DiffEqBase.AbstractDiscretization
2-
dxs::T
3-
order::Int
3+
dxs::T
4+
order::Int
5+
MOLFiniteDifference(args...;order=2) = new{typeof(args[1])}(args[1],order)
46
end
5-
MOLFiniteDifference(args...;order=2) = MOLFiniteDifference(args,order)
67

8+
# Get boundary conditions from an array
9+
function get_bcs(bcs,tdomain,domain)
10+
lhs_deriv_depvars_bcs = Dict()
11+
no_bcs = size(bcs,1)
12+
for i = 1:no_bcs
13+
var = bcs[i].lhs.op
14+
if var isa Variable
15+
var = var.name
16+
if !haskey(lhs_deriv_depvars_bcs,var)
17+
lhs_deriv_depvars_bcs[var] = Array{Expr}(undef,3)
18+
end
19+
j = 0
20+
if isequal(bcs[i].lhs.args[1],tdomain.lower) # u(t=0,x)
21+
j = 1
22+
elseif isequal(bcs[i].lhs.args[2],domain.lower) # u(t,x=x_init)
23+
j = 2
24+
elseif isequal(bcs[i].lhs.args[2],domain.upper) # u(t,x=x_final)
25+
j = 3
26+
end
27+
if bcs[i].rhs isa ModelingToolkit.Constant
28+
lhs_deriv_depvars_bcs[var][j] = :(var=$(bcs[i].rhs.value))
29+
else
30+
lhs_deriv_depvars_bcs[var][j] = Expr(bcs[i].rhs)
31+
end
32+
end
33+
end
34+
return lhs_deriv_depvars_bcs
35+
end
36+
37+
38+
# Recursively traverses the input expression (rhs), replacing derivatives by
39+
# finite difference schemes. It returns a time dependent expression (expr)
40+
# that will be evaluated in the "f" ODE function (in DiffEqBase.discretize),
41+
# Note: 'non-derived' dependent variables are inserted into the diff. equations
42+
# E.g. Dx(u(t,x))=v(t,x)*Dx(u(t,x)), v(t,x)=t*x
43+
# => Dx(u(t,x))=t*x*Dx(u(t,x))
44+
45+
function discretize_2(input,deriv_order,approx_order,dx,X,len,
46+
deriv_var,dep_var_idx,indep_var_idx)
47+
if input isa ModelingToolkit.Constant
48+
return :($(input.value))
49+
elseif input isa Operation
50+
if input.op isa Variable
51+
expr = :(0.0)
52+
var = input.op.name
53+
if haskey(indep_var_idx,var) # ind. var.
54+
if var != :(t)
55+
i = indep_var_idx[var]
56+
expr = :($X[$i][2:$len[$i]-1])
57+
else
58+
expr = :(t)
59+
end
60+
else # dep. var.
61+
# TODO: time and cross derivatives terms
62+
i = indep_var_idx[deriv_var[1]]
63+
j = dep_var_idx[var]
64+
if deriv_order == 0
65+
expr = :(u[:,$j])
66+
elseif deriv_order == 1
67+
# TODO: approx_order and forward/backward should be
68+
# input parameters of each derivative
69+
approx_order = 1
70+
L = UpwindDifference(deriv_order,approx_order,dx[i],len[i]-2,-1)
71+
expr = :(-1*($L*Q[$j]*u[:,$j]))
72+
elseif deriv_order == 2
73+
L = CenteredDifference(deriv_order,approx_order,dx[i],len[i]-2)
74+
expr = :($L*Q[$j]*u[:,$j])
75+
end
76+
end
77+
return expr
78+
elseif input.op isa Differential
79+
var = input.op.x.op.name
80+
push!(deriv_var,var)
81+
return discretize_2(input.args[1],deriv_order+1,approx_order,dx,X,
82+
len,deriv_var,dep_var_idx,indep_var_idx)
83+
pop!(deriv_var,var)
84+
else
85+
if size(input.args,1) == 1
86+
aux = discretize_2(input.args[1],deriv_order,approx_order,dx,X,
87+
len,deriv_var,dep_var_idx,indep_var_idx)
88+
return :(broadcast($(input.op), $aux))
89+
else
90+
aux_1 = discretize_2(input.args[1],deriv_order,approx_order,dx,X,
91+
len,deriv_var,dep_var_idx,indep_var_idx)
92+
aux_2 = discretize_2(input.args[2],deriv_order,approx_order,dx,X,
93+
len,deriv_var,dep_var_idx,indep_var_idx)
94+
return :(broadcast($(input.op), $aux_1, $aux_2))
95+
end
96+
end
97+
end
98+
end
99+
100+
# Convert a PDE problem into an ODE problem
7101
function DiffEqBase.discretize(pdesys::PDESystem,discretization::MOLFiniteDifference)
8-
tdomain = pdesys.domain[1].domain
9-
domain = pdesys.domain[2].domain
10-
@assert domain isa IntervalDomain
11-
len = domain.upper - domain.lower
12-
dx = discretization.dxs[1]
13-
interior = domain.lower+dx:dx:domain.upper-dx
14-
X = domain.lower:dx:domain.upper
15-
L = CenteredDifference(2,2,dx,Int(len/dx)-2)
16-
Q = DirichletBC(0.0,0.0)
17-
function f(du,u,p,t)
18-
mul!(du,L,Array(Q*u))
19-
end
20-
u0 = @. - interior * (interior - 1) * sin(interior)
21-
PDEProblem(ODEProblem(f,u0,(tdomain.lower,tdomain.upper),nothing),Q,X)
102+
103+
# TODO: discretize the following cases
104+
#
105+
# 1) PDE System
106+
# 1.a) Transient
107+
# There is more than one indep. variable, including 't'
108+
# E.g. du/dt = d2u/dx2 + d2u/dy2 + f(t,x,y)
109+
# 1.b) Stationary
110+
# There is more than one indep. variable, 't' is not included
111+
# E.g. 0 = d2u/dx2 + d2u/dy2 + f(x,y)
112+
# 2) ODE System
113+
# 't' is the only independent variable
114+
# The ODESystem is packed inside a PDESystem
115+
# E.g. du/dt = f(t)
116+
#
117+
# Note: regarding input format, lhs must be "du/dt" or "0".
118+
#
119+
120+
# The following code deals with 1.a case for 1D,
121+
# i.e. only considering 't' and 'x'
122+
123+
124+
### Declare and define independent-variable data structures ###############
125+
126+
tdomain = 0.0
127+
indep_var_idx = Dict()
128+
no_indep_vars = size(pdesys.domain,1)
129+
domain = Array{Any}(undef,no_indep_vars)
130+
dx = Array{Any}(undef,no_indep_vars)
131+
X = Array{Any}(undef,no_indep_vars)
132+
len = Array{Any}(undef,no_indep_vars)
133+
k = 0
134+
for i = 1:no_indep_vars
135+
var = pdesys.domain[i].variables.op.name
136+
indep_var_idx[var] = i
137+
domain[i] = pdesys.domain[i].domain
138+
if var != :(t)
139+
dx[i] = discretization.dxs[i-k]
140+
X[i] = domain[i].lower:dx[i]:domain[i].upper
141+
len[i] = size(X[i],1)
142+
else
143+
dx[i] = 0.0
144+
X[i] = 0.0
145+
len[i] = 0.0
146+
tdomain = pdesys.domain[1].domain
147+
@assert tdomain isa IntervalDomain
148+
k = 1
149+
end
150+
end
151+
152+
### Declare and define dependent-variable data structures #################
153+
154+
# TODO: specify order for each derivative
155+
approx_order = discretization.order
156+
157+
lhs_deriv_depvars = Dict()
158+
dep_var_idx = Dict()
159+
dep_var_disc = Dict() # expressions evaluated in the ODE function (f)
160+
161+
# if there is only one equation
162+
if pdesys.eq isa Equation
163+
eqs = [pdesys.eq]
164+
else
165+
eqs = pdesys.eq
166+
end
167+
no_dep_vars = size(eqs,1)
168+
for j = 1:no_dep_vars
169+
input = eqs[j].lhs
170+
if input.op isa Variable
171+
var = input.op.name
172+
else #var isa Differential
173+
var = input.args[1].op.name
174+
lhs_deriv_depvars[var] = j
175+
end
176+
dep_var_idx[var] = j
177+
end
178+
for (var,j) in dep_var_idx
179+
aux = discretize_2( eqs[j].rhs,0,approx_order,dx,X,len,
180+
[],dep_var_idx,indep_var_idx)
181+
# TODO: is there a better way to convert an Expr into a Function?
182+
dep_var_disc[var] = @eval (Q,u,t) -> $aux
183+
end
184+
185+
### Declare and define boundary conditions ################################
186+
187+
# TODO: extend to Neumann BCs and Robin BCs
188+
lhs_deriv_depvars_bcs = get_bcs(pdesys.bcs,tdomain,domain[2])
189+
t = 0.0
190+
u_t0 = Array{Float64}(undef,len[2]-2,no_dep_vars)
191+
u_x0 = Array{Any}(undef,no_dep_vars)
192+
u_x1 = Array{Any}(undef,no_dep_vars)
193+
Q = Array{RobinBC}(undef,no_dep_vars)
194+
195+
for var in keys(dep_var_idx)
196+
j = dep_var_idx[var]
197+
bcs = lhs_deriv_depvars_bcs[var]
198+
199+
g = eval(:((x,t) -> $(bcs[1])))
200+
u_t0[:,j] = @eval $g.($(X[2][2:len[2]-1]),$t)
201+
202+
u_x0[j] = @eval (x,t) -> $(bcs[2])
203+
u_x1[j] = @eval (x,t) -> $(bcs[3])
204+
205+
a = Base.invokelatest(u_x0[j],X[2][1],0.0)
206+
b = Base.invokelatest(u_x1[j],last(X[2]),0.0)
207+
Q[j] = DirichletBC(a,b)
208+
end
209+
210+
### Define the discretized PDE as an ODE function #########################
211+
212+
function f(du,u,p,t)
213+
214+
# Boundary conditions can vary with respect to time
215+
for j in 1:no_dep_vars
216+
a = Base.invokelatest(u_x0[j],X[2][1],t)
217+
b = Base.invokelatest(u_x1[j],last(X[2]),t)
218+
Q[j] = DirichletBC(a,b)
219+
end
220+
221+
for (var,disc) in dep_var_disc
222+
j = dep_var_idx[var]
223+
res = Base.invokelatest(disc,Q,u,t)
224+
if haskey(lhs_deriv_depvars,var)
225+
du[:,j] = res
226+
else
227+
u[:,j] .= res
228+
end
229+
end
230+
231+
end
232+
233+
# Return problem ##########################################################
234+
return PDEProblem(ODEProblem(f,u_t0,(tdomain.lower,tdomain.upper),nothing),Q,X)
22235
end
236+

‎test/MOL_0D_Logistic.jl‎

Lines changed: 53 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,53 @@
1+
# Logistic Equation (test incomplete)
2+
3+
# Packages and inclusions
4+
using ModelingToolkit,DiffEqOperators,DiffEqBase,LinearAlgebra,Test
5+
6+
# Tests
7+
@testset "Test 00: Dt(u(t)) ~ u(t)*(1.0-u(t))" begin
8+
# Parameters, variables, and derivatives
9+
@parameters t x
10+
@variables u(..)
11+
@derivatives Dt'~t
12+
@derivatives Dx'~x
13+
14+
# 1D PDE and boundary conditions
15+
eq = Dt(u(t,x)) ~ 0.0*Dx(u(t,x))+u(t,x)*(1.0-u(t,x))
16+
bcs = [ u(0,x) ~ 0.5+x*0.0,
17+
u(t,0) ~ 0.5,
18+
u(t,1) ~ 0.5]
19+
20+
# Space and time domains
21+
domains = [t ∈ IntervalDomain(0.0,5.0),
22+
x ∈ IntervalDomain(0.0,1.0)]
23+
24+
# PDE system
25+
pdesys = PDESystem(eq,bcs,domains,[t,x],[u])
26+
27+
# Method of lines discretization
28+
dx = 0.1
29+
order = 1
30+
discretization = MOLFiniteDifference(dx,order)
31+
32+
# Convert the PDE problem into an ODE problem
33+
prob = discretize(pdesys,discretization)
34+
35+
# Solve ODE problem
36+
using OrdinaryDiffEq
37+
38+
sol = solve(prob,Euler(),dt=0.01,saveat=0.1)
39+
40+
# Plot and save results
41+
using Plots
42+
time = domains[1].domain.lower:0.1:domains[1].domain.upper
43+
44+
plot(time,sol[5,1,:])
45+
savefig("MOL_0D_Logistic.png")
46+
47+
# Test
48+
# x_interval = domains[2].domain.lower+dx:dx:domains[2].domain.upper-dx
49+
# u = @. (0.5/(0.2*sqrt(2.0*3.1415)))*exp(-(x_interval-(0.75+0.6))^2/(2.0*0.2^2))
50+
# t_f = size(sol,3)
51+
# @test sol[t_f] ≈ u atol = 0.1;
52+
53+
end

0 commit comments

Comments
 (0)