MethodOfLines.jl icon indicating copy to clipboard operation
MethodOfLines.jl copied to clipboard

ExtraVariablesSystemException when u(t) depends on v(r,t) at specific r values

Open dtdowd opened this issue 2 years ago • 2 comments
trafficstars

Hello, I'm trying to solve a system where I have one equation that depends only on time (u(t)) and another that depends on time and one spatial dimension (v(r,t)). u(t) depends on v(r,t) at specific r values r1=1.0, r2=2.0, and r3=3.0.

The first thing I tried was calling v(r1,t), v(r2,t), v(r3,t) in the definition of Dt(u(t)). When I go to discretize the system, I get an error that "The system is unbalanced. There are 39 highest order derivative variables and 36 equations." However, when I scroll up in the REPL to view the output of discretize, right below the printout of the system of equations, it says "There are 36 variables and 36 equations."

I believe I have written the program so that r1, r2, and r3 are among the discretized r values (domains = [r ∈ (0.6, 4.0)...] and MOLFiniteDifference([r=>0.1], t)), so I thought the v(r1,t) call would be replaced by v(t)[index] at whatever the corresponding index for r1 is.

Is there a way to call v(r,t) at a given r value in the definition of Dt(u(t))? If so, what's the best way to do it?

I've tried a few solutions, but none of them seemed to be on the right track, so I've just posted my original code below. Any help is much appreciated!

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets

@parameters r r1 r2 r3 t;
@variables u(..) v(..);
Dt = Differential(t); Dr = Differential(r); Drr = Differential(r)^2;
r0, rmax = 0.6, 4.0; r1 = 1.0; r2 = 2.0; r3 = 3.0;
Q1 = 1.0; Q2 = 2.0; Q3 = 3.0; Q = Q1 + Q2 + Q3;

eq = [
	Dt(v(r,t)) ~ (1/r^2)*((r^2)*Drr(v(r,t)) + 2*r*Dr(v(r,t))) - (r^(-2))*Dr(v(r,t)),
    Dt(u(t)) ~ sum([Q1*v(r1,t), Q2*v(r2,t), Q3*v(r3,t)]) - u(t)]

domains = [r ∈ Interval(r0, rmax), t ∈ Interval(0.0, 5.0)]

bcs = [
    u(0) ~ 20.0,

    v(r,0) ~ 23.5-5.85*r,
    v(r0, t) ~ u(t),
    v(rmax, t) ~ 0.1
]

@named pdesys = PDESystem(eq, bcs, domains, [r,t], [u(t), v(r,t)])

discretization = MOLFiniteDifference([r => 0.1], t)

prob = discretize(pdesys, discretization);

dtdowd avatar Sep 13 '23 16:09 dtdowd

You are getting a common MTK Error, but I'm not certain you can solve this system with MOL. You can use values at a boundary in equations, but unfortunately not in the middle of the domain. You might be able to do this by connecting some separate variables defined on different domains together at boundaries to represent v when joined together but I have not tested this with multiple interfaces.

See the interfaces section here for how to do this: https://docs.sciml.ai/MethodOfLines/stable/boundary_conditions/

xtalax avatar Sep 13 '23 17:09 xtalax

Thanks for the tip

dtdowd avatar Sep 13 '23 20:09 dtdowd