gridap / GridapODEs.jl

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

Added TransientTrialFreeFESpace and corresponding tests #42

Closed oriolcg closed 3 years ago

oriolcg commented 3 years ago

I defined TransientTrialFESpace as an abstract type with the subtypes TransientTrialDirichletFESpace and TransientTrialFreeFESpace. The former corresponds to the original TransientTrialFESpace and the later implements the functions that were assuming the presence of Dirichlet data.

These changes are required to define transient FE spaces without strong Dirichlet data, which is relevant when only imposing weak boundary conditions.

Closes #41

codecov-io commented 3 years ago

Codecov Report

Merging #42 into master will increase coverage by 0.03%. The diff coverage is 100.00%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master      #42      +/-   ##
==========================================
+ Coverage   63.06%   63.09%   +0.03%     
==========================================
  Files          18       18              
  Lines         666      672       +6     
==========================================
+ Hits          420      424       +4     
- Misses        246      248       +2     
Impacted Files Coverage Δ
src/TransientFETools/TransientFESpaces.jl 83.60% <100.00%> (-1.85%) :arrow_down:

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update beb8148...fd00933. Read the comment docs.

fverdugo commented 3 years ago

@oriolcg

the old FE term with a 0-lenght vector of dirichlet dofs is not enough?

oriolcg commented 3 years ago

@oriolcg

the old FE term with a 0-lenght vector of dirichlet dofs is not enough?

I tried to send an empty vector of functions but then it fails at https://github.com/gridap/GridapODEs.jl/blob/beb8148f5c83bcd7bc6e80cb966b454ce02587cf/src/TransientFETools/TransientFESpaces.jl#L27

Is that what you mean?

fverdugo commented 3 years ago

Is that what you mean?

The the case without strong Dirichlet BCs is conceptually just a particular case with 0 Dirichlet DOFs. The code has to be able to handle this, if not, it is telling that there is a small design flaw or a bug. The steady FE spaces are able to deal with 0 Dirichlet DOFs. The transient should be able as well.

I propose to try to understand why we have a problem in TransientTrialFESpace with this particular case and try to fix it.

santiagobadia commented 3 years ago

@oriolcg I am not sure what you say is possible.

In the previous version, you had a HomogeneousTrialFESpace, which is not a sub-type of TransientTrialFESpace and thus, you cannot reach that line in which you say there is an error.

The idea of the original design is that a TrialFESpace already implements the API of the TransientTrialFESpace so you should not need to do the dev in this pull request.

There can be bugs, and we can fix them, but I am not sure we want to add this unless there is a compelling reason. For the TrialFESpace we just do trivial (cheap) implementations of the transient ones.

fverdugo commented 3 years ago

As @santiagobadia you don't need a TrialTransientFESpace for this case since the only time-dependent data in this object are the dirichlet functions.

In any case, I see that it can be problematic to call TrialFESpace!(Ut,objects_at_t) with a vector of empty Dirichlet functions (I have never tested it). If this is the problem, it should be easy to fix it in the Gridap repo.

oriolcg commented 3 years ago

Thanks for the comments @santiagobadia and @fverdugo! I have a better picture now.

(Santi) In the previous version, you had a HomogeneousTrialFESpace, which is not a sub-type of TransientTrialFESpace and thus, you cannot reach that line in which you say there is an error.

True, in the previous version it was not failing there. It was failing at https://github.com/gridap/GridapODEs.jl/blob/beb8148f5c83bcd7bc6e80cb966b454ce02587cf/src/TransientFETools/TransientFESpaces.jl#L70 where it was expected a TransientTrialFESpace and a multifield composed of TrialFESpace (the result of HomogeneousTrialFESpace) was received.

ERROR: MethodError: no method matching time_derivative(::Gridap.MultiField.MultiFieldFESpace{Gridap.MultiField.ConsecutiveMultiFieldStyle,false})
Closest candidates are:
  time_derivative(::GridapODEs.TransientFETools.TransientMultiFieldTrialFESpace) at C:\Users\ocolomesgene\.julia\packages\GridapODEs\JBOmm\src\TransientFETools\TransientFESpaces.jl:153
  time_derivative(::GridapODEs.TransientFETools.TransientTrialFESpace) at C:\Users\ocolomesgene\.julia\packages\GridapODEs\JBOmm\src\TransientFETools\TransientFESpaces.jl:70
  time_derivative(::Function) at C:\Users\ocolomesgene\.julia\packages\GridapODEs\JBOmm\src\ODETools\DiffOperators.jl:1

I believed that this is not the intended behavior and that a transient trial FE space should always be a TransientTrialFESpace, even for the homogeneous case. This is the original reasoning behind this pull request.

(Santi) The idea of the original design is that a TrialFESpace already implements the API of the TransientTrialFESpace so you should not need to do the dev in this pull request. There can be bugs, and we can fix them, but I am not sure we want to add this unless there is a compelling reason. For the TrialFESpace we just do trivial (cheap) implementations of the transient ones.

I see your point.

(Francesc) In any case, I see that it can be problematic to call TrialFESpace!(Ut,objects_at_t) with a vector of empty Dirichlet functions (I have never tested it). If this is the problem, it should be easy to fix it in the Gridap repo.

Yes, that is the main problem. Based on your comments I propose to: 1) Define the constructor TransientTrialFESpace without Dirichlet data as:

function TransientTrialFESpace(space::SingleFieldFESpace)
  emptyVector = Function[]
  TransientTrialFESpace(space, emptyVector)
end

2.1) Modify the evaluate! function to handle this case

function evaluate!(Ut::TrialFESpace,U::TransientTrialFESpace,t::Real)
  if isempty(U.dirichlet_t)
    return Ut
  end
  if isa(U.dirichlet_t,Vector)
    objects_at_t = map( o->o(t), U.dirichlet_t)
  else
    objects_at_t = U.dirichlet_t(t)
  end
  TrialFESpace!(Ut,objects_at_t)
  Ut
end

or 2.2) Modify TrialFESpace! to handle this case

function TrialFESpace!(space::TrialFESpace,objects)
  if isempty(objects)
    dir_values = zero_dirichlet_values(space)
  else
    dir_values = get_dirichlet_values(space)
  end
  dir_values = compute_dirichlet_values_for_tags!(dir_values,space,objects)
  space
end

What do you thing is the best option? Is there a better way?

fverdugo commented 3 years ago

@oriolcg

allowing to pass an empty vector to TrialFESpace and TrialFESpace! in Gridap is ok for me. In fact, I would fix this independently of how we decide to solve this PR, since it is a separate issue.

Some remarks:

fverdugo commented 3 years ago

@oriolcg another remark:

I would keep TrialFESpace! as it is now and fix the problem in the lower level function compute_dirichlet_values_for_tags!.

santiagobadia commented 3 years ago

I would pick 2.2 and change Gridap. If you agree, we can close this PR and the issue.

oriolcg commented 3 years ago

@santiagobadia @fverdugo, going deeper into the code I've realized that TrialFESpace can handle an empty vector of functions, so there is no need for modifications in Gridap. Moreover, I've seen that with the latest version Gridap#master the original error that motivated this PR does not show up. Maybe due to some modifications in MultiFieldFESpace?

If you agree, I'll close the issue and PR and work with the master branch in my application.