the solution looks weird in https://docs.sciml.ai/MethodOfLines/stable/tutorials/heatss/
why its u(1, 1) = 0?
Yeah, that’s odd. @xtalax
I’m not sure if this is the intended behavior (still in the process of digging in the code), but it feels like the equations from symbolic_discretize show the wrong boundary conditions:
prob_sym = symbolic_discretize(pdesys, discretization)
for bc in prob_sym[1].eqs[end-3:end]
@show bc
end
# bc = 0 ~ -u[1, 1]
# bc = 0 ~ -u[11, 1]
# bc = 0 ~ -u[1, 11]
# bc = 0 ~ -u[11, 11]
The last one should be 1 ~ u[11,11] I guess. It also looks the same if one specifies it explicitly as
bcs
bcs = [u(0, y) ~ 0,
u(1, y) ~ y,
u(x, 0) ~ 0,
u(x, 1) ~ x]
Perhaps this is a starting point?
PS: The boundary conditions for the other 36 points look correct, it’s just the u[11,11] corner.
That would definitely be a bug. Please open an issue.
Is it related to this issue?
If I understand it correctly, the corner values are not that important anyways… for (continuous) Dirichlet b.c. their value is known and for Neumann/mixed b.c. their value needs to be consistent (in some sense) with the adjacent parts of the boundary, which would require some extrapolation.
Yes this is related. The corner equations have no effect on the interior solution, and their states do not appear in any other equation. As such they are currently held at 0, but I want to interpolate them using adjacent values so that they appear more in line with the rest of the solution.
1 Like
But this is only true for the “usual” stencils, right?
One could, in principle, come up with some stencils that use the corner points also for an equation in the interior. Not sure if this would be useful for anything though…
In principle yes, and in that case a choice will need to be made about which equation governs the corner states - I don’t know of any discretization where this is the case but if you have literature I’d love to see it.
I was mostly just curious about it, but apparently some people have been thinking about it at least theoretically
The Python package doesn’t really use this as far as I see, but they show an example of taking the derivatives in a 45º rotated basis, then the stencil looks like an “X”. But as I said, I’m not sure if there is any benefit to this and I haven’t seen it anywhere in the wild 
BTW, it seems that MethodOfLines.jl only works when the coefficients are constant, for example:
dT/dt = -alpha*d^2T/dz^2
with alpha being constant. It will not work if the coefficients are functions of the dependent variables, for example:
dT/dt = -alpha(T)*d^2T/dz^2
with alpha being dependent on T
BTW, it seems that
MethodOfLines.jlonly works when the coefficients are constant
No, this should work. Did you @register_symbolic alpha(T)?
Thanks!
Then should I define for example
alpha(T) = a*T+b
first, then
@register_symbolic alpha(T)
, and then define the equation?
Also, which definition is correct (or recommended) ?
eq = Dt(T(x,t)) ~ alpha(T(x,t))*Dxx(T(x,t))
or
eq = Dt(T) ~ alpha(T)*Dxx(T) ?
The first method is generally better, you want to define your T like @variables T(..) to ease defining your bcs and ics, then T on its own is just the operation