`ExtraVariablesSystemException` when using a Robin BC

[SimonEnsemble]

hi! solving the diffusion equation with (i) concentration-dependent diffusion coefficient and (ii) Robin BC.

I get an error:

ExtraVariablesSystemException: The system is unbalanced. There are 102 highest order derivative variables and 101 equations.

More variables than equations, here are the potential extra variable(s):

here is our PDE:

	@parameters x t

	@variables c(..)

	∂t  = Differential(t)
	∂x  = Differential(x)
	∂²x = Differential(x) ^ 2

	diff_eq = ∂t(c(x, t)) ~ ∂x((D₀ + α * c(x, t)) * ∂x(c(x, t)))

	bcs = [c(x, 0) 	 ~ 0.0,  	       # initial condition
                  R * (D₀ + α * c(0, t)) * ∂x(c(0, t)) ~ c(0, t) - cₑ, # Robin
                  ∂x(c(ℓ, t)) ~ 0.0]   # no flux

so my student found that, if we re-arrange the Robin BC to instead look like:

	   R * ∂x(c(0, t)) + (cₑ - c(0, t)) / (D₀ + α * c(0, t)) ~ 0.0,

it works! maybe something to add to the docs?

(for novices like me, would be great to have c[1](t) and c[end](t) inferred by default instead of cutting off the solution at the boundaries.) thanks!
EDIT: nevermind, it doesn't start at c[3](t). I had an additional BC in there by accident. sorry. I can figure out how to infer c[1](t) from c[2](t). I'll edit this post when I figure it out. :)

[valentinsulzer]

Yeah the idea is that sol[c] will automatically reconstruct the full solution. I think for now something like

@variables c[1:n](t)
sol[c[1](t)]

should work (but having to redefine that first line is annoying)

[ChrisRackauckas]

@unpack c = sys
sol[c[1](t)]

if you did the symbolic discretize to get sys. We should make that a bit nicer.

[SimonEnsemble]

Valentin's way

@variables c[1:Nₓ](t)
sol[c[1](t)]

gives the error Sym c[1](t) is not callable. Use @syms c[1](t)(var1, var2,...) to create it as a callable.

changing to @syms c[1:Nₓ](t) gives MethodError: Cannot convert an object of type Expr to an object of type Symbol

Chris's way

@named pdesys = PDESystem(...)
@unpack c = pdesys # very much later

fails with type PDESystem has no field systems

what would be awesome

is if c[1](t) and c[Nₓ](t) appeared in the data frame! we are using the dataframe output.

my hack to account for the BCs and put them in the data frame output

(here, data is the solution to the PDE)

	#      Robin BC (write out first order finite difference)
	insertcols!(data, 2, "c[1](t)" => zeros(nrow(data)))
	for i = 2:nrow(data)
		c₂ = data[i, "c[2](t)"]
		c₁ = Polynomial([0, 1])
		poly = -α * c₁^2 + (-D₀ + α * c₂ - Δx / R) * c₁ + Δx / R * cₑ + D₀ * c₂
		poly_rts = roots(poly)
		data[i, "c[1](t)"] = poly_rts[isreal.(poly_rts)][1]
	end
	#    no flux BC
	insertcols!(data, "c[$(settings.Nₓ)](t)" => 
		        data[:, "c[$(settings.Nₓ-1)](t)"])

[valentinsulzer]

Mixed up the two ways, you can either do

@variables c[1:n](..)
sol[c[1](t)]

or

@variables c[1:n](t)
sol[c[1]]

[SimonEnsemble]

🙇‍♂️ thank you!

another weird thing is that for my diffusion eqn. with Robin BC and concentration-dependent diffusion, the "c1" IS automatically appended at the end! but not the end point corresponding to no-flux.

function toy_pde()
	@parameters x t

	@variables c(..)

	∂t  = Differential(t)
	∂x  = Differential(x)
	∂²x = Differential(x) ^ 2
	
	D₀ = 1.0
	α = 0.15
	R = 0.1
	cₑ = 2.0
	ℓ = 1.0
	Δx = 0.05

	diff_eq = ∂t(c(x, t)) ~ ∂x((D₀ + α * c(x, t)) * ∂x(c(x, t)))
	
	bcs = [
		# initial condition
		c(x, 0) ~ 0.0,
		# Robin BC
		∂x(c(0.0, t)) + (cₑ - c(0.0, t)) / (D₀ + α * c(0.0, t)) / R ~ 0.0, 
		# no flux BC
		∂x(c(ℓ, t)) ~ 0.0]
	

	# define space-time plane
	domains = [x ∈ Interval(0.0, ℓ), t ∈ Interval(0.0, 5.0)]

	# put it all together into a PDE system
	@named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]);

	# discretize space
	discretization = MOLFiniteDifference([x=>Δx], t)

	# convert the PDE into a system of ODEs via method of lines
	prob = discretize(pdesys, discretization)
	
	# compute number of spatial discretization points, boundaries included.
	Nₓ = length(prob.u0) + 1 # change to +2 for Dirichlet
	@assert (length(0:Δx:ℓ) == Nₓ)

	# solve system of ODEs in time.
	sol = solve(prob, saveat=0.01)
	
	# convert solution to data frame. solution at boundaries must be inferred.
	data = DataFrame(sol)
	@variables c[1:Nₓ](t)
	data[:, "c[$Nₓ](t)"] = sol[c[Nₓ]]
	
	ordered_cols = ["c[$i](t)" for i = 1:Nₓ]
	pushfirst!(ordered_cols, "timestamp")
	
	return data[:, ordered_cols]
end