From 840493caf53cb4f47bdbbd9b8f0672d667b88c52 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Mon, 2 Dec 2024 14:56:17 +1100 Subject: [PATCH] Minor refactor to accomodate HDiv-conforming velocity --- src/Applications/channel.jl | 2 +- src/Applications/hunt.jl | 4 +- src/GridapMHD.jl | 1 + src/Solvers/li2019.jl | 2 +- src/geometry.jl | 148 ++++++++++++++++++++++++++++++++++ src/gridap_extras.jl | 35 --------- src/main.jl | 49 +----------- src/parameters.jl | 5 +- src/weakforms.jl | 153 ++++++++++++++++++++++++------------ test/seq/hunt_tests.jl | 1 + 10 files changed, 262 insertions(+), 138 deletions(-) create mode 100644 src/geometry.jl diff --git a/src/Applications/channel.jl b/src/Applications/channel.jl index 8d83462..72a0ce1 100644 --- a/src/Applications/channel.jl +++ b/src/Applications/channel.jl @@ -160,7 +160,7 @@ function _channel(; # Fluid parameters params[:fluid] = Dict( - :domain=>nothing, + :domain => nothing, :α=>α, :β=>β, :γ=>γ, diff --git a/src/Applications/hunt.jl b/src/Applications/hunt.jl index c274715..4cc9bee 100644 --- a/src/Applications/hunt.jl +++ b/src/Applications/hunt.jl @@ -287,9 +287,9 @@ function hunt_mesh( ranks_per_level = nothing, adaptivity_method = 0 ) if isnothing(ranks_per_level) # Single grid - model = Meshers.hunt_generate_base_mesh(parts,np,nc,L,tw,Ha,kmap_x,kmap_y,BL_adapted,mesh_postpro) + model = Meshers.hunt_generate_base_mesh(parts,np,nc,L,tw,Ha,kmap_x,kmap_y,BL_adapted) else # Multigrid - mh = Meshers.hunt_generate_mesh_hierarchy(parts,ranks_per_level,nc,L,tw,Ha,kmap_x,kmap_y,BL_adapted,mesh_postpro) + mh = Meshers.hunt_generate_mesh_hierarchy(parts,ranks_per_level,nc,L,tw,Ha,kmap_x,kmap_y,BL_adapted) params[:multigrid] = Dict{Symbol,Any}( :mh => mh, :num_refs_coarse => 0, diff --git a/src/GridapMHD.jl b/src/GridapMHD.jl index 954c277..ab3bc42 100644 --- a/src/GridapMHD.jl +++ b/src/GridapMHD.jl @@ -40,6 +40,7 @@ include("Solvers/badia2024.jl") # Main driver include("gridap_extras.jl") include("utils.jl") +include("geometry.jl") include("parameters.jl") include("weakforms.jl") include("main.jl") diff --git a/src/Solvers/li2019.jl b/src/Solvers/li2019.jl index 8ff0d7a..584fcb1 100644 --- a/src/Solvers/li2019.jl +++ b/src/Solvers/li2019.jl @@ -16,7 +16,7 @@ function Li2019Solver(op::FEOperator,params) # Preconditioner model = params[:model] k = max(params[:fespaces][:order_u],params[:fespaces][:order_j]) - Ωf = _interior(model,params[:fluid][:domain]) + Ωf = params[:Ωf] dΩ = Measure(Ωf,2*k) a_Dj = li2019_Dj(dΩ,params) a_Fk = li2019_Fk(dΩ,params) diff --git a/src/geometry.jl b/src/geometry.jl new file mode 100644 index 0000000..105f3dd --- /dev/null +++ b/src/geometry.jl @@ -0,0 +1,148 @@ + +const DiscreteModelTypes = Union{Gridap.DiscreteModel,GridapDistributed.DistributedDiscreteModel} +const TriangulationTypes = Union{Gridap.Triangulation,GridapDistributed.DistributedTriangulation} + +function setup_geometry!(params) + params[:geometry] = Dict{Symbol,Any}( + :interiors => Dict{UInt64,Any}(), + :boundaries => Dict{UInt64,Any}(), + :skeletons => Dict{UInt64,Any}(), + :normals => Dict{UInt64,Any}(), + :measures => Dict{UInt64,Any}(), + :other => Dict{UInt64,Any}(), + ) + + Ω = interior(params,nothing) + params[:Ω] = Ω + + fluid = params[:fluid] + Ωf = interior(params,fluid[:domain]) + params[:Ωf] = Ωf + + solid = params[:solid] + Ωs = !isnothing(solid) ? interior(params,solid[:domain]) : nothing + params[:Ωs] = Ωs + + if uses_multigrid(params[:solver]) + params[:multigrid][:Ω] = params[:multigrid][:mh] + params[:multigrid][:Ωf] = params[:multigrid][:mh] + params[:multigrid][:Ωs] = params[:multigrid][:mh] + end +end + +domain_hash(domain::Union{Nothing,DiscreteModelTypes,TriangulationTypes}) = objectid(domain) +domain_hash(domain::Union{String,Vector{String}}) = hash(join(domain)) + +domain_hash(model,domain::Union{String,Vector{String}}) = hash(join(domain),objectid(model)) +domain_hash(model,domain) = domain_hash(domain) + +interior(params::Dict,domain) = interior(params,params[:model],domain) + +function interior(params::Dict,model,domain) + interiors = params[:geometry][:interiors] + key = domain_hash(model,domain) + if haskey(interiors,key) + trian = interiors[key] + else + trian = _interior(model,domain) + interiors[domain_hash(trian)] = trian + interiors[key] = trian + end + return trian +end + +boundary(params::Dict,domain) = boundary(params,params[:model],domain) + +function boundary(params::Dict,model,domain) + boundaries = params[:geometry][:boundaries] + key = domain_hash(model,domain) + if haskey(boundaries,key) + trian = boundaries[key] + else + trian = _boundary(model,domain) + boundaries[domain_hash(trian)] = trian + boundaries[key] = trian + end + return trian +end + +skeleton(params::Dict,domain) = skeleton(params,params[:model],domain) + +function skeleton(params::Dict,model,domain) + skeletons = params[:geometry][:skeletons] + key = domain_hash(model,domain) + if haskey(skeletons,key) + trian = skeletons[key] + else + trian = _skeleton(model,domain) + skeletons[domain_hash(trian)] = trian + skeletons[key] = trian + end + return trian +end + +function normal_vector(params::Dict,trian::TriangulationTypes) + normals = params[:geometry][:normals] + key = domain_hash(trian) + if haskey(normals,key) + n = normals[key] + else + n = get_normal_vector(trian) + normals[key] = n + end + return n +end + +function measure(params::Dict,trian::TriangulationTypes) + measures = params[:geometry][:measures] + key = domain_hash(trian) + if haskey(measures,key) + m = measures[key] + else + d = num_cell_dims(trian) + q = params[:fespaces][:quadratures][d] + m = Measure(trian,q) + measures[key] = m + end + return m +end + +_interior(model,domain::DiscreteModelTypes) = Triangulation(domain) +_interior(model,domain::TriangulationTypes) = domain +_interior(model,domain::Nothing) = Triangulation(model) +_interior(model,domain) = Triangulation(model,tags=domain) + +_boundary(model,domain::TriangulationTypes) = domain +_boundary(model,domain::Nothing) = Boundary(model) +_boundary(model,domain) = Boundary(model,tags=domain) + +_skeleton(model,domain::TriangulationTypes) = Skeleton(domain) +_skeleton(model,domain::Nothing) = Skeleton(model) +_skeleton(model,domain) = Skeleton(_interior(model,domain)) + +# Face/cell sizes + +get_cell_size(t::TriangulationTypes) = CellField(_get_cell_size(t),t) + +get_cell_size(m::DiscreteModelTypes) = get_cell_size(Triangulation(m)) +_get_cell_size(m::DiscreteModelTypes) = _get_cell_size(Triangulation(m)) + +function _get_cell_size(t::Triangulation) :: Vector{Float64} + if iszero(num_cells(t)) + return Float64[] + end + meas = get_cell_measure(t) + d = num_dims(t) + return collect(Float64, meas .^ (1/d)) +end + +function _get_cell_size(t::GridapDistributed.DistributedTriangulation) + map(_get_cell_size,local_views(t)) +end + +get_mesh_size(m::DiscreteModel) = minimum(_get_cell_size(m)) + +function get_mesh_size(m::GridapDistributed.DistributedDiscreteModel) + h = map(get_mesh_size,local_views(m)) + return reduce(min,h;init=one(eltype(h))) +end diff --git a/src/gridap_extras.jl b/src/gridap_extras.jl index 8c999c2..5df8bc9 100644 --- a/src/gridap_extras.jl +++ b/src/gridap_extras.jl @@ -1,9 +1,4 @@ -# Quality of life definitions - -const DiscreteModelTypes = Union{Gridap.DiscreteModel,GridapDistributed.DistributedDiscreteModel} -const TriangulationTypes = Union{Gridap.Triangulation,GridapDistributed.DistributedTriangulation} - # Get polytopes function Geometry.get_polytopes(model::GridapDistributed.DistributedDiscreteModel) @@ -21,33 +16,6 @@ function Geometry.get_polytopes(trian::GridapDistributed.DistributedTriangulatio return getany(polys) end -# Face/cell sizes - -get_cell_size(t::TriangulationTypes) = CellField(_get_cell_size(t),t) - -get_cell_size(m::DiscreteModelTypes) = get_cell_size(Triangulation(m)) -_get_cell_size(m::DiscreteModelTypes) = _get_cell_size(Triangulation(m)) - -function _get_cell_size(t::Triangulation) :: Vector{Float64} - if iszero(num_cells(t)) - return Float64[] - end - meas = get_cell_measure(t) - d = num_dims(t) - return collect(Float64, meas .^ (1/d)) -end - -function _get_cell_size(t::GridapDistributed.DistributedTriangulation) - map(_get_cell_size,local_views(t)) -end - -get_mesh_size(m::DiscreteModel) = minimum(_get_cell_size(m)) - -function get_mesh_size(m::GridapDistributed.DistributedDiscreteModel) - h = map(get_mesh_size,local_views(m)) - return reduce(min,h;init=one(eltype(h))) -end - # MacroReferenceFEs function conformity_from_symbol(sym::Symbol) @@ -91,6 +59,3 @@ function Gridap.Adaptivity.MacroReferenceFE( reffes = Fill(reffe,Gridap.Adaptivity.num_subcells(rrule)) return Gridap.Adaptivity.MacroReferenceFE(rrule,reffes;macro_kwargs...) end - -# Local projection operators - diff --git a/src/main.jl b/src/main.jl index c161e14..11de219 100644 --- a/src/main.jl +++ b/src/main.jl @@ -124,8 +124,7 @@ function main(_params::Dict;output::Dict=Dict{Symbol,Any}()) params = add_default_params(_params) t = params[:ptimer] - # Compute triangulations only once for performance - _setup_trians!(params) + setup_geometry!(params) # FESpaces tic!(t;barrier=true) @@ -256,7 +255,6 @@ function _fe_space(::Val{:p},params) end function _fe_space(::Val{:j},params) - k = params[:fespaces][:order_j] uses_mg = space_uses_multigrid(params[:solver])[3] model = uses_mg ? params[:multigrid][:mh] : params[:model] @@ -355,51 +353,6 @@ function _continuation_fe_operator(mfs,U,V,params) return op end -# Sub-triangulations - -_interior(model,domain::DiscreteModelTypes) = Interior(domain) -_interior(model,domain::TriangulationTypes) = domain -_interior(model,domain::Nothing) = Triangulation(model) -_interior(model,domain) = Interior(model,tags=domain) - -_boundary(model,domain::TriangulationTypes) = domain -_boundary(model,domain::Nothing) = Boundary(model) -_boundary(model,domain) = Boundary(model,tags=domain) - -_skeleton(model,domain::TriangulationTypes) = SkeletonTriangulation(domain) -_skeleton(model,domain::Nothing) = SkeletonTriangulation(model) -_skeleton(model,domain) = _skeleton(model,_interior(model,domain)) - -function _setup_trians!(params) - Ω = _interior(params[:model],nothing) - - fluid = params[:fluid] - Ωf = _interior(params[:model],fluid[:domain]) - - solid = params[:solid] - Ωs = !isnothing(solid) ? _interior(params[:model],solid[:domain]) : nothing - - #if !uses_multigrid(params[:solver]) - params[:Ω] = Ω - params[:Ωf] = Ωf - params[:Ωs] = Ωs - #else - # params[:multigrid][:Ω] = Ω - # params[:multigrid][:Ωf] = Ωf - # params[:multigrid][:Ωs] = Ωs - # params[:Ω] = Ω[1] - # params[:Ωf] = Ωf[1] - # params[:Ωs] = Ωs[1] - #end - - if uses_multigrid(params[:solver]) - params[:multigrid][:Ω] = params[:multigrid][:mh] - params[:multigrid][:Ωf] = params[:multigrid][:mh] - params[:multigrid][:Ωs] = params[:multigrid][:mh] - end - -end - # Random vector generation function _rand(vt::Type{<:Vector{T}},r::AbstractUnitRange) where T diff --git a/src/parameters.jl b/src/parameters.jl index 66cc074..a4ba2d1 100644 --- a/src/parameters.jl +++ b/src/parameters.jl @@ -372,6 +372,8 @@ const FLUID_DISCRETIZATIONS = (; ) ) +has_hdiv_fluid_disc(params) = params[:fespaces][:fluid_disc] ∈ (:RT, :BDM) + function fluid_discretization(disc::Symbol,poly::Polytope,feparams) D = num_dims(poly) k = feparams[:order_u] @@ -509,11 +511,10 @@ function rt_scaling(poly,feparams) end function generate_quadratures(poly::Polytope{D},feparams) where D - fluid_disc = feparams[:fluid_disc] qdegree = feparams[:q] fpolys = ReferenceFEs.get_reffaces(poly)[2:end] # Skip 0-dfaces - if fluid_disc == :SV + if haskey(feparams,:rrule) # TODO: This should be made more general, to ensure all d-faces are properly # refined if need be. In the case of SV, there is not subdivision of the # facets and edges, so its fine. diff --git a/src/weakforms.jl b/src/weakforms.jl index 2d8a815..e4b0090 100644 --- a/src/weakforms.jl +++ b/src/weakforms.jl @@ -17,12 +17,10 @@ function _weak_form(params,k) dΩ = Measure(Ω,2*k) fluid = params[:fluid] - Ωf, dΩf, α, β, γ, σf, f, B, ζ, g = retrieve_fluid_params(params,k) - - solid = params[:solid] - Ωs, dΩs, σs = retrieve_solid_params(params,k) - - bcs_params = retrieve_bcs_params(params,k) + dΩf, α, β, γ, σf, f, B, ζ, g = retrieve_fluid_params(params) + dΩs, σs = retrieve_solid_params(params) + bcs_params = retrieve_bcs_params(params) + hdiv_params = retrieve_hdiv_fluid_params(params) params_φ, params_thin_wall, params_f, params_B, params_Λ = bcs_params Πp = local_projection_operator(params,k) @@ -41,6 +39,9 @@ function _weak_form(params,k) if has_solid(params) r = r + a_solid(x,dy,σs,dΩs) end + if !isnothing(hdiv_params) + r = r + a_HDiv(x,dy,hdiv_params...) + end if abs(ζ) > eps(typeof(ζ)) r = r + a_al(x,dy,ζ,Πp,dΩf,dΩ) end @@ -58,6 +59,9 @@ function _weak_form(params,k) for p in params_f r = r + ℓ_f(dy,p...) end + if !isnothing(hdiv_params) + r = r + ℓ_HDiv(dy,hdiv_params...) + end r end @@ -96,12 +100,10 @@ end function _ode_weak_form(params,k) fluid = params[:fluid] - Ωf, dΩf, α, β, γ, σf, f, B, ζ, g = retrieve_fluid_params(params,k) - - solid = params[:solid] - Ωs, dΩs, σs = retrieve_solid_params(params,k) - - bcs_params = retrieve_bcs_params(params,k) + dΩf, α, β, γ, σf, f, B, ζ, g = retrieve_fluid_params(params) + dΩs, σs = retrieve_solid_params(params) + bcs_params = retrieve_bcs_params(params) + hdiv_params = retrieve_hdiv_fluid_params(params) params_φ, params_thin_wall, params_f, params_B, params_Λ = bcs_params Πp = local_projection_operator(params,k) @@ -124,6 +126,9 @@ function _ode_weak_form(params,k) if has_solid(params) r = r + a_solid(x,dy,σs,dΩs) end + if !isnothing(hdiv_params) + r = r + a_HDiv(x,dy,hdiv_params...) + end if abs(ζ) > eps(typeof(ζ)) r = r + a_al(x,dy,ζ,Πp,dΩf,dΩ) end @@ -141,6 +146,9 @@ function _ode_weak_form(params,k) for p in params_f r = r + ℓ_f(dy,time_eval(p,t)...) end + if !isnothing(hdiv_params) + r = r + ℓ_HDiv(x,dy,hdiv_params...) + end r end @@ -168,50 +176,72 @@ end ############################################################################################ # Parameter retrieval -retrieve_fluid_params(params,k) = retrieve_fluid_params(params[:model],params,k) +retrieve_fluid_params(params) = retrieve_fluid_params(params[:model],params) -function retrieve_fluid_params(model,params,k) +function retrieve_fluid_params(model,params) fluid = params[:fluid] - quad3D = params[:fespaces][:quadratures][3] Ωf = params[:Ωf] - dΩf = Measure(Ωf,quad3D) + dΩf = measure(params,Ωf) α, β, γ, σf = fluid[:α], fluid[:β], fluid[:γ], fluid[:σ] - f = fluid[:f] - B = fluid[:B] - ζ = fluid[:ζ] - g = fluid[:g] - return Ωf, dΩf, α, β, γ, σf, f, B, ζ, g + f, B, ζ, g = fluid[:f], fluid[:B], fluid[:ζ], fluid[:g] + return dΩf, α, β, γ, σf, f, B, ζ, g +end + +retrieve_hdiv_fluid_params(params) = retrieve_hdiv_fluid_params(params[:model],params) + +function retrieve_hdiv_fluid_params(model,params) + Ωf = params[:Ωf] + + if has_hdiv_fluid_disc(params) + Γ = boundary(params,Ωf,nothing) + Λ = skeleton(params,Ωf,nothing) + Γ_D = boundary(params,Ωf,params[:bcs][:u][:tags]) + + h_Γ = get_cell_size(Γ) + h_Λ = get_cell_size(Λ) + n_Γ_D = normal_vector(params,Γ_D) + n_Λ = normal_vector(params,Λ) + + dΓ = measure(params,Γ) + dΓ_D = measure(params,Γ_D) + dΛ = measure(params,Λ) + + μ = 1.0 + u_D = params[:bcs][:u][:values] + + return (μ,h_Γ,h_Λ,n_Γ_D,n_Λ,u_D,dΓ,dΓ_D,dΛ) + else + return nothing + end + return hdiv_params end -retrieve_solid_params(params,k) = retrieve_solid_params(params[:model],params,k) +retrieve_solid_params(params) = retrieve_solid_params(params[:model],params) -function retrieve_solid_params(model,params,k) - quad3D = params[:fespaces][:quadratures][3] +function retrieve_solid_params(model,params) solid = params[:solid] if solid !== nothing Ωs = params[:Ωs] - dΩs = Measure(Ωs,quad3D) + dΩs = measure(params,Ωs) σs = solid[:σ] - return Ωs, dΩs, σs + return dΩs, σs else - return nothing, nothing, nothing + return nothing, nothing end end -retrieve_bcs_params(params,k) = retrieve_bcs_params(params[:model],params,k) +retrieve_bcs_params(params) = retrieve_bcs_params(params[:model],params) -function retrieve_bcs_params(model,params,k) - quad3D = params[:fespaces][:quadratures][3] - quad2D = params[:fespaces][:quadratures][2] +function retrieve_bcs_params(model,params) bcs = params[:bcs] params_φ = [] for i in 1:length(bcs[:φ]) φ_i = bcs[:φ][i][:value] - Γ = _boundary(model,bcs[:φ][i][:domain]) - dΓ = Measure(Γ,quad2D) - n_Γ = get_normal_vector(Γ) + Γ = boundary(params,bcs[:φ][i][:domain]) + dΓ = measure(params,Γ) + n_Γ = normal_vector(params,Γ) push!(params_φ,(φ_i,n_Γ,dΓ)) end @@ -220,9 +250,9 @@ function retrieve_bcs_params(model,params,k) τ_i = bcs[:thin_wall][i][:τ] cw_i = bcs[:thin_wall][i][:cw] jw_i = bcs[:thin_wall][i][:jw] - Γ = _boundary(model,bcs[:thin_wall][i][:domain]) - dΓ = Measure(Γ,quad2D) - n_Γ = get_normal_vector(Γ) + Γ = boundary(params,bcs[:thin_wall][i][:domain]) + dΓ = measure(params,Γ) + n_Γ = normal_vector(params,Γ) push!(params_thin_wall,(τ_i,cw_i,jw_i,n_Γ,dΓ)) end @@ -233,23 +263,23 @@ function retrieve_bcs_params(model,params,k) params_f = [] for i in 1:length(bcs[:f]) f_i = bcs[:f][i][:value] - Ω_i = _interior(model,bcs[:f][i][:domain]) - dΩ_i = Measure(Ω_i,quad3D) + Ω_i = interior(params,bcs[:f][i][:domain]) + dΩ_i = measure(params,Ω_i) push!(params_f,(f_i,dΩ_i)) end params_B = [] for i in 1:length(bcs[:B]) B_i = bcs[:B][i][:value] - Ω_i = _interior(model,bcs[:B][i][:domain]) - dΩ_i = Measure(Ω_i,quad3D) + Ω_i = interior(params,bcs[:B][i][:domain]) + dΩ_i = measure(params,Ω_i) push!(params_f,(γ,B_i,dΩ_i)) end params_Λ = [] for i in 1:length(params[:bcs][:stabilization]) - Λ = _skeleton(model,params[:bcs][:stabilization][i][:domain]) - dΛ = Measure(Λ,quad2D) + Λ = skeleton(params,params[:bcs][:stabilization][i][:domain]) + dΛ = measure(params,Λ) h = get_cell_size(Λ) μ = params[:bcs][:stabilization][i][:μ] push!(params_Λ,(μ,h,dΛ)) @@ -286,12 +316,6 @@ a_mhd_j_j(j,v_j,dΩ) = ∫( j⋅v_j )*dΩ a_mhd_j_φ(φ,v_j,σ,dΩ) = ∫( -σ*φ*(∇⋅v_j) )*dΩ a_mhd_φ_j(j,v_φ,dΩ) = ∫( -(∇⋅j)*v_φ )*dΩ -function a_Λ(x,dy,μ,h,dΛ) - u, p, j, φ = x - v_u, v_p, v_j, v_φ = dy - ∫( (1/2) * μ * (h*h) * jump( ∇(u) ) ⊙ jump( ∇(v_u) ))*dΛ -end - function ℓ_mhd(dy,f,dΩ) v_u, v_p, v_j, v_φ = dy ℓ_mhd_u(v_u,f,dΩ) @@ -327,6 +351,14 @@ function p_dc_mhd(x,dx,dy,α,dΩ) end p_dc_mhd_u_u(u,du,v_u,α,dΩ) = ∫( α*v_u⋅( conv∘(u,∇(du)) ) ) * dΩ +# Stabilisation + +function a_Λ(x,dy,μ,h,dΛ) + u, p, j, φ = x + v_u, v_p, v_j, v_φ = dy + ∫( (1/2) * μ * (h*h) * jump( ∇(u) ) ⊙ jump( ∇(v_u) ))*dΛ +end + # Augmented lagrangian function a_al(x,dy,ζ,Πp,dΩf,dΩ) @@ -353,6 +385,29 @@ function local_projection_operator(params,k) return Πp end +# Div-Conforming laplacian terms (parts integration + tangent component penalty) + +function a_HDiv(x,y,μ,h_Γ,h_Λ,n_Γ_D,n_Λ,u_D,dΓ,dΓ_D,dΛ) + u, _, _, _ = x + v, _, _, _ = y + + αΛ, αΓ = μ/h_Λ, μ/h_Γ + ∇u, ∇v = ∇(u), ∇(v) + uᵗ, vᵗ = jump(u⊗n_Λ), jump(v⊗n_Λ) + + c = ∫( αΛ*vᵗ⊙uᵗ - vᵗ⊙mean(∇u) - mean(∇v)⊙uᵗ)dΛ + c -= ∫(v⋅(∇u⋅n_Γ_D) + (∇v⋅n_Γ_D)⋅u)dΓ_D + c += ∫(αΓ*v⋅u)dΓ + return c +end + +function ℓ_HDiv(dy,μ,h_Γ,h_Λ,n_Γ_D,n_Λ,u_D,dΓ,dΓ_D,dΛ) + v, _, _, _ = dy + αΓ = μ/h_Γ + c = ∫(αΓ*v⋅u_D)dΓ - ∫((∇(v)⋅n_Γ_D)⋅u_D)dΓ_D + return c +end + # Solid equations function a_solid(x,dy,σ,dΩ) diff --git a/test/seq/hunt_tests.jl b/test/seq/hunt_tests.jl index 3f9cf1d..00df2e1 100644 --- a/test/seq/hunt_tests.jl +++ b/test/seq/hunt_tests.jl @@ -21,6 +21,7 @@ hunt( vtk=true, title="hunt", solver=:julia, + fluid_disc = :RT, ) hunt(