Visualization of the solution of a multifield time dependant problem

[Mahdimasmoudi0]

hi,
I am solving an electrophysiology problem using Gridap.
I have solved the problem, the issue is that while trying to write my solution in vtu file, I am getting this : ERROR: This function belongs to an interface definition and cannot be used.
The issue is within these lines :
createpvd("Cardiac_Electrophysiology_transient_solution") do pvd for (uₕ,t) in uₕₜ pvd[t] = createvtk(Ω,"Cardiac_Electrophysiology_transient_solution_$t"*".vtu",cellfields=["uh"=>uₕ[1],"ph"=>uₕ[2]]) end end

This is the complete code :
`using Gridap

mesh generation

n = 100
domain = (0,1,0,1)
partition = (n,n)
model = CartesianDiscreteModel(domain,partition)

Define test and trial spaces

order = 1
reffeᵤ = ReferenceFE(lagrangian,Float64,order)
V = TestFESpace(model,reffeᵤ,conformity=:H1)
reffeₚ = ReferenceFE(lagrangian,Float64,order)
Q = TestFESpace(model,reffeₚ,conformity=:L2)
U = TransientTrialFESpace(V)
P = TransientTrialFESpace(Q)
Y = MultiFieldFESpace([V, Q])
X = TransientMultiFieldFESpace([U, P])

Triangulation and integration quadrature

degree = order
Ωₕ = Triangulation(model)
dΩ = Measure(Ωₕ,degree)

weak formulation

α=0.01
c1= 8
c2=8
γ=0.002
μ1=0.2
μ2=0.3
bcnst=0.15
D = 0.1
f(ϕ,r)=c1ϕ(ϕ−α)(1-ϕ)-c2rϕ
g(ϕ,r)=( γ + μ1r*(1/(μ2+ϕ)))(-r-c2ϕ*(ϕ-bcnst-1))
res(t,(ϕ,r),(ν,η)) = ∫(∂t(ϕ)ν +∂t(r)η+ D( ∇(ν)⋅(∇(ϕ))) - νf(ϕ,r) - η*g(ϕ,r))*dΩ
op_AD = TransientFEOperator(res,X,Y)

Transient solver

linear_solver = LUSolver()
Δt = 0.05
θ = 0.5
ode_solver = ThetaMethod(linear_solver,Δt,θ)
t₀ = 0.0
T = 10.0
Uh0 = interpolate_everywhere(0.0, U(t₀)) # not like tutorial
Ph0 = interpolate_everywhere(0.0, P(t₀)) # not like tutorial
xh0 = interpolate_everywhere([Uh0, Ph0], X(t₀))
xhs0 = (xh0,)
uₕₜ = solve(ode_solver,op_AD,xhs0,t₀,T)
createpvd("Cardiac_Electrophysiology_transient_solution") do pvd
for (uₕ,t) in uₕₜ
pvd[t] = createvtk(Ω,"Cardiac_Electrophysiology_transient_solution_$t"*".vtu",cellfields=["uh"=>uₕ[1],"ph"=>uₕ[2]])
end
end`