Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Shape derivative #653

Merged
merged 44 commits into from
Feb 27, 2022
Merged

Conversation

ConnorMallon
Copy link

Modifications to allow propogation of duals for autodiff ( to resolve gridap/GridapEmbedded.jl#51 ).

The main problem I had here was with the appended triangulation object. When using a cut triangulation, some of the points come from the background triangulation and some from the cut section of the triangulation. With a pertubation of the geometry, dual numbers appear in the points around the cut areas but the unaffected points are still built on primitive types. Since these points get " apppended " we end up with a mixed type array. I have got around this but it is probably not the cleanest way. Let me know what I can do to make this better. Note that I made the concrete type check look 1 layer deeper to cater for the appended array.

Revert "change to FLoat 32 + cellpoint struct orig"

This reverts commit f4b09ee.

reverting
working for autodiff - needs ironing(checksdeleted

more for allowing piping duals (dont need the SMA)

get rid of show

compat w/ reversediff

get rid assmebler generaliser

type stability

fixing monomial basis to pass tests

get rid comments
@fverdugo
Copy link
Member

fverdugo commented Sep 7, 2021

Hi @ConnorMallon. Test are not passing. Try to pull the latest version in master and push again to this branch.

@ConnorMallon
Copy link
Author

Tests are now passing. Please have a look and let me know what I can do to fix anything or make it a bit more clean.

src/Arrays/AppendedArrays.jl Outdated Show resolved Hide resolved
@ConnorMallon
Copy link
Author

Hey @fverdugo, What is purpose of doing a collect in the line T = eltype(collect((ai,bi))) when we are going to have a type stability issue if ai and bi are of different types? In the latest commit, I just take the element type from one of the objects and it seems to work ( all tests pass and it works for doing a shape derivative ) and in that case there is no type instability.

@ConnorMallon
Copy link
Author

ConnorMallon commented Sep 22, 2021

I have smuggled in some more changes here. These are for something slightly different but still related. I noticed an issue when trying to use Gridaps built in Autodiff on a bimaterial problem. Here we have a problem where a FEFunction only exists in certain parts of the domain. It means that if we do something like get_cell_dof_values we get back an array for which some of the entries are empty. The issue arises then when we try to find the element type of this object. Since the type inference occurs by querying the type of the first element then we might get back the type representing an empty array in the case which the FEFunction does not exist in the first cell. This then causes problems down the line. I fix this by avoiding relying on the lazy_map type inference and instead doing something more manual. Here is a MWE which works only after the addition of the last few commits:

module BimaterialCutFEMTestsModifiedforMWE
using Gridap
import Gridap: ∇
using GridapEmbedded
R = 0.7
geo1 = disk(R)
n = 30
domain = (-1,1,-1,1)
partition = (n,n)
bgmodel = simplexify(CartesianDiscreteModel(domain,partition))
cutgeo = cut(bgmodel,geo1)
model1 = DiscreteModel(cutgeo,geo1)
Ω1 = Triangulation(cutgeo,geo1)
order = 1
degree = 2*order
dΩ1 = Measure(Ω1,degree)
V1 = TestFESpace(model1,ReferenceFE(lagrangian,Float64,order),conformity=:H1)
uh = interpolate(1,V1)
a(uh) =∫(uh)dΩ1
∇(a)(uh)
end # module

@ConnorMallon
Copy link
Author

ConnorMallon commented Sep 27, 2021

Also added some support for multifield but the following MWE does not work:

module BimaterialPoissonCutFEMTests

using Gridap
import Gridap: ∇
using GridapEmbedded
using Test

# Select geometry
const R = 0.7
geo1 = disk(R)
geo2 =  !geo1
n = 10
domain = (-1,1,-1,1)
partition = (n,n)

# Setup background model
bgmodel = simplexify(CartesianDiscreteModel(domain,partition))

# Cut the background model
cutgeo = cut(bgmodel,union(geo1,geo2))

# Setup models
model1 = DiscreteModel(cutgeo,geo1)
model2 = DiscreteModel(cutgeo,geo2)

# Setup integration meshes
Ω1 = Triangulation(cutgeo,geo1)
Ω2 = Triangulation(cutgeo,geo2)

# Setup Lebesgue measures
order = 1
degree = 2*order
dΩ1 = Measure(Ω1,degree)
dΩ2 = Measure(Ω2,degree)

# Setup FESpace

V1 = TestFESpace(model1,ReferenceFE(lagrangian,Float64,order),conformity=:H1)
V2 = TestFESpace(model2,ReferenceFE(lagrangian,Float64,order),conformity=:H1)

U1 = TrialFESpace(V1)
U2 = TrialFESpace(V2)

V = MultiFieldFESpace([V1,V2])
U = MultiFieldFESpace([U1,U2])

uh=FEFunction(V,ones(num_free_dofs(V)))
j((u1,u2)) = ∫(u1)dΩ1 + ∫(u2)dΩ2
djdu=∇(j)(uh)
assemble_vector(djdu,V)

end # module

@ConnorMallon
Copy link
Author

ConnorMallon commented Sep 27, 2021

Reverted changes related to Gridaps Autodiff (derivative w.r.t a FEfunction) . This PR is now only for making the shape derivative possible (derivative w.r.t the goemetry)

@fverdugo
Copy link
Member

And this seems what @fverdugo wanted to avoid in his initial comment here: #653 (comment)

Thus, now, I am little bit lost/confused. @fverdugo can you please help to clarify the situation?

Yes, AppendedArray is general enough to handle cases with different eltypes. But, even in this situation, it is preferable to feed it with two arrays with the same eltype. You can evaluate with a profiler if this is a hotspot or not and then decide.

For me, the only important request in this PR is to not introduce type instabilities in the code that was already working. If you end up with changes that introduce type instabilities for your new use case, but you do not introduce type instabilities for the previous use cases, then it is rather acceptable.

@amartinhuertas
Copy link
Member

Yes, AppendedArray is general enough to handle cases with different eltypes. But, even in this situation, it is preferable to feed it with two arrays with the same eltype. You can evaluate with a profiler if this is a hotspot or not and then decide.
For me, the only important request in this PR is to not introduce type instabilities in the code that was already working. If you end up with changes that introduce type instabilities for your new use case, but you do not introduce type instabilities for the previous use cases, then it is rather acceptable.

Ok, thanks for your answer. We are on the same page. The question that I have now is why @ConnorMallon used Union instead of eltype((ai,bi)). Is there any reason @ConnorMallon ?

@ConnorMallon
Copy link
Author

ConnorMallon commented Feb 24, 2022

Ok, thanks for your answer. We are on the same page. The question that I have now is why @ConnorMallon used Union instead of eltype((ai,bi)). Is there any reason @ConnorMallon ?

Because when a was a vector of floats and b a vector of duals, eltype((ai,bi)) was not giving me back a concrete type. it was giving me for the AppendedArray something like Vector{T} where T. Thats why i changed it to promote_type where we can explicitly define which type to choose. Union also solves this problem but as I have learnt it introduces type instability.

So in essence, I was trying to avoid the abstract types given by eltype((ai,bi)) as you mentioned.

@amartinhuertas
Copy link
Member

amartinhuertas commented Feb 24, 2022

So in essence, I was trying to avoid the abstract types given by eltype((ai,bi)) as you mentioned.

Ok. Thus, it was not actually something that you needed to avoid a code crash or something like that, right? In such a case, let use live with the current solution. It does not introduce type instability for the most common case in which both arrays are of the same eltype.

@ConnorMallon
Copy link
Author

ConnorMallon commented Feb 24, 2022

So in essence, I was trying to avoid the abstract types given by eltype((ai,bi)) as you mentioned.

Ok. Thus, it was not actually something that you needed to avoid a code crash or something like that, right? In such a case, let use live with the current solution. It does not introduce type instability for the most common case in which both arrays are of the same eltype.

no it was actually causing an error down the line somewhere T not defined if i remember correctly. I will reproduce it and put the whole stacktrace.

@amartinhuertas
Copy link
Member

no it was actually causing an error down the line somewhere T not defined if i remember correctly. I will reproduce it and put the whole stacktrace.

Ok!!!!!!! At last!!!! :-) Everything starts making sense to me now. I was wondering how you realized of this very subtle detail without an error.

@ConnorMallon
Copy link
Author

Ok!!!!!!! At last!!!! :-) Everything starts making sense to me now. I was wondering how you realized of this very subtle detail without an error.

Yes it would of probably been more helpful if I provided a bit more explanation on each of the commits explaining why exactly the changes were made including the errors in which they were made to get around.

@ConnorMallon
Copy link
Author

Hey @amartinhuertas, I have made the changes as requested in the last commit.

@amartinhuertas
Copy link
Member

Hey @amartinhuertas, I have made the changes as requested in the last commit.

Ok. Thanks for your work. We have resolved one of the two main issues in this PR (frequent dynamic memory allocation of small arrays). There is still the other that has to be solved. I will wait for your reproducer.

@ConnorMallon
Copy link
Author

ConnorMallon commented Feb 24, 2022

Here is the error:

ERROR: LoadError: MethodError: no method matching Gridap.CellData.CellPoint(::Gridap.Arrays.AppendedArray{Vector{Gridap.TensorValues.VectorValue{2, Float64}}, Gridap.Arrays.CompressedArray{Vector{Gridap.TensorValues.VectorValue{2, Float64}}, 1, Vector{Vector{Gridap.TensorValues.VectorValue{2, Float64}}}, Vector{Int8}}, FillArrays.Fill{Vector{Gridap.TensorValues.VectorValue{2, Float64}}, 1, Tuple{Base.OneTo{Int64}}}}, ::Gridap.Arrays.AppendedArray{Vector,Gridap.Arrays.LazyArray{FillArrays.Fill{Gridap.Arrays.Broadcasting{Gridap.Arrays.Reindex{Vector{Gridap.TensorValues.VectorValue{2, ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}}}}}}, 1, Tuple{Base.OneTo{Int64}}}, Vector{Gridap.TensorValues.VectorValue{2, ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}}}}, 1, Tuple{Gridap.Arrays.Table{Int32, Vector{Int32}, Vector{Int32}}}}, Gridap.Arrays.LazyArray{FillArrays.Fill{Gridap.Arrays.Broadcasting{Gridap.Arrays.Reindex{Gridap.Geometry.CartesianCoordinates{2, Float64, typeof(identity)}}}, 1, Tuple{Base.OneTo{Int64}}}, Vector{Gridap.TensorValues.VectorValue{2, Float64}}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Gridap.Arrays.Reindex{Gridap.Geometry.CartesianCellNodes{2}}, 1, Tuple{Base.OneTo{Int64}}}, Vector{Int32}, 1, Tuple{Vector{Int32}}}}}}, ::Gridap.Geometry.AppendedTriangulation{2, 2, GridapEmbedded.Interfaces.SubCellTriangulation{2, 2, ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}}, Gridap.Geometry.CartesianDiscreteModel{2, Float64, typeof(identity)}}, Gridap.Geometry.BodyFittedTriangulation{2, 2, Gridap.Geometry.CartesianDiscreteModel{2, Float64, typeof(identity)}, Gridap.Geometry.GridPortion{2, 2, Gridap.Geometry.CartesianGrid{2, Float64, typeof(identity)}}, Vector{Int32}}}, ::Gridap.CellData.ReferenceDomain)
Closest candidates are:
  Gridap.CellData.CellPoint(::AbstractArray{<:Union{AbstractArray{<:Gridap.TensorValues.VectorValue}, Gridap.TensorValues.VectorValue}}, ::AbstractArray{<:Union{AbstractArray{<:Gridap.TensorValues.VectorValue}, Gridap.TensorValues.VectorValue}}, ::Gridap.Geometry.Triangulation, ::DS) where DS at ~/.julia/dev/Gridap/src/CellData/CellFields.jl:6
  Gridap.CellData.CellPoint(::AbstractArray{<:Union{AbstractArray{<:Gridap.TensorValues.VectorValue}, Gridap.TensorValues.VectorValue}}, ::Gridap.Geometry.Triangulation, ::Gridap.CellData.PhysicalDomain) at ~/.julia/dev/Gridap/src/CellData/CellFields.jl:22
  Gridap.CellData.CellPoint(::AbstractArray{<:Union{AbstractArray{<:Gridap.TensorValues.VectorValue}, Gridap.TensorValues.VectorValue}}, ::Gridap.Geometry.Triangulation, ::Gridap.CellData.ReferenceDomain) at ~/.julia/dev/Gridap/src/CellData/CellFields.jl:12

Note the

::Gridap.Arrays.AppendedArray{
- Vector,
Gridap.Arr............................}

buried in there (the second argument we are trying to give to the function). It is a vector without any eltype defined which results from trying to do T = eltype(collect((ai,bi)))

@amartinhuertas
Copy link
Member

amartinhuertas commented Feb 24, 2022

Here is the error:

For clarity in the discussion, let me separate the 4 arguments to the CellField constructor into separate boxes

1st argument

::Gridap.Arrays.AppendedArray{Vector{Gridap.TensorValues.VectorValue{2, Float64}},
Gridap.Arrays.CompressedArray{Vector{Gridap.TensorValues.VectorValue{2, Float64}}, 1, Vector{Vector{Gridap.TensorValues.VectorValue{2, Float64}}}, Vector{Int8}}, FillArrays.Fill{Vector{Gridap.TensorValues.VectorValue{2, Float64}}, 1, Tuple{Base.OneTo{Int64}}}},

2nd argument

::Gridap.Arrays.AppendedArray{Vector,Gridap.Arrays.LazyArray{FillArrays.Fill{Gridap.Arrays.Broadcasting{Gridap.Arrays.Reindex{Vector{Gridap.TensorValues.VectorValue{2, ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}}}}}}, 1, Tuple{Base.OneTo{Int64}}}, Vector{Gridap.TensorValues.VectorValue{2, ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}}}}, 1, Tuple{Gridap.Arrays.Table{Int32, Vector{Int32}, Vector{Int32}}}}, Gridap.Arrays.LazyArray{FillArrays.Fill{Gridap.Arrays.Broadcasting{Gridap.Arrays.Reindex{Gridap.Geometry.CartesianCoordinates{2, Float64, typeof(identity)}}}, 1, Tuple{Base.OneTo{Int64}}}, Vector{Gridap.TensorValues.VectorValue{2, Float64}}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Gridap.Arrays.Reindex{Gridap.Geometry.CartesianCellNodes{2}}, 1, Tuple{Base.OneTo{Int64}}}, Vector{Int32}, 1, Tuple{Vector{Int32}}}}}},

3rd argument

::Gridap.Geometry.AppendedTriangulation{2, 2, GridapEmbedded.Interfaces.SubCellTriangulation{2, 2, ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}}, Gridap.Geometry.CartesianDiscreteModel{2, Float64, typeof(identity)}}, Gridap.Geometry.BodyFittedTriangulation{2, 2, Gridap.Geometry.CartesianDiscreteModel{2, Float64, typeof(identity)}, Gridap.Geometry.GridPortion{2, 2, Gridap.Geometry.CartesianGrid{2, Float64, typeof(identity)}}, Vector{Int32}}})

4th argument

::Gridap.CellData.ReferenceDomain

These 4th arguments must match the type annotations of the member variables in the type definition statement

struct CellPoint{DS} <: CellDatum
  cell_ref_point::AbstractArray{<:Union{Point,AbstractArray{<:Point}}}
  cell_phys_point::AbstractArray{<:Union{Point,AbstractArray{<:Point}}}
  trian::Triangulation
  domain_style::DS
end

Vector <: Union{Point,AbstractArray{<:Point}} is false in the cell_phys_point member variable, and this causes the error.

@amartinhuertas
Copy link
Member

Now that the problem and its cause are clear, we can start talking about possible solutions.

@fverdugo ... can we remove the constraint from the CellPoint member variables? This is one possible solution.

@amartinhuertas
Copy link
Member

@fverdugo ... can we remove the constraint from the CellPoint member variables? This is one possible solution.

@ConnorMallon ... if you remove the constraint from the CellPoint member variables, does your code run smoothly, or there are other issues to be solved?

@ConnorMallon
Copy link
Author

ConnorMallon commented Feb 24, 2022

@ConnorMallon ... if you remove the constraint from the CellPoint member variables, does your code run smoothly, or there are other issues to be solved?

It does run smoothly. However,I am going to guess that this change is going to break stability for regular Gridap functions that were type stable. I have tried to avoid this approach after being warned about this by francesc here

@amartinhuertas
Copy link
Member

Ok, I was probably loose when using the expression remove the constraint. By this I mean to annotate cell_ref and cell_phys member variables with Tr, and Tp types that become type parameters right after dc in the type declaration. I don't see why this introdces type instability. Indeed the current member variables are polymorphic abstractarray

@amartinhuertas
Copy link
Member

Ok, I was probably loose when using the expression remove the constraint. By this I mean to annotate cell_ref and cell_phys member variables with Tr, and Tp types that become type parameters right after dc in the type declaration. I don't see why this introdces type instability. Indeed the current member variables are polymorphic abstractarray

HI @ConnorMallon ! I have implemented this here ConnorMallon@56dd373

Can you please double check that your code runs with this mod? If affirmative, and the tests pass, I will accept this PR.

@amartinhuertas
Copy link
Member

and the tests pass

@ConnorMallon tests are passing, good news. Once you confirm that your code runs succesfully, I will accept the PR.

@ConnorMallon
Copy link
Author

Once you confirm that your code runs succesfully, I will accept the PR

code runs successfully. thanks Alberto

@amartinhuertas amartinhuertas merged commit 36bc2a7 into gridap:master Feb 27, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Shape derivative
3 participants