How to iterate through transient solutions

[nickbrady]

This is bit similar to the solution iterpolation discussion, but instead of looking at the spatial solution, I am trying to iterate over different temporal solutions.

I am working through one of the notebooks: t017_transient_linear.ipynb

The notebooks predominantly offload the visualization to something like Paraview, which is fine, but I am trying do some of the visualizations, at least for the 2D problems within the Jupyter environment natively, using Julia.

Mostly I have been successful, but I am struggling to visualize solutions at different timesteps of GenericTransientFESolution().

The solution is generated here

t0, tF = 0.0, 10.0
uh0 = interpolate_everywhere(g(t0), Ug(t0))
uh = solve(solver, op, t0, tF, uh0)

The last cell in t017_transient_linear.ipynb is able to iterate through different time steps and save the solutions and I was able to verify in paraview that the solutions mimic the gif that is embedded.

results_tmp = joinpath(results_dir, "tmp")

if !isdir(results_tmp)
  mkdir(results_tmp)
end

createpvd("results") do pvd
  pvd[0] = createvtk(Ω, joinpath(results_tmp, "results_0" * ".vtu"), cellfields=["u" => uh0])
  for (tn, uhn) in uh
    pvd[tn] = createvtk(Ω, joinpath(results_tmp, "results_$tn" * ".vtu"), cellfields=["u" => uhn])
  end
end

Although, uh is iteratable, one cannot access different solution slices using the notation uh[1], so instead I tried to make it into a vector

uhn = collect(uh)

uhn[1] then yields (0.05, SingleFieldFEFunction()).

The problem is that the SingleFieldFEFunction()'s seem to be identical at every time step and all seem to reflect the solution from the final time step (tF = 10.0). I seem to be missing something, any advice would be appreciated.

using Gridap
domain = (-1, +1, -1, +1)
partition = (20, 20)
model = CartesianDiscreteModel(domain, partition)

using IterTools

# Create linspaces for x and y
x = LinRange(domain[1], domain[2], partition[1])
y = LinRange(domain[3], domain[4], partition[2])

# Create grid points as Cartesian product
gridmap = collect(product(x, y))

# Convert tuples to Points
gridmap_points = [Point(p...) for p in gridmap];

# 
gridmap_points_vec = vec(gridmap_points)
x_vec = [pt[1] for pt in gridmap_points_vec]
y_vec = [pt[2] for pt in gridmap_points_vec];

# plotting the heatmap of the solution, uh
using CairoMakie; cm = CairoMakie 

uhn = collect(uh)

f_uh = uhn[1][2]
uh_points_1 = f_uh(gridmap_points_vec);
cm.heatmap(x_vec, y_vec, vec(uh_points_1))

f_uh = uhn[200][2]
uh_points_2 = f_uh(gridmap_points_vec);
cm.heatmap(x_vec, y_vec, vec(uh_points_2))

[nickbrady]

Update: I was able to confirm that uhn[1][2] basically just tracks the solution at tF. I confirmed this by changing the value tF and seeing that it tracked with the corresponding .vtu file result.

[JordiManyer]

From what I understand, we reuse almost everything from one timestep to the next. Otherwise memory would become a problem. This includes the array used to store the dof values of the solution. This is why it is being overwritten.

I would plot the solution within the loop, or do some kind of deep copy. For instance, you could store the free values by doing something like

saved_sols = []
  for (tn, uhn) in uh
     push!(saved_sols,copy(get_free_dof_values(uhn)))
  end

Notice the copy, otherwise you will have the same problem.

[nickbrady]

Thank you, this was helpful.

A follow-up question though because even with your answer, I am not able to accomplish my goal.

saved_sols[1] gives me a 361-element Vector{Float64}, but these are just values, uncoupled to the coordinates so I can't make the heatmap.

I realize the grid is static so doesn't need to be updated at each iteration. I tried to recover the coordinates in two ways:

uhn = collect(uh)
get_cell_dof_values(uhn[1][2])

this gives 400-element Gridap.Arrays.LazyArray.... Not only do the sizes mismatch, but each element also contains 4 values:

4-element Vector{Float64}:
 -0.0
 -0.0
  0.009705231466233432
  0.008639616642907466

so I cannot tell which values are the coordinates and which are the solution values.

another method I tried is:
Gridap.ReferenceFEs.get_node_coordinates(model) which produces 21×21 Matrix{VectorValue{2, Float64}} which does indeed seem to correspond to the coordinates, but the dimensions are not compatible.

I realize the the difference here is likely get_cell_dof_values does not grab the values imposed by the Dirichlet boundary conditions, but I have not been able to resolve how to grab those and then how to properly merge them with the results from get_cell_dof_values.

[JordiManyer]

As the name suggests, get_cell_dof_values(uh) returns the dof values in each cell. If you want the coordinates that go with it, I would suggest get_cell_coordinates(model). Then zip(cell_coords,cell_vals) would return an iterator that for each cell has the pair (coords,dofs).

If you have dirichlet bcs, this will indeed not work. Dirichlet nodes are not free, therefore do not appear as degrees of freedom. I would just interpolate your solution in whatever points you want (uh(pts), where pts is either a vector of points (slow) or a CellPoint object (fast).

If you are looking for plotting capabilities, I would have a look at GridapMakie.jl. If you are not familiar with how Gridap works, I recommend tutorial 13 to get down to the dev level.

[nickbrady]

Ok, thanks for pointing me in the right direction. I will consult those resources and then hopefully accomplish my goal.