Error with ODEProblem and sparse solution storage

[j-fu]

No description provided.

[JanSuntajs]

@j-fu, what is the status of this and what is it about?

I am asking, since I am bumping into this (or perhaps a similar) problem when trying to use ODEProblem with more than one species. For instance, if I try to adapt example 115 to use an ODEProblem, the error I get is the following:

AssertionError: length(v) == num_dof(system)

which is related to the way the initial condition is provided to the solver. Is there a way around this or one simply cannot use
external ODE solvers in this case?

Regards,
Jan

--

Edit: upon closer examination, I saw that the changing mass example in 1D works as expected, so I should probably be more specific in the description of my problem. I am trying to introduce an additional variable at the boundary. It should be related to the boundary flux $j_B$ through the relation $j_B = f(u_1, u_2)$, where $u_1$ is the first variable (for instance concentration in a diffusion problem); the value of $u_2$ should in general be determined from the relation just mentioned. The closest physical analog would be the galvanostatic battery operation, in which we prescribe some boundary flux (current) for the concentration, and the Butler-Volmer equation would determine the potential on the battery.

Would this be possible within the current VoronoiFVM framework?

[j-fu]

I consider the failing ODEProblem as a bug, will look into this.

A Butler-Volmer type boundary condition in the potentiostatic case shoul work. I need to think about galvanostatic. It might be possible that the package needs to be extended. Could you write up a simple example ?

[JanSuntajs]

@j-fu sure, this is a working example for the galvanostatic mode, which is quite easy, since it requires no additional modifications, just the driving flux should be voltage dependent. For the other case, I've had some ideas, such as defining an additional subgrid with only one element, but to no success so far. In Gridap, I have this implemented via a variable of a ConstantFESPace type which is something I'd probably need here as well. For now, the main problem remains the above mentioned assertion error - the number of degrees of freedom is correct, but length(v) equals the number of nodes, which would probably introduce additional problems in the assembly stage. Anyhow, here is the MWE in which a homogenous diffusion equation in spherical symmetry is driven in the potentiostatic mode; before solving the system, all the terms are cast in the dimensionless form; if plotting=true, an animation is drawn. This was tested with the latest (2.0.1) version of VoronoiFVM.

# imports
using Printf
using VoronoiFVM
using Plots
using DifferentialEquations
using ExtendableGrids
# using OrdinaryDiffEq
using Suppressor

# params - some physically relevant values
# physical parameters
const D0 = 1e-014;  # Diffusion coefficient [m^2/s]
const Cmax = 5.0 * 1e04; # Maximal cconcentration [mol/m^3]
const C0 = 0.5 # initial concentration in dimensionless units 
const R0 = 1e-05; # Particle radius [m]
const θ0 = 1.0 # stoichiometry at SOC=0
const θ1 = 0.0 # stoichiometry at SOC=1
# flux parameters
const C_rate = 0.001 #0.001 # C rate
const j0 = D0 * Cmax / R0 # characteristic flux
const τ = R0^2 / D0 # characteristic time
# we determine the applied flux from the condition that a sphere
# should be fully charged in an hour (3600 seconds).
const j_app = R0 * Cmax * (θ0 - θ1) / (3 * 3600) # applied flux
# driving parameters
const freq = 0.001 #.1 # frequency of the driving [hz]
# frequency in the dimensionless units [how many oscillations during ones
# characteristic time]
const freq_dimless = freq * τ # dimensionless frequency
const n_osc = 20 # number of oscillations during the simulation
const tf = n_osc / freq # final time in physical units
const RT = 8.3144621 * 300. # gas constant times temperature
const F = 96485.3365 # Faraday constant
const i0 = 2. # exchange current density
const α = 0.5 
const UOCV = 4.0 # open circuit voltage
const U1 = 0.001 # amplitude of the driving voltage
const C0 = 0.5 # initial concentration
const plotting = true
# note: we cast everything in dimensionless units to aid the performance
# of the solver

# boundary flux function divided by the characteristic flux to
# make it dimensionless
_voltage(t; U0=UOCV - C0, U1=-U1) = U0 + U1 * sin(2π * freq_dimless * t)

"""
Butler-Volmer current which also depends on the surface concentration
in this case.
"""
_flux_bv(t, c; U0=UOCV - C0, U1=-U1 ) = (i0 / F) * (exp(-α * F * (_voltage(t; U0=U0, U1=U1) + c - UOCV) / RT) - exp(α * F * (_voltage(t; U0=U0, U1=U1) + c - UOCV) / RT)) / j0;

# build params for the assembly
const unknown_storage::Symbol = :sparse
const assembly::Symbol = :edgewise
const geometry::Symbol = :spherical_symmetry;

# Solver params
const solver_params = Dict(
    :ti => 0.0,
    :tf => tf / τ,
)



const ode_params = Dict(
    :abstol => 1e-12,
    :reltol => 1e-12,
    :saveat => collect(solver_params[:ti]:0.01:solver_params[:tf]),
)

# boundary condition
#
# We drive the system using a periodically changing voltage
# that drives the Butler-Volmer current
#
function bcond_driving!(f, u, bnode, data)

	# _bflux(t) = sin(2π * t)
	boundary_neumann!(f, u, bnode, species = 1, region = 2, value = _flux_bv(bnode.time, u[1]))
	boundary_neumann!(f, u, bnode, species = 1, region = 1, value = 0.0)

end

"""
The flux term for the PDE
"""
function flux!(f, u, edge, data)
    _D = data.D
    f[1] = _D * (u[1, 1] - u[1, 2])
end

function storage!(f, u, node, data)
    f[1] = u[1]
end
# define the grid
# Define the grid

const n = 200
h = 1.0 / convert(Float64, n)
X = collect(0:h:1)

_grid = simplexgrid(X)
spherical_symmetric!(_grid)

# we set D equal to 1. assuming the equation is dimensionless
physics = VoronoiFVM.Physics(;flux=flux!, storage=storage!, data = (D=1., ), breaction = bcond_driving!);
sys = VoronoiFVM.System(_grid, physics; 
unknown_storage=unknown_storage, assembly=assembly)

# add species 1, run the simulation
enable_species!(sys, 1, [1])

# prepare the initial condition
inival = unknowns(sys)
# inival .= 0.
inival[1, :] .= map(x -> 0.5, X)

# solve the system: use the Forward-Backward Differentiation Formula
# along with the DifferentialEquations.jl tools.
problem = ODEProblem(sys, inival,
    [solver_params[:ti], solver_params[:tf]])
odesol = solve(problem, FBDF(); ode_params...)
tsol_lin = reshape(odesol, sys);

if plotting
    default(titlefont=font("Computer Modern", 12), guidefont=font("Computer Modern", 11), 
    tickfont=font("Computer Modern", 10), legendfont=font("Computer Modern", 10))
    
    
    surf_conc = [tsol_lin.u[i][end] for i in eachindex(tsol_lin.t)]
    
    @suppress begin
        anim = @animate for i in eachindex(tsol_lin.t[1:200])
    
            # Define a 1x2 layout (two panels side by side)
            l = @layout [a{0.9h} b{0.9h} c{0.9h}]
    
            # First panel: linear diffusion plot
            p1 = plot(X, tsol_lin.u[i][:], title="Concentration profile",
                label="linear diffusion", lw=2, color=:blue, linestyle=:dash, ylim=(0.499, .501))
            xlabel!(p1, "r")
            ylabel!(p1, "c")
            plot!(p1, legend=:topright)
    
            # Second panel: plot other data (replace `other_data[i]` with your desired data)
            # This could be, for example, another solution, a reference line, etc.
            p2 = hline([_voltage(tsol_lin.t[i])] , title="Voltage",
                label="voltage", lw=2, color=:red, ylim=(3.499, 3.501))  # Customize ylim as needed
            xlabel!(p2, "r")
            ylabel!(p2, "\$\\Delta \\phi\$")
            plot!(p2, legend=:topright)
    
            # p3 = plot(tsol_lin.t, _voltage.(tsol_lin.t),
            #     label="voltage", lw=2, color=:red,)  # Customize ylim as needed
            p3 = hline([_flux_bv(tsol_lin.t[i], surf_conc[i])], color=:blue, lw=2, label="flux",
            ylim = (-0.008, 0.006))
            plot!(p3, legend=:bottomright)
            # Combine both panels in a single plot
            plot(p1, p2, p3,  layout=l, size=(1200, 400))  # Adjust size as needed
    
        end
    
        # Display the animation
        gif(anim, "galvanostatic_driving_example.gif", fps=5)  # Save it as a gif

end

end

[j-fu]

Let us continue the discussion in #131