Impedance with multiple species

[sezerh24]

Hi @j-fu,

I have tested the ODE solver in the impedance examples and it works. thank you for fixing this. However, when I have a system that have multiple species, the transient solvers (ODEs and built in) do not work with the params (in example 151). I am trying to calculate the impedance of SOFCs (half cell) by using your embedded impedance method. However, in my system, I have two species for the electrical potential (i.e. ionic and electrical conductive phases). also, I have two gas species (i.e. H2 and H2O). So the transient solvers can not pass params = [1.0]. below, I had created a hypothetical case where there are 2 species. this is not working as well.

`using VoronoiFVM
using ExtendableGrids: geomspace, simplexgrid
using GridVisualize
using OrdinaryDiffEqSDIRK

using GLMakie

nref = 0
Plotter = GLMakie
verbose = false
unknown_storage = :sparse
assembly = :edgewise
time_embedding = :none

L = 1.0
R = 1.0
D = 1.0
C = 1.0
ω0 = 1.0e-3
ω1 = 5.0e1

Create array which is refined close to 0

h0 = 0.005 / 2.0^nref
h1 = 0.1 / 2.0^nref

X = geomspace(0, L, h0, h1)

Create discretitzation grid

grid = simplexgrid(X)

Create and fill data

data = (R = R, D = D, C = C)

Declare constitutive functions

flux = function (f, u, edge, data)
f[1] = data.D * (u[1, 1]*u[2, 1] - u[1, 2]*u[2, 2])
f[2] = data.D * (u[2, 1] - u[2, 2])
end

storage = function (f, u, node, data)
f[1] = data.C * u[1]
f[2] = data.C * u[2]
end

reaction = function (f, u, node, data)
f[1] = data.R * u[1]u[2] - u[2]
f[2] = data.R * u[2] u[1]
end

excited_bc= 1
excited_bcval = 1.0
excited_spec = 1
meas_bc = 2

bc = function (f, u, node, data)
p = parameters(u)
boundary_dirichlet!(f, u, node; species = 2, region = 1, value = 1.0)
boundary_dirichlet!(f, u, node; species = 1, region = excited_bc, value = p[1])
boundary_dirichlet!(f, u, node; species = 1, region = meas_bc, value = 0.0)
end

sys = VoronoiFVM.System(grid; unknown_storage = unknown_storage,
data = data,
flux = flux,
storage = storage,
reaction = reaction,
bcondition = bc,
nparams = 1,
assembly = assembly)

enable_species!(sys, 1, [1])
enable_species!(sys, 2, [1])
factory = TestFunctionFactory(sys)
measurement_testfunction = testfunction(factory, [excited_bc], [meas_bc])

tend=1.0

Transient solvers

#tsol=solve(sys; inival = 0.0, params = [1.0], times=(0.0,tend),force_first_step=true)

inival=unknowns(sys,inival=0)
problem = ODEProblem(sys,inival,(0,tend); params = [1.0])
odesol = solve(problem,ImplicitEuler())
tsol=reshape(odesol,sys)
steadystate=tsol.u[end]`