Skip to content

Commit 3c406c7

Browse files
termi-officialtermi-official
authored
Add Land et al. verification (#197)
* Fix pressure assembly bug * Fix nonlinear callback * Add first verification problem by Land et al. 2015. * Add weak bc consistency check * Merge solid tests * Fix 3D0D init * Add option to ignore monotonic convergence in Newton
1 parent 3271f71 commit 3c406c7

13 files changed

Lines changed: 490 additions & 253 deletions

File tree

‎src/discretization/operator.jl‎

Lines changed: 2 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -201,9 +201,7 @@ function _update_linearization_on_subdomain_J!(assembler, sdh, element_cache, fa
201201
@timeit_debug "assemble element" assemble_element!(Jₑ, uₑ, cell, element_cache, time)
202202
# TODO maybe it makes sense to merge this into the element routine in a modular fasion?
203203
# TODO benchmark against putting this into the FacetIterator
204-
@timeit_debug "assemble faces" for local_face_index ∈ 1:nfacets(cell)
205-
assemble_face!(Jₑ, uₑ, cell, local_face_index, face_cache, time)
206-
end
204+
@timeit_debug "assemble faces" assemble_element!(Jₑ, uₑ, cell, face_cache, time)
207205
@timeit_debug "assemble tying" assemble_tying!(Jₑ, uₑ, uₜ, cell, tying_cache, time)
208206
assemble!(assembler, celldofs(cell), Jₑ)
209207
end
@@ -242,9 +240,7 @@ function _update_linearization_on_subdomain_Jr!(assembler, sdh, element_cache, f
242240
@timeit_debug "assemble element" assemble_element!(Jₑ, rₑ, uₑ, cell, element_cache, time)
243241
# TODO maybe it makes sense to merge this into the element routine in a modular fasion?
244242
# TODO benchmark against putting this into the FacetIterator
245-
@timeit_debug "assemble faces" for local_face_index ∈ 1:nfacets(cell)
246-
assemble_face!(Jₑ, rₑ, uₑ, cell, local_face_index, face_cache, time)
247-
end
243+
@timeit_debug "assemble faces" assemble_element!(Jₑ, rₑ, uₑ, cell, face_cache, time)
248244
@timeit_debug "assemble tying" assemble_tying!(Jₑ, rₑ, uₑ, uₜ, cell, tying_cache, time)
249245
assemble!(assembler, dofs, Jₑ, rₑ)
250246
end

‎src/modeling/core/composite_elements.jl‎

Lines changed: 34 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -7,29 +7,29 @@ struct CompositeVolumetricElementCache{CacheTupleType <: Tuple} <: AbstractVolum
77
inner_caches::CacheTupleType
88
end
99
# Main entry point for bilinear operators
10-
assemble_element!(Kₑ::AbstractMatrix, cell::CellCache, element_cache::CompositeVolumetricElementCache, time) = assemble_element!(Kₑ, cell, element_cache.inner_caches, time)
11-
@unroll function assemble_element!(Kₑ::AbstractMatrix, cell::CellCache, inner_caches::CacheTupleType, time) where CacheTupleType <: Tuple
10+
assemble_element!(Kₑ::AbstractMatrix, cell::CellCache, element_cache::CompositeVolumetricElementCache, time) = assemble_composite_element!(Kₑ, cell, element_cache.inner_caches, time)
11+
@unroll function assemble_composite_element!(Kₑ::AbstractMatrix, cell::CellCache, inner_caches::CacheTupleType, time) where CacheTupleType <: Tuple
1212
@unroll for inner_cache ∈ inner_caches
1313
assemble_element!(Kₑ, cell, inner_cache, time)
1414
end
1515
end
1616
# Update element matrix in nonlinear operators
17-
assemble_element!(Kₑ::AbstractMatrix, uₑ::AbstractVector, cell::CellCache, element_cache::CompositeVolumetricElementCache, time) = assemble_element!(Kₑ, uₑ, cell, element_cache.inner_caches, time)
18-
@unroll function assemble_element!(Kₑ::AbstractMatrix, uₑ::AbstractVector, cell::CellCache, inner_caches::CacheTupleType, time) where CacheTupleType <: Tuple
17+
assemble_element!(Kₑ::AbstractMatrix, uₑ::AbstractVector, cell::CellCache, element_cache::CompositeVolumetricElementCache, time) = assemble_composite_element!(Kₑ, uₑ, cell, element_cache.inner_caches, time)
18+
@unroll function assemble_composite_element!(Kₑ::AbstractMatrix, uₑ::AbstractVector, cell::CellCache, inner_caches::CacheTupleType, time) where CacheTupleType <: Tuple
1919
@unroll for inner_cache ∈ inner_caches
2020
assemble_element!(Kₑ, uₑ, cell, inner_cache, time)
2121
end
2222
end
2323
# Update element matrix and residual in nonlinear operators
24-
assemble_element!(Kₑ::AbstractMatrix, residualₑ::AbstractVector, uₑ::AbstractVector, cell::CellCache, element_cache::CompositeVolumetricElementCache, time) = assemble_element!(Kₑ, residualₑ, uₑ, cell, element_cache.inner_caches, time)
25-
@unroll function assemble_element!(Kₑ::AbstractMatrix, residualₑ::AbstractVector, uₑ::AbstractVector, cell::CellCache, inner_caches::CacheTupleType, time) where CacheTupleType <: Tuple
24+
assemble_element!(Kₑ::AbstractMatrix, residualₑ::AbstractVector, uₑ::AbstractVector, cell::CellCache, element_cache::CompositeVolumetricElementCache, time) = assemble_composite_element!(Kₑ, residualₑ, uₑ, cell, element_cache.inner_caches, time)
25+
@unroll function assemble_composite_element!(Kₑ::AbstractMatrix, residualₑ::AbstractVector, uₑ::AbstractVector, cell::CellCache, inner_caches::CacheTupleType, time) where CacheTupleType <: Tuple
2626
@unroll for inner_cache ∈ inner_caches
2727
assemble_element!(Kₑ, residualₑ, uₑ, cell, inner_cache, time)
2828
end
2929
end
3030
# Update residual in nonlinear operators
31-
assemble_element!(residualₑ::AbstractVector, uₑ::AbstractVector, cell::CellCache, element_cache::CompositeVolumetricElementCache, time) = assemble_element!(residualₑ, uₑ, cell, element_cache.inner_caches, time)
32-
@unroll function assemble_element!(residualₑ::AbstractVector, uₑ::AbstractVector, cell::CellCache, inner_caches::CacheTupleType, time) where CacheTupleType <: Tuple
31+
assemble_element!(residualₑ::AbstractVector, uₑ::AbstractVector, cell::CellCache, element_cache::CompositeVolumetricElementCache, time) = assemble_composite_element!(residualₑ, uₑ, cell, element_cache.inner_caches, time)
32+
@unroll function assemble_composite_element!(residualₑ::AbstractVector, uₑ::AbstractVector, cell::CellCache, inner_caches::CacheTupleType, time) where CacheTupleType <: Tuple
3333
@unroll for inner_cache ∈ inner_caches
3434
assemble_element!(residualₑ, uₑ, cell, inner_cache, time)
3535
end
@@ -67,34 +67,41 @@ struct CompositeSurfaceElementCache{CacheTupleType <: Tuple} <: AbstractSurfaceE
6767
inner_caches::CacheTupleType
6868
end
6969
# Main entry point for bilinear operators
70-
assemble_face!(Kₑ::AbstractMatrix, cell::CellCache, local_facet_index::Int, surface_cache::CompositeSurfaceElementCache, time) = assemble_face!(Kₑ, cell, surface_cache.inner_caches, local_facet_index, time)
71-
@unroll function assemble_face!(Kₑ::AbstractMatrix, cell::CellCache, local_facet_index::Int, inner_caches::CacheTupleType, time) where CacheTupleType <: Tuple
70+
assemble_face!(Kₑ::AbstractMatrix, cell::CellCache, local_facet_index::Int, surface_cache::CompositeSurfaceElementCache, time) = assemble_composite_face!(Kₑ, cell, surface_cache.inner_caches, local_facet_index, time)
71+
@unroll function assemble_composite_face!(Kₑ::AbstractMatrix, cell::CellCache, local_facet_index::Int, inner_caches::CacheTupleType, time) where CacheTupleType <: Tuple
7272
@unroll for inner_cache ∈ inner_caches
7373
assemble_face!(Kₑ, cell, local_facet_index, inner_cache, time)
7474
end
7575
end
7676
# Update element matrix in nonlinear operators
77-
assemble_face!(Kₑ::AbstractMatrix, uₑ::AbstractVector, cell::CellCache, local_facet_index::Int, surface_cache::CompositeSurfaceElementCache, time) = assemble_face!(Kₑ, uₑ, cell, local_facet_index, surface_cache.inner_caches, time)
78-
@unroll function assemble_face!(Kₑ::AbstractMatrix, uₑ::AbstractVector, cell::CellCache, local_facet_index::Int, inner_caches::CacheTupleType, time) where CacheTupleType <: Tuple
77+
assemble_face!(Kₑ::AbstractMatrix, uₑ::AbstractVector, cell::CellCache, local_facet_index::Int, surface_cache::CompositeSurfaceElementCache, time) = assemble_composite_face!(Kₑ, uₑ, cell, local_facet_index, surface_cache.inner_caches, time)
78+
@unroll function assemble_composite_face!(Kₑ::AbstractMatrix, uₑ::AbstractVector, cell::CellCache, local_facet_index::Int, inner_caches::CacheTupleType, time) where CacheTupleType <: Tuple
7979
@unroll for inner_cache ∈ inner_caches
8080
assemble_face!(Kₑ, uₑ, cell, local_facet_index, inner_cache, time)
8181
end
8282
end
8383
# Update element matrix and residual in nonlinear operators
84-
assemble_face!(Kₑ::AbstractMatrix, residualₑ::AbstractVector, uₑ::AbstractVector, cell::CellCache, local_facet_index::Int, surface_cache::CompositeSurfaceElementCache, time) = assemble_face!(Kₑ, residualₑ, uₑ, cell, local_facet_index, surface_cache.inner_caches, time)
85-
@unroll function assemble_face!(Kₑ::AbstractMatrix, residualₑ::AbstractVector, uₑ::AbstractVector, cell::CellCache, local_facet_index::Int, inner_caches::CacheTupleType, time) where CacheTupleType <: Tuple
84+
assemble_face!(Kₑ::AbstractMatrix, residualₑ::AbstractVector, uₑ::AbstractVector, cell::CellCache, local_facet_index::Int, surface_cache::CompositeSurfaceElementCache, time) = assemble_composite_face!(Kₑ, residualₑ, uₑ, cell, local_facet_index, surface_cache.inner_caches, time)
85+
@unroll function assemble_composite_face!(Kₑ::AbstractMatrix, residualₑ::AbstractVector, uₑ::AbstractVector, cell::CellCache, local_facet_index::Int, inner_caches::CacheTupleType, time) where CacheTupleType <: Tuple
8686
@unroll for inner_cache ∈ inner_caches
8787
assemble_face!(Kₑ, residualₑ, uₑ, cell, local_facet_index, inner_cache, time)
8888
end
8989
end
9090
# Update residual in nonlinear operators
91-
assemble_face!(residualₑ::AbstractVector, uₑ::AbstractVector, cell::CellCache, local_facet_index::Int, surface_cache::CompositeSurfaceElementCache, time) = assemble_face!(residualₑ, uₑ, cell, local_facet_index, surface_cache.inner_caches, time)
92-
@unroll function assemble_face!(residualₑ::AbstractVector, uₑ::AbstractVector, cell::CellCache, local_facet_index::Int, inner_caches::CacheTupleType, time) where CacheTupleType <: Tuple
91+
assemble_face!(residualₑ::AbstractVector, uₑ::AbstractVector, cell::CellCache, local_facet_index::Int, surface_cache::CompositeSurfaceElementCache, time) = assemble_composite_face!(residualₑ, uₑ, cell, local_facet_index, surface_cache.inner_caches, time)
92+
@unroll function assemble_composite_face!(residualₑ::AbstractVector, uₑ::AbstractVector, cell::CellCache, local_facet_index::Int, inner_caches::CacheTupleType, time) where CacheTupleType <: Tuple
9393
@unroll for inner_cache ∈ inner_caches
9494
assemble_face!(residualₑ, uₑ, cell, local_facet_index, inner_cache, time)
9595
end
9696
end
9797

98+
# If we compose a face cache into an element cache, then we loop over the faces of the elements and try to assemble
99+
# Update element matrix in nonlinear operators
100+
assemble_element!(Kₑ::AbstractMatrix, uₑ::AbstractVector, cell::CellCache, surface_cache::CompositeSurfaceElementCache, time) = assemble_composite_element!(Kₑ, uₑ, cell, surface_cache.inner_caches, time)
101+
# Update element matrix and residual in nonlinear operators
102+
assemble_element!(Kₑ::AbstractMatrix, residualₑ::AbstractVector, uₑ::AbstractVector, cell::CellCache, surface_cache::CompositeSurfaceElementCache, time) = assemble_composite_element!(Kₑ, residualₑ, uₑ, cell, surface_cache.inner_caches, time)
103+
# Update residual in nonlinear operators
104+
assemble_element!(residualₑ::AbstractVector, uₑ::AbstractVector, cell::CellCache, surface_cache::CompositeSurfaceElementCache, time) = assemble_composite_element!(residualₑ, uₑ, cell, surface_cache.inner_caches, time)
98105

99106

100107
"""
@@ -104,30 +111,30 @@ struct CompositeInterfaceElementCache{CacheTupleType <: Tuple} <: AbstractInterf
104111
inner_caches::CacheTupleType
105112
end
106113
# Main entry point for bilinear operators
107-
assemble_interface!(Kₑ::AbstractMatrix, interface, interface_cache::CompositeInterfaceElementCache, time) = assemble_interface!(Kₑ, interface, interface_cache.inner_caches, time)
108-
@unroll function assemble_interface!(Kₑ::AbstractMatrix, interface, inner_caches::CacheTupleType, time) where CacheTupleType <: Tuple
114+
assemble_interface!(Kₑ::AbstractMatrix, interface, interface_cache::CompositeInterfaceElementCache, time) = assemble_composite_interface!(Kₑ, interface, interface_cache.inner_caches, time)
115+
@unroll function assemble_composite_interface!(Kₑ::AbstractMatrix, interface, inner_caches::CacheTupleType, time) where CacheTupleType <: Tuple
109116
@unroll for inner_cache ∈ inner_caches
110117
assemble_interface!(Kₑ, interface, inner_cache, time)
111118
end
112119
end
113120
# Update element matrix in nonlinear operators
114-
assemble_interface!(Kₑ::AbstractMatrix, uₑ::AbstractVector, interface, interface_cache::CompositeInterfaceElementCache, time) = assemble_interface!(Kₑ, uₑ, interface, interface_cache.inner_caches, time)
115-
@unroll function assemble_interface!(Kₑ::AbstractMatrix, uₑ::AbstractVector, interface, inner_caches::CacheTupleType, local_facet_index::Int, time) where CacheTupleType <: Tuple
121+
assemble_interface!(Kₑ::AbstractMatrix, uₑ::AbstractVector, interface, interface_cache::CompositeInterfaceElementCache, time) = assemble_composite_interface!(Kₑ, uₑ, interface, interface_cache.inner_caches, time)
122+
@unroll function assemble_composite_interface!(Kₑ::AbstractMatrix, uₑ::AbstractVector, interface, inner_caches::CacheTupleType, time) where CacheTupleType <: Tuple
116123
@unroll for inner_cache ∈ inner_caches
117-
assemble_interface!(Kₑ, uₑ, interface, inner_cache, local_facet_index, time)
124+
assemble_interface!(Kₑ, uₑ, interface, inner_cache, time)
118125
end
119126
end
120127
# Update element matrix and residual in nonlinear operators
121-
assemble_interface!(Kₑ::AbstractMatrix, residualₑ::AbstractVector, uₑ::AbstractVector, interface, interface_cache::CompositeInterfaceElementCache, time) = assemble_interface!(Kₑ, residualₑ, uₑ, interface, interface_cache.inner_caches, time)
122-
@unroll function assemble_interface!(Kₑ::AbstractMatrix, residualₑ::AbstractVector, uₑ::AbstractVector, interface, inner_caches::CacheTupleType, local_facet_index::Int, time) where CacheTupleType <: Tuple
128+
assemble_interface!(Kₑ::AbstractMatrix, residualₑ::AbstractVector, uₑ::AbstractVector, interface, interface_cache::CompositeInterfaceElementCache, time) = assemble_composite_interface!(Kₑ, residualₑ, uₑ, interface, interface_cache.inner_caches, time)
129+
@unroll function assemble_composite_interface!(Kₑ::AbstractMatrix, residualₑ::AbstractVector, uₑ::AbstractVector, interface, inner_caches::CacheTupleType, time) where CacheTupleType <: Tuple
123130
@unroll for inner_cache ∈ inner_caches
124-
assemble_interface!(Kₑ, residualₑ, uₑ, interface, inner_cache, local_facet_index, time)
131+
assemble_interface!(Kₑ, residualₑ, uₑ, interface, inner_cache, time)
125132
end
126133
end
127134
# Update residual in nonlinear operators
128-
assemble_interface!(residualₑ::AbstractVector, uₑ::AbstractVector, interface, interface_cache::CompositeInterfaceElementCache, time) = assemble_interface!(Kₑ, interface, interface_cache.inner_caches, time)
129-
@unroll function assemble_interface!(residualₑ::AbstractVector, uₑ::AbstractVector, interface, inner_caches::CacheTupleType, local_facet_index::Int, time) where CacheTupleType <: Tuple
135+
assemble_interface!(residualₑ::AbstractVector, uₑ::AbstractVector, interface, interface_cache::CompositeInterfaceElementCache, time) = assemble_composite_interface!(Kₑ, interface, interface_cache.inner_caches, time)
136+
@unroll function assemble_composite_interface!(residualₑ::AbstractVector, uₑ::AbstractVector, interface, inner_caches::CacheTupleType, time) where CacheTupleType <: Tuple
130137
@unroll for inner_cache ∈ inner_caches
131-
assemble_interface!(residualₑ, uₑ, interface, inner_cache, local_facet_index, time)
138+
assemble_interface!(residualₑ, uₑ, interface, inner_cache, time)
132139
end
133140
end

‎src/modeling/core/weak_boundary_conditions.jl‎

Lines changed: 70 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -345,10 +345,12 @@ function assemble_face_pressure_qp!(Kₑ::AbstractMatrix, residualₑ::AbstractV
345345
invF = inv(F)
346346
cofF = transpose(invF)
347347
J = det(F)
348-
neumann_term = p * J * cofF
348+
# @info qp, J, cofF ⋅ n₀
349+
neumann_term = p * J * cofF ⋅ n₀
350+
# neumann_term = p * n₀
349351
for i in 1:ndofs_face
350352
δuᵢ = shape_value(fv, qp, i)
351-
residualₑ[i] += neumann_term ⋅ n₀ ⋅ δuᵢ * dΓ
353+
residualₑ[i] += neumann_term ⋅ δuᵢ * dΓ
352354

353355
for j in 1:ndofs_face
354356
∇δuⱼ = shape_gradient(fv, qp, j)
@@ -372,11 +374,11 @@ function assemble_face_pressure_qp!(Kₑ::AbstractMatrix, uₑ::AbstractVector,
372374

373375
∇u = function_gradient(fv, qp, uₑ)
374376
F = one(∇u) + ∇u
375-
377+
376378
invF = inv(F)
377379
cofF = transpose(invF)
378380
J = det(F)
379-
# neumann_term = p * J * cofF
381+
# neumann_term = p * J * cofF ⋅ n₀
380382
for i in 1:ndofs_face
381383
δuᵢ = shape_value(fv, qp, i)
382384

@@ -406,10 +408,11 @@ function assemble_face_pressure_qp!(residualₑ::AbstractVector, uₑ::AbstractV
406408
invF = inv(F)
407409
cofF = transpose(invF)
408410
J = det(F)
409-
neumann_term = p * J * cofF
411+
neumann_term = p * J * cofF ⋅ n₀
412+
# neumann_term = p * n₀
410413
for i in 1:ndofs_face
411414
δuᵢ = shape_value(fv, qp, i)
412-
residualₑ[i] += neumann_term ⋅ n₀ ⋅ δuᵢ * dΓ
415+
residualₑ[i] += neumann_term ⋅ δuᵢ * dΓ
413416
end
414417
end
415418

@@ -485,3 +488,64 @@ function assemble_face!(residualₑ::AbstractVector, uₑ::AbstractVector, cell,
485488
assemble_face_pressure_qp!(residualₑ, uₑ, p, qp, fv)
486489
end
487490
end
491+
492+
493+
# We can use this to debug weak BCs for their consistency
494+
struct ConsistencyCheckWeakBoundaryCondition{BC} <: AbstractWeakBoundaryCondition
495+
bc::BC
496+
Δ::Float64
497+
end
498+
499+
struct ConsistencyCheckWeakBoundaryConditionCache{IC} <: AbstractSurfaceElementCache
500+
inner_cache::IC
501+
Kₑfd::Matrix{Float64}
502+
uₑfd::Vector{Float64}
503+
residualₑfd::Vector{Float64}
504+
residualₑref::Vector{Float64}
505+
Δ::Float64
506+
end
507+
@inline is_facet_in_cache(facet::FacetIndex, cell::CellCache, face_cache::ConsistencyCheckWeakBoundaryConditionCache) = is_facet_in_cache(facet, cell, face_cache.inner_cache)
508+
@inline getboundaryname(face_cache::ConsistencyCheckWeakBoundaryConditionCache) = getboundaryname(face_cache.inner_cache)
509+
@inline getboundaryname(check::ConsistencyCheckWeakBoundaryCondition) = getboundaryname(check.bc)
510+
511+
function setup_boundary_cache(ccc::ConsistencyCheckWeakBoundaryCondition, qr::FacetQuadratureRule, sdh::SubDofHandler)
512+
N = ndofs_per_cell(sdh)
513+
return ConsistencyCheckWeakBoundaryConditionCache(
514+
setup_boundary_cache(ccc.bc, qr, sdh),
515+
zeros(N, N),
516+
zeros(N),
517+
zeros(N),
518+
zeros(N),
519+
ccc.Δ,
520+
)
521+
end
522+
523+
function assemble_face!(Kₑ::AbstractMatrix, residualₑ::AbstractVector, uₑ::AbstractVector, cell, local_face_index::Int, cache::ConsistencyCheckWeakBoundaryConditionCache, time)
524+
(; Δ, inner_cache, Kₑfd, uₑfd, residualₑfd, residualₑref) = cache
525+
526+
# The incoming element matrix might be non-empty, so we need to start by storing the offset.
527+
Kₑfd .= Kₑ
528+
529+
# The actual assembly is happening here
530+
assemble_face!(Kₑ, residualₑ, uₑ, cell, local_face_index, inner_cache, time)
531+
532+
# Now we get a fresh reference state to pull the differences
533+
fill!(residualₑref, 0.0)
534+
assemble_face!(residualₑref, uₑ, cell, local_face_index, inner_cache, time)
535+
# Here we actually compute teh finite difference
536+
for i in 1:length(uₑfd)
537+
fill!(residualₑfd, 0.0)
538+
uₑfd .= uₑ
539+
uₑfd[i] += Δ
540+
assemble_face!(residualₑfd, uₑfd, cell, local_face_index, inner_cache, time)
541+
residualₑfd .-= residualₑref
542+
residualₑfd /= Δ
543+
Kₑfd[:,i] .+= residualₑfd
544+
end
545+
546+
# Finally we check for consistency
547+
if maximum(abs.(Kₑfd .- Kₑ)) > Δ
548+
@warn "Inconsistent element $(cellid(cell)) face $(local_face_index)! Jacobian difference: $(maximum(abs.(Kₑfd .- Kₑ)))"
549+
@info uₑ
550+
end
551+
end

0 commit comments

Comments
 (0)