Added a new, more general inner constructor to CartesianDiscreteModel #245

amartinhuertas · 2020-05-06T13:18:53Z

Add the following functionality into Gridap with the following API:

CartesianDiscreteModel(gdesc::CartesianDescriptor) # Existing API
CartesianDiscreteModel(gdesc::CartesianDescriptor,cmin::CartesianIndex,cmax::CartesianIndex) # New API

This requirement arised from gridap/GridapDistributed.jl#10

Please note that once we merge this PR into Gridap.jl master, we should pull from master into GridapDistributed branch

fverdugo · 2020-05-06T14:05:42Z

src/Geometry/CartesianDiscreteModels.jl

+function _find_ncube_face_neighbor_delta(p::Polytope{D}, face_lid) where {D}
+  @assert is_n_cube(p)
+  result = fill(0, D)
+  face = p.dface.nfaces[face_lid]


This function is not valid for any Polytope, only for ExtrussionPolytope (note that you are accessing the field dface)

Change the function signature to

function _find_ncube_face_neighbor_delta(p::ExtrusionPolytope{D}, face_lid) where {D}

and implement a fallback function

function _find_ncube_face_neighbor_delta(p::Polytope{D}, face_lid) where {D} @notimplemented "Only implemented for ExtrusionPolytope" end

Please check if other parts of the new code are also specific for ExtrusionPolytope

This function is not valid for any Polytope, only for ExtrussionPolytope (note that you are accessing the field dface)

Ok. Got it!

Change the function signature to

function _find_ncube_face_neighbor_delta(p::ExtrusionPolytope{D}, face_lid) where {D}

Done. e4eda3c

and implement a fallback function

function _find_ncube_face_neighbor_delta(p::Polytope{D}, face_lid) where {D} @notimplemented "Only implemented for ExtrusionPolytope" end

I think there is no need to implemented a fallback function. Please note that the functions (now) restricted to ExtrusionPolytope{D} I have added are helper functions which are only used in the domain space of CartesianDiscreteModels.jl. In other words, these are NOT methods bound to Polytope{D}'s API. Agreed?

Please check if other parts of the new code are also specific for ExtrusionPolytope

Did it. The aforementioned two helper functions are the only ones which are specific to ExtrusionPolytope.

type ExtrusionPolytope{D}

codecov-io · 2020-05-07T00:54:03Z

Codecov Report

Merging #245 into master will increase coverage by 0.08%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master     #245      +/-   ##
==========================================
+ Coverage   88.77%   88.86%   +0.08%     
==========================================
  Files         146      146              
  Lines        9096     9169      +73     
==========================================
+ Hits         8075     8148      +73     
  Misses       1021     1021

Impacted Files	Coverage Δ
src/Geometry/CartesianDiscreteModels.jl	`100.00% <100.00%> (ø)`

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 07bd53c...a9e3d46. Read the comment docs.

fverdugo · 2020-05-07T05:26:42Z

src/Geometry/CartesianDiscreteModels.jl

+    cell_to_faces = get_faces(topo, D, d)
+    face_to_geolabel = face_labeling.d_to_dface_to_entity[d+1]
+    nfaces = length(face_to_geolabel)
+    for face_gid = 1:nfaces


@amartinhuertas Could you put long loops like this one (i.e. loops over cells or over faces) inside a function barrier?
This is a way to ensure type stability in computationally critical parts

Just to double-check, you are asking that I move the loop in lines 272-315 to a separate function, right? I will definitely do it.

On the other hand, while I clearly understand the issue in the Julia Manual's example, I fail to see why here there is a type stability issue (I did not check it with @code_typed yet, I will do it). Where does the type instability stems from? face_to_cells, cell_to_faces ? anything else?

The problem is that some lines at the beginning of the function

function _fill_subgrid_cartesian_entities!(labels, topo, subdesc, desc, cmin)

are potentially type unstable, like

polytope = first(get_polytopes(topo))

For performance reasons, a common pattern in Gridap is to split functions into a small type unstable pre-process and a type-stable computationally expensive part.

The problem is that some lines at the beginning of the function
function _fill_subgrid_cartesian_entities!(labels, topo, subdesc, desc, cmin)
are potentially type unstable, like
polytope = first(get_polytopes(topo))

Ok. I see now the point after using the @code_warntype macro:

julia> @code_warntype Gridap.Geometry._fill_subgrid_cartesian_entities!(model.face_labeling, model.grid_topology, desc, desc, CartesianIndex(2,2,2)) Variables #self#::Core.Compiler.Const(Gridap.Geometry._fill_subgrid_cartesian_entities!, false) labels::Gridap.Geometry.FaceLabeling topo::Gridap.Geometry.UnstructuredGridTopology{3,3,Float64,true} subdesc::Gridap.Geometry.CartesianDescriptor{3,Float64,typeof(identity)} desc::Gridap.Geometry.CartesianDescriptor{3,Float64,typeof(identity)} cmin::CartesianIndex{3} D::Int64 d_to_dface_to_entity::Array{Array{Int32,1},1} gcis::CartesianIndices{3,Tuple{Base.OneTo{Int64},Base.OneTo{Int64},Base.OneTo{Int64}}} subcis::CartesianIndices{3,Tuple{Base.OneTo{Int64},Base.OneTo{Int64},Base.OneTo{Int64}}} polytope::Any face_labeling::Gridap.Geometry.FaceLabeling offsets::Any interior_id::Any boundary_id::Int64 minus_one_ci::CartesianIndex{3} polytope_d_face_to_jfaces::Array{Array{Array{Int64,1},1},2} @_18::Union{Nothing, Tuple{Int64,Int64}} face_deltas::Array{CartesianIndex,1} @_20::Union{Nothing, Tuple{Int64,Int64}} d@_21::Int64 @_22::Union{Nothing, Tuple{Int64,Int64}} j@_23::Int64 face_to_cells::Gridap.Arrays.Table{Int64,Int32} cell_to_faces::Gridap.Arrays.Table{Int64,Int32} face_to_geolabel::Array{Int32,1} nfaces::Int64 @_28::Union{Nothing, Tuple{Int64,Int64}} d@_29::Int64 face_gid::Int64 cell_gid::Int64 a::Int32 b::Int64 gci::CartesianIndex{3} @_35::Union{Nothing, Tuple{Int64,Int64}} face_lid::Any is_assigned_face_delta::Bool cell_found::Bool @_39::Union{Nothing, Tuple{Int64,Int64}} j@_40::Int64 dface_to_jfaces::Any j@_42::Int64 @_43::Union{Missing, Bool} Body::Array{Int32,1} 1 ── Core.NewvarNode(:(face_deltas)) │ Core.NewvarNode(:(@_20)) │ (D = Gridap.Geometry.num_cell_dims(topo)) │ (d_to_dface_to_entity = Base.getproperty(labels, :d_to_dface_to_entity)) │ %5 = Base.getproperty(desc, :partition)::Gridap.TensorValues.MultiValue{Tuple{3},Int64,1,3} │ %6 = Gridap.Geometry.Tuple(%5)::Tuple{Int64,Int64,Int64} │ (gcis = Gridap.Geometry.CartesianIndices(%6)) │ %8 = Base.getproperty(subdesc, :partition)::Gridap.TensorValues.MultiValue{Tuple{3},Int64,1,3} │ %9 = Gridap.Geometry.Tuple(%8)::Tuple{Int64,Int64,Int64} │ (subcis = Gridap.Geometry.CartesianIndices(%9)) │ %11 = Gridap.Geometry.get_polytopes(topo)::Array{_A,1} where _A │ (polytope = Gridap.Geometry.first(%11)) │ (face_labeling = labels) │ (offsets = Gridap.Geometry.get_offsets(polytope)) │ (interior_id = Gridap.Geometry.num_faces(polytope)) │ (boundary_id = -1) │ %17 = Core.apply_type(Gridap.Geometry.Val, D::Core.Compiler.Const(3, false))::Core.Compiler.Const(Val{3}, false) │ %18 = (%17)()::Core.Compiler.Const(Val{3}(), false) │ %19 = Gridap.Geometry.tfill(-1, %18)::Core.Compiler.Const((-1, -1, -1), false) │ (minus_one_ci = Gridap.Geometry.CartesianIndex(%19)) │ %21 = Core.apply_type(Gridap.Geometry.Vector, Gridap.Geometry.Int)::Core.Compiler.Const(Array{Int64,1}, false) │ %22 = Core.apply_type(Gridap.Geometry.Vector, %21)::Core.Compiler.Const(Array{Array{Int64,1},1}, false) │ %23 = Core.apply_type(Gridap.Geometry.Matrix, %22)::Core.Compiler.Const(Array{Array{Array{Int64,1},1},2}, false) │ %24 = Core.tuple(D::Core.Compiler.Const(3, false), D::Core.Compiler.Const(3, false))::Core.Compiler.Const((3, 3), false) │ (polytope_d_face_to_jfaces = (%23)(Gridap.Geometry.undef, %24)) │ %26 = (D::Core.Compiler.Const(3, false) - 1)::Core.Compiler.Const(2, false) │ %27 = (0:%26)::Core.Compiler.Const(0:2, false) │ (@_18 = Base.iterate(%27)) │ %29 = (@_18::Core.Compiler.Const((0, 0), false) === nothing)::Core.Compiler.Const(false, false) │ %30 = Base.not_int(%29)::Core.Compiler.Const(true, false) └─── goto #7 if not %30 2 ┄─ %32 = @_18::Tuple{Int64,Int64}::Tuple{Int64,Int64} │ (d@_21 = Core.getfield(%32, 1)) │ %34 = Core.getfield(%32, 2)::Int64 │ %35 = (d@_21 + 1)::Int64 │ %36 = (D::Core.Compiler.Const(3, false) - 1)::Core.Compiler.Const(2, false) │ %37 = (%35:%36)::UnitRange{Int64} │ (@_22 = Base.iterate(%37)) │ %39 = (@_22 === nothing)::Bool │ %40 = Base.not_int(%39)::Bool └─── goto #5 if not %40 3 ┄─ %42 = @_22::Tuple{Int64,Int64}::Tuple{Int64,Int64} │ (j@_23 = Core.getfield(%42, 1)) │ %44 = Core.getfield(%42, 2)::Int64 │ %45 = Gridap.Geometry.get_faces(polytope, d@_21, j@_23)::Any │ %46 = polytope_d_face_to_jfaces::Array{Array{Array{Int64,1},1},2} │ %47 = (d@_21 + 1)::Int64 │ %48 = (j@_23 + 1)::Int64 │ Base.setindex!(%46, %45, %47, %48) │ (@_22 = Base.iterate(%37, %44)) │ %51 = (@_22 === nothing)::Bool │ %52 = Base.not_int(%51)::Bool └─── goto #5 if not %52 4 ── goto #3 5 ┄─ (@_18 = Base.iterate(%27, %34)) │ %56 = (@_18 === nothing)::Bool │ %57 = Base.not_int(%56)::Bool └─── goto #7 if not %57 6 ── goto #2 7 ┄─ (face_deltas = Gridap.Geometry._find_ncube_face_neighbor_deltas(polytope)) │ %61 = (D::Core.Compiler.Const(3, false) - 1)::Core.Compiler.Const(2, false) │ %62 = (0:%61)::Core.Compiler.Const(0:2, false) │ (@_20 = Base.iterate(%62)) │ %64 = (@_20::Core.Compiler.Const((0, 0), false) === nothing)::Core.Compiler.Const(false, false) │ %65 = Base.not_int(%64)::Core.Compiler.Const(true, false) └─── goto #35 if not %65 8 ┄─ %67 = @_20::Tuple{Int64,Int64}::Tuple{Int64,Int64} │ (d@_29 = Core.getfield(%67, 1)) │ %69 = Core.getfield(%67, 2)::Int64 │ (face_to_cells = Gridap.Geometry.get_faces(topo, d@_29, D::Core.Compiler.Const(3, false))) │ (cell_to_faces = Gridap.Geometry.get_faces(topo, D::Core.Compiler.Const(3, false), d@_29)) │ %72 = Base.getproperty(face_labeling, :d_to_dface_to_entity)::Array{Array{Int32,1},1} │ %73 = (d@_29 + 1)::Int64 │ (face_to_geolabel = Base.getindex(%72, %73)) │ (nfaces = Gridap.Geometry.length(face_to_geolabel)) │ %76 = (1:nfaces)::Core.Compiler.PartialStruct(UnitRange{Int64}, Any[Core.Compiler.Const(1, false), Int64]) │ (@_28 = Base.iterate(%76)) │ %78 = (@_28 === nothing)::Bool │ %79 = Base.not_int(%78)::Bool └─── goto #33 if not %79 9 ┄─ Core.NewvarNode(:(is_assigned_face_delta)) │ Core.NewvarNode(:(cell_found)) │ Core.NewvarNode(:(@_39)) │ %84 = @_28::Tuple{Int64,Int64}::Tuple{Int64,Int64} │ (face_gid = Core.getfield(%84, 1)) │ %86 = Core.getfield(%84, 2)::Int64 │ %87 = Base.getproperty(face_to_cells, :data)::Array{Int64,1} │ %88 = Base.getproperty(face_to_cells, :ptrs)::Array{Int32,1} │ %89 = Base.getindex(%88, face_gid)::Int32 │ (cell_gid = Base.getindex(%87, %89)) │ %91 = Base.getproperty(cell_to_faces, :ptrs)::Array{Int32,1} │ (a = Base.getindex(%91, cell_gid)) │ %93 = Base.getproperty(cell_to_faces, :ptrs)::Array{Int32,1} │ %94 = (cell_gid + 1)::Int64 │ %95 = Base.getindex(%93, %94)::Int32 │ (b = %95 - 1) │ %97 = Base.getindex(subcis, cell_gid)::CartesianIndex{3} │ %98 = (%97 + minus_one_ci::Core.Compiler.Const(CartesianIndex(-1, -1, -1), false))::CartesianIndex{3} │ (gci = %98 + cmin) │ (face_lid = -1) │ %101 = (a:b)::UnitRange{Int64} │ (@_35 = Base.iterate(%101)) │ %103 = (@_35 === nothing)::Bool │ %104 = Base.not_int(%103)::Bool └─── goto #14 if not %104 10 ┄ %106 = @_35::Tuple{Int64,Int64}::Tuple{Int64,Int64} │ (j@_40 = Core.getfield(%106, 1)) │ %108 = Core.getfield(%106, 2)::Int64 │ %109 = Base.getproperty(cell_to_faces, :data)::Array{Int64,1} │ %110 = Base.getindex(%109, j@_40)::Int64 │ %111 = (%110 == face_gid)::Bool └─── goto #12 if not %111 11 ─ %113 = (j@_40 - a)::Int64 │ (face_lid = %113 + 1) └─── goto #14 12 ─ (@_35 = Base.iterate(%101, %108)) │ %117 = (@_35 === nothing)::Bool │ %118 = Base.not_int(%117)::Bool └─── goto #14 if not %118 13 ─ goto #10 14 ┄ %121 = (face_lid::Int64 != -1)::Bool └─── goto #16 if not %121 15 ─ goto #17 16 ─ %124 = Base.AssertionError("face_lid != -1")::AssertionError └─── Base.throw(%124) 17 ┄ %126 = face_lid::Int64::Int64 │ %127 = offsets::Any │ %128 = (d@_29 + 1)::Int64 │ %129 = Base.getindex(%127, %128)::Any │ (face_lid = %126 + %129) │ (is_assigned_face_delta = Gridap.Geometry.isassigned(face_deltas, face_lid)) └─── goto #19 if not is_assigned_face_delta 18 ─ %133 = gci::CartesianIndex{3} │ %134 = Base.getindex(face_deltas, face_lid)::Any │ %135 = (%133 + %134)::Any │ (@_43 = %135 in gcis) └─── goto #20 19 ─ (@_43 = false) 20 ┄ goto #22 if not @_43 21 ─ Base.setindex!(face_to_geolabel, interior_id, face_gid) └─── goto #31 22 ─ (cell_found = false) │ %143 = (d@_29 + 1)::Int64 │ %144 = (D::Core.Compiler.Const(3, false) - 1)::Core.Compiler.Const(2, false) │ %145 = (%143:%144)::UnitRange{Int64} │ (@_39 = Base.iterate(%145)) │ %147 = (@_39 === nothing)::Bool │ %148 = Base.not_int(%147)::Bool └─── goto #28 if not %148 23 ┄ %150 = @_39::Tuple{Int64,Int64}::Tuple{Int64,Int64} │ (j@_42 = Core.getfield(%150, 1)) │ %152 = Core.getfield(%150, 2)::Int64 │ %153 = polytope_d_face_to_jfaces::Array{Array{Array{Int64,1},1},2} │ %154 = (d@_29 + 1)::Int64 │ %155 = (j@_42 + 1)::Int64 │ %156 = Base.getindex(%153, %154, %155)::Array{Array{Int64,1},1} │ %157 = face_lid::Any │ %158 = offsets::Any │ %159 = (d@_29 + 1)::Int64 │ %160 = Base.getindex(%158, %159)::Any │ %161 = (%157 - %160)::Any │ (dface_to_jfaces = Base.getindex(%156, %161)) │ %163 = dface_to_jfaces::Any │ %164 = offsets::Any │ %165 = (j@_42 + 1)::Int64 │ %166 = Base.getindex(%164, %165)::Any │ %167 = gcis::CartesianIndices{3,Tuple{Base.OneTo{Int64},Base.OneTo{Int64},Base.OneTo{Int64}}} │ %168 = gci::CartesianIndex{3} │ (cell_found = Gridap.Geometry._is_there_interior_cell_across_higher_dim_faces(%163, %166, %167, %168, face_deltas)) └─── goto #26 if not cell_found 24 ─ goto #28 25 ─ Core.Compiler.Const(:(goto %173), false) 26 ┄ (@_39 = Base.iterate(%145, %152)) │ %174 = (@_39 === nothing)::Bool │ %175 = Base.not_int(%174)::Bool └─── goto #28 if not %175 27 ─ goto #23 28 ┄ goto #30 if not cell_found 29 ─ Base.setindex!(face_to_geolabel, boundary_id::Core.Compiler.Const(-1, false), face_gid) └─── goto #31 30 ─ Base.setindex!(face_to_geolabel, face_lid, face_gid) 31 ┄ (@_28 = Base.iterate(%76, %86)) │ %183 = (@_28 === nothing)::Bool │ %184 = Base.not_int(%183)::Bool └─── goto #33 if not %184 32 ─ goto #9 33 ┄ (@_20 = Base.iterate(%62, %69)) │ %188 = (@_20 === nothing)::Bool │ %189 = Base.not_int(%188)::Bool └─── goto #35 if not %189 34 ─ goto #8 35 ┄ Gridap.Geometry._fix_geolabels(D::Core.Compiler.Const(3, false), topo, d_to_dface_to_entity, interior_id, boundary_id::Core.Compiler.Const(-1, false)) │ %193 = d_to_dface_to_entity::Array{Array{Int32,1},1} │ %194 = Base.lastindex(d_to_dface_to_entity)::Int64 │ %195 = Base.getindex(%193, %194)::Array{Int32,1} │ %196 = Gridap.Geometry.fill!(%195, interior_id)::Array{Int32,1} └─── return %196

Hopefully the function barrier is type stable. Is it possible to check this with the debugger?
Does this work? moving to just before calling the barrier, and call @call_warntype ... in the debugger prompt.

Great! This function looks green now! (except for the unions with nothing that are inevitable and unproblematic)

@fverdugo ... I already finished addressing your third round of comments.

But still the following line "%59 = Base.getindex(face_deltas, face_lid)::CartesianIndex" is shown with "CartesianIndex" in red. Any clue?

This means that the compiler was not able to determine the D parameter of CartesianIndex{D}

It looks like you are mixing different values of D in face_deltas

@fverdugo ... I already finished addressing your third round of comments.

There is still room for a 4th round 😄 if you want to fix the CartesianIndex in red!

There is still room for a 4th round smile if you want to fix the CartesianIndex in red!

Solved in 🤘 aa9e320

(d + 1)::Int64 │ %53 = Base.getindex(offsets, %52)::Int64 │ (face_lid = %51 + %53) │ %55 = gci::CartesianIndex{3} │ %56 = Base.getindex(face_deltas, face_lid)::CartesianIndex{3} │ %57 = (%55 + %56)::CartesianIndex{3} │ %58 = (%57 in gcis)::Bool └─── goto #12 if not %58

…ntities! introducing a function barrier (https://docs.julialang.org/en/v1/manual/performance-tips/index.html#kernel-functions-1)

Added a new, more general inner constructor to CartesianDiscreteModel

79d3403

amartinhuertas requested a review from fverdugo May 6, 2020 13:18

amartinhuertas mentioned this pull request May 6, 2020

Extend cartesian model face labels to distributedcase gridap/GridapDistributed.jl#10

Merged

fverdugo requested changes May 6, 2020

View reviewed changes

Restricted _find_ncube_face_neighbor_deltas helper function arguments to

e4eda3c

type ExtrusionPolytope{D}

fverdugo requested changes May 7, 2020

View reviewed changes

amartinhuertas added 2 commits May 7, 2020 20:24

Mitigated the effect of type instability in _fill_subgrid_cartesian_e…

a9e3d46

…ntities! introducing a function barrier (https://docs.julialang.org/en/v1/manual/performance-tips/index.html#kernel-functions-1)

Minor in _generate_subgrid_pregeo_label!

7fb3ae8

fverdugo approved these changes May 7, 2020

View reviewed changes

Type instability fix in _find_ncube_face_neighbor_deltas

aa9e320

fverdugo merged commit bb1f937 into master May 7, 2020

fverdugo deleted the adding_inner_constructor_to_cartesian_discrete_model branch May 7, 2020 12:47

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Added a new, more general inner constructor to CartesianDiscreteModel #245

Added a new, more general inner constructor to CartesianDiscreteModel #245

amartinhuertas commented May 6, 2020 •

edited

Loading

fverdugo May 6, 2020 •

edited

Loading

amartinhuertas May 7, 2020

codecov-io commented May 7, 2020 •

edited

Loading

fverdugo May 7, 2020

amartinhuertas May 7, 2020

fverdugo May 7, 2020 •

edited

Loading

amartinhuertas May 7, 2020

fverdugo May 7, 2020

fverdugo May 7, 2020

amartinhuertas May 7, 2020

fverdugo May 7, 2020

fverdugo May 7, 2020

amartinhuertas May 7, 2020

Added a new, more general inner constructor to CartesianDiscreteModel #245

Added a new, more general inner constructor to CartesianDiscreteModel #245

Conversation

amartinhuertas commented May 6, 2020 • edited Loading

fverdugo May 6, 2020 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

codecov-io commented May 7, 2020 • edited Loading

Codecov Report

Choose a reason for hiding this comment

Choose a reason for hiding this comment

fverdugo May 7, 2020 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

amartinhuertas commented May 6, 2020 •

edited

Loading

fverdugo May 6, 2020 •

edited

Loading

codecov-io commented May 7, 2020 •

edited

Loading

fverdugo May 7, 2020 •

edited

Loading