gridap / GridapODEs.jl

Time stepping for Gridap
MIT License
33 stars 7 forks source link

Error with Gridap@0.14 #37

Closed oriolcg closed 3 years ago

oriolcg commented 3 years ago

Hi @fverdugo,

The following assert breaks when using TransientFESpaces with the latest version of Gridap: https://github.com/gridap/Gridap.jl/blob/5d8ba281f8e96c63235be59a5821c7eda23ac5bf/src/CellData/CellFields.jl#L252

The new developments in Gridap are supposed to work with GridapODEs or there is something that has to be accommodated in GridapODEs?

I don't know if this issue is also related to issue https://github.com/gridap/Gridap.jl/issues/380

fverdugo commented 3 years ago

Do you have a reproducer?

oriolcg commented 3 years ago

Yes, this simple code fails at the same point:

module Reproducer
using Gridap
using GridapODEs
using GridapODEs.TransientFETools

domain = (-1,1,-1,1)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)
trian = Triangulation(model)
quad = CellQuadrature(trian,2)

V = TestFESpace(
  model=model,
  valuetype=VectorValue{2,Float64},
  reffe=:Lagrangian,
  order=2,
  conformity =:H1,
)
U = TransientTrialFESpace(V)

res(t,x,xt,y) = xt⋅y + ∇(x)⋅∇(y)
jac(t,x,xt,dx,y) = ∇(dx)⋅∇(y)
jac_t(t,x,xt,dxt,y) = dxt⋅y

terms = FETerm(res,jac,jac_t,trian,quad)
op = TransientFEOperator(U,V,terms)

end
oriolcg commented 3 years ago

My understanding of this error is that when calling TransientFEOperator(U,V,terms), an HomogeneousTrialFESpace is created for the time derivative of the Trial space U. In the process of creating this new trial space, a TestFESpace is expected, but we are sending a TrialFESpace.

Before, in the TrialFESpace constructor we were calling the function https://github.com/gridap/Gridap.jl/blob/670dabda49345e518541279a71cda09df281ab77/src/FESpaces/TrialFESpaces.jl#L59, which does not check if the input space is TestFESpace. Now we call the function https://github.com/gridap/Gridap.jl/blob/35c4c6377c3b3b9cf85b1d8566730b2f3d43eec5/src/CellData/CellFields.jl#L251, which has the assert that fails.

In the case that concerns this issue, that is creating a new TrialFESpace from a TransientTrialFESpace, I think we shouldn't expect to receive a TestFESpace to create the new TrialFESpace, just create an homogeneous version of the one that's already existing.

What do you think @fverdugo ?

fverdugo commented 3 years ago

Yes, as of Gridap 0.14, you can only construct a TrialFESpace from a test space. The call to trialize_cell_basis in the inner constructor makes the restriction: https://github.com/gridap/Gridap.jl/blob/35c4c6377c3b3b9cf85b1d8566730b2f3d43eec5/src/FESpaces/TrialFESpaces.jl#L8

This was not the case in Gridap 0.13, and therefore the error. We need to recover the old behaviour. It should be trivial:

function trialize_cell_basis(test::CellField)
  if is_trial(test)
    return test 
  end
  @assert is_test(test)
  array = trialize_array_of_bases(get_array(test))
  axs = apply(Fields._add_singleton_block,get_cell_axes(test))
  similar_object(test,array,axs,Val((1,:)))
end

or using traits:

function trialize_cell_basis(test::CellField)
  _trialize_cell_basis(test,MetaSizeStyle(test))
end

function _trialize_cell_basis(test,metasize)
  @unreachable
end

function _trialize_cell_basis(trial,metasize::Val{(1,:)})
  trial
end

function _trialize_cell_basis(test::CellField,metasize::Val{(:,)})
  array = trialize_array_of_bases(get_array(test))
  axs = apply(Fields._add_singleton_block,get_cell_axes(test))
  similar_object(test,array,axs,Val((1,:)))
end

Perhaps the second option is more elegant. @oriolcg can you open a PR with this and with some test that checks that we can indeedd build an HomogeneousTrialFESpace from another TrialFESpace ?

oriolcg commented 3 years ago

Thanks @fverdugo. I'll do the PR with the fix.

fverdugo commented 3 years ago

@oriolcg can we close the issue?

oriolcg commented 3 years ago

Yes, it's solved. Closing issue