Issues developing FSI

[oriolcg]

Here I'll put some issues that I'm facing on the implementation of a transient FSI tutorial, so I can track all the missing points. @fverdugo @santiagobadia feel free to add comments.

• We are not able to integrate terms of the form: q * (∇⋅( ( J(∇(u)) * Finv(∇(u)) )⋅v )) for two main reasons:
  • Second derivatives not ready
  • Chain rule (this could be done directly at the driver level, but we might want a general implementation in Gridap)
    -->Alternative approach: by the Piola identity we have that ∇⋅(JF^{-1}v) = J∇v:F^{-T}
• Entry not found in the sparsity pattern when the temporal jacobian has entries that are not in the spatial jacobian.
  -->Temporary fix: add an entry on the corresponding spatial jacobian multiplied by 0.
• The following error appears when integrating on an interface: ERROR: LoadError: MethodError: no method matching array_cache(::Gridap.Geometry.SkeletonPair{... (still trying to discover the cause) --> Wrong usage of interface terms.

[santiagobadia]

Hi @oriolcg

I would probably not implement chain rule yet. I think it is not that common in FE codes. We can do the analytical work a priori. It would be interesting in the future though.

You could probably try to implement second order derivatives. I don't think it is complicated. @fverdugo is the expert here and can provide guidance.

Can you provide a small code that shows the error in GridapODEs due to the problem you mention and create an issue there? This problem is probably due to the way things are implemented at the compressed array level. Most likely, the fix will be in Gridap anyway.

I think that the problem is here in GridapODEs

function allocate_jacobian(op::TransientFEOperatorFromTerms,uh,cache)
  Uh = evaluate(op.trial,nothing)
  @assert is_a_fe_function(uh)
  du = get_cell_basis(Uh)
  v = get_cell_basis(op.test)
  matdata = collect_cell_jacobian(0.0,uh,uh,du,v,op.terms)
  allocate_matrix(op.assem_t,matdata)
end

We are not using collect_cell_jacobian_t in the allocation of the jacobian. But I am not sure we can do it with the current functionality. @fverdugo what do you think?

[fverdugo]

We are not using collect_cell_jacobian_t in the allocation of the jacobian. But I am not sure we can do it with the current functionality. @fverdugo what do you think?

It can be done (by implementing a small new function).

  matdata_j = collect_cell_jacobian(0.0,uh,uh,du,v,op.terms)
  matdata_jt = collect_cell_jacobian_t(0.0,uh,uh,du,v,op.terms)
  matdata = append_matdata(matdata_j,matdata_jt ) 
  allocate_matrix(op.assem_t,matdata)

You only need to implement append_matdata (perhaps in Gridap).

[fverdugo]

function append_matdata(matdata_j,matdata_jt ) 
term_to_cellmat_j, term_to_cellidsrows_j, term_to_cellidscols_j = matdata_j
term_to_cellmat_jt, term_to_cellidsrows_jt, term_to_cellidscols_jt = matdata_jt

term_to_cellmat = vcat(term_to_cellmat_j,term_to_cellmat_jt)
# idem for the others

matdata = (term_to_cellmat,term_to_cellidsrows, term_to_cellidscols)
matdata
end

[oriolcg]

Hi @fverdugo, I was checking the fix you propose and I have some questions:

1. The arguments that you suggested for the fix on collect_cell_jacobian_t(0.0,uh,uh,du,v,op.terms) do not match the arguments of the implemented function collect_cell_jacobian_t(t::Real,uh,uh_t,du,v,duht_du::Real,terms). I think we need to modify the current allocate_jacobian function such that it takes all the required arguments (as in jacobian_t!), but I'm not sure, could you confirm that?

2. The new append_matdata function should be located in TransientFEOperators.jl, right?

[santiagobadia]

I have done this changes in GridapODEs, but not yet checked it works for your case.
Where is the FSI transient tutorial? I cannot find it in tour GridapFSI repos

[oriolcg]

You can check with the debugging tutorial I have in the transient_fsi branch (https://github.com/gridap/Tutorials/blob/transient_fsi/src/transient_fsi_debug.jl). It's not the latest version because I switched to https://github.com/gridapapps/GridapFSI.jl, but there I already had this issue.

In particular, you can easily chek it by removing the term 0.0*(du⋅ϕ) in

Tutorials/src/transient_fsi_debug.jl

Line 399 in 43caad8

0.0*(du⋅ϕ) - (ϕ⋅dv)