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

Commit dc7b601

Browse files
emmanuellujanemmanuellujan
committed
first steps in multidimensions
1 parent beddc6c commit dc7b601

3 files changed

Lines changed: 135 additions & 113 deletions

File tree

‎src/MOL_discretization.jl‎

Lines changed: 92 additions & 68 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,3 @@
1-
#using Reduce
2-
31
# Method of lines discretization scheme
42
struct MOLFiniteDifference{T} <: DiffEqBase.AbstractDiscretization
53
dxs::T
@@ -10,10 +8,11 @@ end
108
# Get boundary conditions from an array
119
function get_bcs(bcs,tdomain,domain)
1210
lhs_deriv_depvars_bcs = Dict()
13-
n = size(bcs,1)
14-
for i = 1:n
11+
no_bcs = size(bcs,1)
12+
for i = 1:no_bcs
1513
var = bcs[i].lhs.op
1614
if var isa Variable
15+
var = var.name
1716
if !haskey(lhs_deriv_depvars_bcs,var)
1817
lhs_deriv_depvars_bcs[var] = Array{Expr}(undef,3)
1918
end
@@ -43,37 +42,55 @@ end
4342
# E.g. Dx(u(t,x))=v(t,x)*Dx(u(t,x)), v(t,x)=t*x
4443
# => Dx(u(t,x))=t*x*Dx(u(t,x))
4544

46-
function discretize_2(input,deriv_order,approx_order,dx,X,len,index)
45+
function discretize_2(input,deriv_order,approx_order,dx,X,len,
46+
deriv_var,dep_var_idx,indep_var_idx)
4747
if input isa ModelingToolkit.Constant
4848
return :($(input.value))
4949
elseif input isa Operation
5050
if input.op isa Variable
5151
expr = :(0.0)
52-
if !haskey(index,input.op) # ind. var.
53-
expr = :($X)
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
5460
else # dep. var.
55-
j = index[input.op]
61+
# TODO: time and cross derivatives terms
62+
i = indep_var_idx[deriv_var[1]]
63+
j = dep_var_idx[var]
5664
if deriv_order == 0
5765
expr = :(u[:,$j])
5866
elseif deriv_order == 1
67+
# TODO: approx_order and forward/backward should be
68+
# input parameters of each derivative
5969
approx_order = 1
60-
L = UpwindDifference(deriv_order,approx_order,dx[1],len,-1)
70+
L = UpwindDifference(deriv_order,approx_order,dx[i],len[i]-2,-1)
6171
expr = :(-1*($L*Q[$j]*u[:,$j]))
6272
elseif deriv_order == 2
63-
L = CenteredDifference(deriv_order,approx_order,dx[1],len)
73+
L = CenteredDifference(deriv_order,approx_order,dx[i],len[i]-2)
6474
expr = :($L*Q[$j]*u[:,$j])
6575
end
6676
end
6777
return expr
6878
elseif input.op isa Differential
69-
return discretize_2(input.args[1],deriv_order+1,approx_order,dx,X,len,index)
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)
7084
else
7185
if size(input.args,1) == 1
72-
aux = discretize_2(input.args[1],deriv_order,approx_order,dx,X,len,index)
86+
aux = discretize_2(input.args[1],deriv_order,approx_order,dx,X,
87+
len,deriv_var,dep_var_idx,indep_var_idx)
7388
return :(broadcast($(input.op), $aux))
7489
else
75-
aux_1 = discretize_2(input.args[1],deriv_order,approx_order,dx,X,len,index)
76-
aux_2 = discretize_2(input.args[2],deriv_order,approx_order,dx,X,len,index)
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)
7794
return :(broadcast($(input.op), $aux_1, $aux_2))
7895
end
7996
end
@@ -103,99 +120,106 @@ function DiffEqBase.discretize(pdesys::PDESystem,discretization::MOLFiniteDiffer
103120
# The following code deals with 1.a case for 1D,
104121
# i.e. only considering 't' and 'x'
105122

106-
### Get domains (typically temporal and spatial) ###########################
107-
# TODO: here it is assumed that the time domain is the first in the array.
108-
# It can be in any order.
109-
110-
tdomain = pdesys.domain[1].domain
111-
@assert tdomain isa IntervalDomain
112-
113-
no_iv = size(pdesys.domain,1)
114-
domain = []
115-
dx = []
116-
X = []
117-
xx = []
118-
for i = 1:no_iv-1
119-
domain = vcat(domain,pdesys.domain[i+1].domain)
120-
dx = vcat(dx,discretization.dxs)
121-
X = vcat(X,domain[i].lower:dx[i]:domain[i].upper)
122-
xx = vcat(xx,size(X,1)-2)
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
123150
end
124-
interior = domain[1].lower+dx[1]:dx[1]:domain[1].upper-dx[1]
151+
152+
### Declare and define dependent-variable data structures #################
125153

126154
# TODO: specify order for each derivative
127155
approx_order = discretization.order
128156

129-
### Calculate discretization expression ####################################
130-
# The discretization is an expression which is then evaluated
131-
# in the ODE function (f)
132-
133-
discretization = Dict()
134157
lhs_deriv_depvars = Dict()
135-
index = Dict()
136-
158+
dep_var_idx = Dict()
159+
dep_var_disc = Dict() # expressions evaluated in the ODE function (f)
160+
137161
# if there is only one equation
138162
if pdesys.eq isa Equation
139163
eqs = [pdesys.eq]
140164
else
141165
eqs = pdesys.eq
142166
end
143-
144-
n_eqs = size(eqs,1)
145-
for j = 1:n_eqs
167+
no_dep_vars = size(eqs,1)
168+
for j = 1:no_dep_vars
146169
input = eqs[j].lhs
147170
if input.op isa Variable
148-
var = input.op
171+
var = input.op.name
149172
else #var isa Differential
150-
var = input.args[1].op
173+
var = input.args[1].op.name
151174
lhs_deriv_depvars[var] = j
152175
end
153-
index[var] = j
176+
dep_var_idx[var] = j
154177
end
155-
for (var,j) in index
156-
aux = discretize_2( eqs[j].rhs,0,approx_order,
157-
dx[1],interior,xx[1],index)
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)
158181
# TODO: is there a better way to convert an Expr into a Function?
159-
discretization[var] = @eval (Q,u,t) -> $aux
182+
dep_var_disc[var] = @eval (Q,u,t) -> $aux
160183
end
161184

162-
### Get boundary conditions ################################################
185+
### Declare and define boundary conditions ################################
186+
163187
# TODO: extend to Neumann BCs and Robin BCs
164-
lhs_deriv_depvars_bcs = get_bcs(pdesys.bcs,tdomain,domain[1])
188+
lhs_deriv_depvars_bcs = get_bcs(pdesys.bcs,tdomain,domain[2])
165189
t = 0.0
166-
u_t0 = Array{Float64}(undef,length(interior),length(discretization))
167-
u_x0 = Array{Any}(undef,length(discretization))
168-
u_x1 = Array{Any}(undef,length(discretization))
169-
Q = Array{RobinBC}(undef,length(discretization))
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)
170194

171-
for var in keys(discretization)
172-
j = index[var]
195+
for var in keys(dep_var_idx)
196+
j = dep_var_idx[var]
173197
bcs = lhs_deriv_depvars_bcs[var]
174198

175199
g = eval(:((x,t) -> $(bcs[1])))
176-
u_t0[:,j] = @eval $g.($interior,$t)
200+
u_t0[:,j] = @eval $g.($(X[2][2:len[2]-1]),$t)
177201

178202
u_x0[j] = @eval (x,t) -> $(bcs[2])
179203
u_x1[j] = @eval (x,t) -> $(bcs[3])
180204

181-
a = Base.invokelatest(u_x0[j],X[1],0.0)
182-
b = Base.invokelatest(u_x1[j],last(X),0.0)
205+
a = Base.invokelatest(u_x0[j],X[2][1],0.0)
206+
b = Base.invokelatest(u_x1[j],last(X[2]),0.0)
183207
Q[j] = DirichletBC(a,b)
184-
185208
end
186209

187-
### Define the discretized PDE as an ODE function ##########################
210+
### Define the discretized PDE as an ODE function #########################
188211

189212
function f(du,u,p,t)
190213

191-
for j in 1:length(discretization)
192-
a = Base.invokelatest(u_x0[j],X[1],t)
193-
b = Base.invokelatest(u_x1[j],last(X),t)
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)
194218
Q[j] = DirichletBC(a,b)
195219
end
196220

197-
for (var,disc) in discretization
198-
j = index[var]
221+
for (var,disc) in dep_var_disc
222+
j = dep_var_idx[var]
199223
res = Base.invokelatest(disc,Q,u,t)
200224
if haskey(lhs_deriv_depvars,var)
201225
du[:,j] = res
@@ -206,7 +230,7 @@ function DiffEqBase.discretize(pdesys::PDESystem,discretization::MOLFiniteDiffer
206230

207231
end
208232

209-
# Return problem ###########################################################
233+
# Return problem ##########################################################
210234
return PDEProblem(ODEProblem(f,u_t0,(tdomain.lower,tdomain.upper),nothing),Q,X)
211235
end
212236

‎test/MOL_1D_Linear_Convection.jl‎

Lines changed: 31 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -39,12 +39,12 @@ using ModelingToolkit,DiffEqOperators,DiffEqBase,LinearAlgebra,Test
3939

4040
# Plot and save results
4141
# using Plots
42-
# plot(prob.space,Array(prob.extrapolation[1]*sol[:,1,1]))
43-
# plot!(prob.space,Array(prob.extrapolation[1]*sol[:,1,2]))
44-
# plot!(prob.space,Array(prob.extrapolation[1]*sol[:,1,3]))
45-
# plot!(prob.space,Array(prob.extrapolation[1]*sol[:,1,4]))
46-
# plot!(prob.space,Array(prob.extrapolation[1]*sol[:,1,5]))
47-
# plot!(prob.space,Array(prob.extrapolation[1]*sol[:,1,6]))
42+
# plot(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,1]))
43+
# plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,2]))
44+
# plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,3]))
45+
# plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,4]))
46+
# plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,5]))
47+
# plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,6]))
4848
# savefig("MOL_1D_Linear_Convection_Test00.png")
4949

5050
# Test
@@ -92,12 +92,12 @@ end
9292

9393
# Plot and save results
9494
# using Plots
95-
# plot(prob.space,Array(prob.extrapolation[1]*sol[:,1,1]))
96-
# plot!(prob.space,Array(prob.extrapolation[1]*sol[:,1,2]))
97-
# plot!(prob.space,Array(prob.extrapolation[1]*sol[:,1,3]))
98-
# plot!(prob.space,Array(prob.extrapolation[1]*sol[:,1,4]))
99-
# plot!(prob.space,Array(prob.extrapolation[1]*sol[:,1,5]))
100-
# plot!(prob.space,Array(prob.extrapolation[1]*sol[:,1,6]))
95+
# plot(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,1]))
96+
# plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,2]))
97+
# plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,3]))
98+
# plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,4]))
99+
# plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,5]))
100+
# plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,6]))
101101
# savefig("MOL_1D_Linear_Convection_Test01.png")
102102

103103
# Test
@@ -144,12 +144,12 @@ end
144144

145145
# Plot and save results
146146
# using Plots
147-
# plot(prob.space,Array(prob.extrapolation[1]*sol[:,1,1]))
148-
# plot!(prob.space,Array(prob.extrapolation[1]*sol[:,1,2]))
149-
# plot!(prob.space,Array(prob.extrapolation[1]*sol[:,1,3]))
150-
# plot!(prob.space,Array(prob.extrapolation[1]*sol[:,1,4]))
151-
# plot!(prob.space,Array(prob.extrapolation[1]*sol[:,1,5]))
152-
# plot!(prob.space,Array(prob.extrapolation[1]*sol[:,1,6]))
147+
# plot(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,1]))
148+
# plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,2]))
149+
# plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,3]))
150+
# plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,4]))
151+
# plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,5]))
152+
# plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,6]))
153153
# savefig("MOL_1D_Linear_Convection_Test02.png")
154154

155155
# Test
@@ -198,12 +198,12 @@ end
198198

199199
#Plot and save results
200200
# using Plots
201-
# plot(prob.space,Array(prob.extrapolation[1]*sol[:,1,1]))
202-
# plot!(prob.space,Array(prob.extrapolation[1]*sol[:,1,2]))
203-
# plot!(prob.space,Array(prob.extrapolation[1]*sol[:,1,3]))
204-
# plot!(prob.space,Array(prob.extrapolation[1]*sol[:,1,4]))
205-
# plot!(prob.space,Array(prob.extrapolation[1]*sol[:,1,5]))
206-
# plot!(prob.space,Array(prob.extrapolation[1]*sol[:,1,6]))
201+
# plot(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,1]))
202+
# plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,2]))
203+
# plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,3]))
204+
# plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,4]))
205+
# plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,5]))
206+
# plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,6]))
207207
# savefig("MOL_1D_Linear_Convection_Test03.png")
208208

209209
# Test
@@ -250,14 +250,14 @@ end
250250
using OrdinaryDiffEq
251251
sol = solve(prob,Euler(),dt=.025,saveat=0.1)
252252

253-
#Plot and save results
253+
# Plot and save results
254254
# using Plots
255-
# plot(prob.space,Array(prob.extrapolation[1]*sol[:,1,1]))
256-
# plot!(prob.space,Array(prob.extrapolation[1]*sol[:,1,2]))
257-
# plot!(prob.space,Array(prob.extrapolation[1]*sol[:,1,3]))
258-
# plot!(prob.space,Array(prob.extrapolation[1]*sol[:,1,4]))
259-
# plot!(prob.space,Array(prob.extrapolation[1]*sol[:,1,5]))
260-
# plot!(prob.space,Array(prob.extrapolation[1]*sol[:,1,6]))
255+
# plot(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,1]))
256+
# plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,2]))
257+
# plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,3]))
258+
# plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,4]))
259+
# plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,5]))
260+
# plot!(prob.space[2],Array(prob.extrapolation[1]*sol[:,1,6]))
261261
# savefig("MOL_1D_Linear_Convection_Test04.png")
262262

263263
# Test

0 commit comments

Comments
 (0)