SciML / MethodOfLines.jl

Automatic Finite Difference PDE solving with Julia SciML
https://docs.sciml.ai/MethodOfLines/stable/
MIT License
157 stars 27 forks source link

@parameters work unexpectedly #286

Open ctessum opened 1 year ago

ctessum commented 1 year ago

Hello,

The code below is a version of the Brusselator example, modified only by changing alpha from a numeric value to a @parameter:

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets

@parameters x y t
@variables u(..) v(..)
Dt = Differential(t)
Dx = Differential(x)
Dy = Differential(y)
Dxx = Differential(x)^2
Dyy = Differential(y)^2

∇²(u) = Dxx(u) + Dyy(u)

brusselator_f(x, y, t) = (((x-0.3)^2 + (y-0.6)^2) <= 0.1^2) * (t >= 1.1) * 5.

x_min = y_min = t_min = 0.0
x_max = y_max = 1.0
t_max = 11.5

@parameters α = 10.

u0(x,y,t) = 22(y*(1-y))^(3/2)
v0(x,y,t) = 27(x*(1-x))^(3/2)

eq = [Dt(u(x,y,t)) ~ 1. + v(x,y,t)*u(x,y,t)^2 - 4.4*u(x,y,t) + α*∇²(u(x,y,t)) + brusselator_f(x, y, t),
       Dt(v(x,y,t)) ~ 3.4*u(x,y,t) - v(x,y,t)*u(x,y,t)^2 + α*∇²(v(x,y,t))]

domains = [x ∈ Interval(x_min, x_max),
              y ∈ Interval(y_min, y_max),
              t ∈ Interval(t_min, t_max)]

# Periodic BCs
bcs = [u(x,y,0) ~ u0(x,y,0),
       u(0,y,t) ~ u(1,y,t),
       u(x,0,t) ~ u(x,1,t),

       v(x,y,0) ~ v0(x,y,0),
       v(0,y,t) ~ v(1,y,t),
       v(x,0,t) ~ v(x,1,t)] 

@named pdesys = PDESystem(eq,bcs,domains,[x,y,t],[u(x,y,t),v(x,y,t)],[α])
discretization = MOLFiniteDifference([x=>6, y=>6], t)
@time prob = discretize(pdesys,discretization)

This code gives the error: ERROR: type Num has no field first.

I can kind of fix this by specifying the default value of the parameter:

@named pdesys = PDESystem(eq,bcs,domains,[x,y,t],[u(x,y,t),v(x,y,t)],[α => 10.0])
@time prob = discretize(pdesys,discretization)

The code now runs, but it gives a warning: Warning: : no method matching get_unit for arguments (Pair{Num, Float64},). I guess this warning occurs because the parameter is now a pair rather than a single variable. This also goes against the PDESystem documentation, which specifies a separate argument for parameter defaults.

It seems like this might be a bug. Let me know if you have any thoughts!

ChrisRackauckas commented 1 year ago

@xtalax ?

xtalax commented 1 year ago

For some reason, PDESystems in the MOL test set have always been specified with param defaults supplied at this stage. I suppose this is because they were written before you could do this with @parameters. It is inconsistent with the rest of the ecosystem, but it has been this way for a while. I will add a path that uses default values from the definition of the variable like you suggest.

ctessum commented 1 year ago

Thanks!