TransientFESolution modifies the cellfields previously yielded

[Antoinemarteau]

When iterating on timesteps solutions of TransientFESolution, new iteration modifies the previously yielded solution.

Reproduce by replacing poisson_transient.jl tutorial postprocessing by :

# ## Postprocessing
(uh1, t1), state1  = iterate(uₕₜ)
writevtk(Ω, "uh1",    cellfields=[ "uh1"=>uh1])    # exports uh1 as expected

(uh2, t2), state2  = iterate(uₕₜ, state1)
writevtk(Ω, "uh1_KO", cellfields=[ "uh1_KO"=>uh1]) # exports something diff from uh1
writevtk(Ω, "uh2",    cellfields=[ "uh2"=>uh2])    # exports uh2 as expected

The second export of 1st time step solution uh1 is different from the fist one. This prevents me from visualizing $\partial_t(uh)$ as I never have two consecutive solutions available simultaneously.

[Antoinemarteau]

I see there is an allocate_trial_space in iterate(sol) but not in iterate(sol, state) https://github.com/gridap/Gridap.jl/blob/f4236fb4b89f61f20570973d8105dc85e7b1b1d4/src/ODEs/TransientFETools/TransientFESolutions.jl#L73

Would adding the allocation in the second iterate method solve it? Would it be better to implement copy(cellfield::CellDatum) to allow the user to save a yielded solution (and avoid allocation when user don't use it) ?

[fverdugo]

@Antoinemarteau I think this is the expected behaviour. It gets modified because we allocate a vector of free dof values and reuse it between time steps. Otherwise, allocating a new vector for each time step would be overkill in a large computation.

If you can allocate a new solution in your script and use it to store the previous step.

[oriolcg]

Hi @Antoinemarteau , you might be able to get the previous solution from state: https://github.com/gridap/Gridap.jl/blob/f4236fb4b89f61f20570973d8105dc85e7b1b1d4/src/ODEs/ODETools/ODESolutions.jl#L91

[Antoinemarteau]

@fverdugo I agree the iterator shouldn't allocate a dof buffer each iteration, but shouldn't it initially allocate N buffer for each time step used by the timesteping method ? (so N=2 for theta method)

I managed to achieve the copy with this code :

# ## Postprocessing
Up = Gridap.ODEs.TransientFETools.allocate_trial_space(uₕₜ.trial)

(uh1, t1), state1  = iterate(uₕₜ)
uhp = FEFunction(Up, copy(uh1.free_values))

(uh2, t2), state2  = iterate(uₕₜ, state1)
writevtk(Ω, "uhp1", cellfields=[ "uhp1"=>uhp]) # exports uh1 as expected
copyto!(uhp.free_values, uh2.free_values)

iterate(uₕₜ, state2)
writevtk(Ω, "uhp2", cellfields=[ "uhp2"=>uhp]) # exports uh2 as expected

Is this the best way to do this ? It would be very nice to have uhp=copy(uh) and an in-place copyto!(uhp, uh) for this use :P

[Antoinemarteau]

@oriolcg that would be perfect ! But I strugle understanding how, I got this

# ## Postprocessing
iter_result = iterate(uₕₜ)
while iter_result !== nothing
    (uh,t), state = iter_result

    Uh, (uf, u0, tf, cache) = state # from iterate in ODESolution.jl
    println(uf == u0)               # always true
    ode_cache, vθ, nl_cache = cache # from next_step! in ThetaMethod.jl
    Us, Uts, fecache = ode_cache    # from update_cache! in ODEOperatorInterfaces.jl
    println(length(Us))             # 2 (tuple of TrialFESpace)
    writevtk(Ω, "uh_$t", cellfields=[ "uh_$t"=>Us[1]])  # Method error (ambiguous)

    global iter_result = iterate(uₕₜ, state)
end

The error with writevtk is

ERROR: LoadError: MethodError: CellField(::TrialFESpace{Gridap.FESpaces.UnconstrainedFESpace{Vector{Float64}, Gridap.FESpaces.NodeToDofGlue{Int32}}}, ::Gridap.Geometry.BodyFittedTriangulation{2, 2, CartesianDiscreteModel{2, Float64, typeof(identity)}, CartesianGrid{2, Float64, typeof(identity)}, Gridap.Arrays.IdentityVector{Int64}}) is ambiguous. Candidates:
  CellField(fs::Gridap.FESpaces.SingleFieldFESpace, cell_vals) in Gridap.FESpaces at /home/antoine/.julia/packages/Gridap/TyUsh/src/FESpaces/SingleFieldFESpaces.jl:131
  CellField(fe::FESpace, cell_vals) in Gridap.FESpaces at /home/antoine/.julia/packages/Gridap/TyUsh/src/FESpaces/FESpaceInterface.jl:103
  CellField(f, trian::Triangulation) in Gridap.CellData at /home/antoine/.julia/packages/Gridap/TyUsh/src/CellData/CellFields.jl:103
Possible fix, define
  CellField(::Gridap.FESpaces.SingleFieldFESpace, ::Triangulation)
Stacktrace:

What can I do with this ?

[fverdugo]

@Antoinemarteau

It would be very nice to have uhp=copy(uh) and an in-place copyto!(uhp, uh) for this use :P

This would be indeed very useful if not yet available! Please feel free to add it via a PR.

[Antoinemarteau]

Yes Base.copy has to be specialized. I will try with ::FEFunction this WE!