From 221af4903fce7cde689838b0349466f5d6df36ba Mon Sep 17 00:00:00 2001 From: Eric Neiva Date: Mon, 20 Nov 2023 12:09:56 +0100 Subject: [PATCH] Refactored way to compute signed distance --- src/AlgoimUtils/AlgoimUtils.jl | 55 ++++- src/Exports.jl | 1 + test/AlgoimUtilsTests/DualQuadraturesTests.jl | 2 +- .../TwoPhaseStokesAlgoimTests.jl | 215 +++++++++++++++++- 4 files changed, 262 insertions(+), 11 deletions(-) diff --git a/src/AlgoimUtils/AlgoimUtils.jl b/src/AlgoimUtils/AlgoimUtils.jl index c5831e6..3e8f64a 100644 --- a/src/AlgoimUtils/AlgoimUtils.jl +++ b/src/AlgoimUtils/AlgoimUtils.jl @@ -44,6 +44,7 @@ export init_bboxes export fill_cpp_data export fill_cpp_data_raw export compute_closest_point_projections +export compute_displacement export compute_normal_displacement export compute_distance_fe_function export delaunaytrian @@ -290,7 +291,6 @@ function compute_closest_point_projections(model::AdaptedDiscreteModel, limitstol::Float64=1.0e-8) reffe = ReferenceFE(lagrangian,Float64,order) rfespace = TestFESpace(model,reffe) - _φ = φ.φ _rφ = interpolate_everywhere(φ.φ,rfespace) rφ = AlgoimCallLevelSetFunction(_rφ,∇(_rφ)) cdesc = get_cartesian_descriptor(get_model(model)) @@ -395,9 +395,54 @@ function compute_normal_displacement( disps end -_dist(x,y) = begin - sign = norm(x) > norm(y) ? -1 : 1 - sign * √( (x[1]-y[1])^2 + (x[2]-y[2])^2 ) + +function compute_displacement( + cps::AbstractVector{<:Point}, + phi::AlgoimCallLevelSetFunction, + fun, + dt::Float64, + Ω::Triangulation) + # Note that cps must be (scalar) DoF-numbered, not lexicographic-numbered + searchmethod = KDTreeSearch() + cache1 = _point_to_cell_cache(searchmethod,Ω) + x_to_cell(x) = _point_to_cell!(cache1, x) + point_to_cell = lazy_map(x_to_cell, cps) + cell_to_points, _ = make_inverse_table(point_to_cell, num_cells(Ω)) + cell_to_xs = lazy_map(Broadcasting(Reindex(cps)), cell_to_points) + cell_point_xs = CellPoint(cell_to_xs, Ω, PhysicalDomain()) + fun_xs = evaluate(fun,cell_point_xs) + ∇φ_xs = evaluate(phi.∇φ,cell_point_xs) + cell_point_disp = lazy_map(Broadcasting(⋅),fun_xs,∇φ_xs) + cache_vals = array_cache(cell_point_disp) + cache_ctop = array_cache(cell_to_points) + disps = zeros(Float64,length(cps)) + for cell in 1:length(cell_to_points) + pts = getindex!(cache_ctop,cell_to_points,cell) + vals = getindex!(cache_vals,cell_point_disp,cell) + for (i,pt) in enumerate(pts) + val = vals[i] + disps[pt] = dt * val + end + end + disps +end + +signed_distance(φ::Function,x,y) = sign(φ(y))*norm(x-y) + +signed_distance(φ::T,x,y) where {T<:Number} = sign(φ)*norm(x-y) + +function _compute_signed_distance( + φ::AlgoimCallLevelSetFunction{<:Function,<:Function}, + cps::Vector{<:Point},cos::Matrix{<:Point}) + _dist(x,y) = signed_distance(φ.φ,x,y) + lazy_map(_dist,cps,cos) +end + +function _compute_signed_distance( + φ::AlgoimCallLevelSetFunction{<:CellField,<:CellField}, + cps::Vector{<:Point},cos::Matrix{<:Point}) + φs = get_free_dof_values(φ.φ) + lazy_map(signed_distance,φs,cps,cos) end function compute_distance_fe_function( @@ -411,7 +456,7 @@ function compute_distance_fe_function( rmodel = refine(bgmodel,order) cos = get_node_coordinates(rmodel) cos = node_to_dof_order(cos,fespace,rmodel,order) - dists = lazy_map(_dist,cps,cos) + dists = _compute_signed_distance(φ,cps,cos) FEFunction(fespace,dists) end diff --git a/src/Exports.jl b/src/Exports.jl index 0b22391..37a6fcc 100644 --- a/src/Exports.jl +++ b/src/Exports.jl @@ -58,6 +58,7 @@ end @publish AlgoimUtils fill_cpp_data @publish AlgoimUtils fill_cpp_data_raw @publish AlgoimUtils compute_closest_point_projections +@publish AlgoimUtils compute_displacement @publish AlgoimUtils compute_normal_displacement @publish AlgoimUtils compute_distance_fe_function @publish AlgoimUtils delaunaytrian diff --git a/test/AlgoimUtilsTests/DualQuadraturesTests.jl b/test/AlgoimUtilsTests/DualQuadraturesTests.jl index 6d5b61c..df0b269 100644 --- a/test/AlgoimUtilsTests/DualQuadraturesTests.jl +++ b/test/AlgoimUtilsTests/DualQuadraturesTests.jl @@ -1,4 +1,4 @@ -module AlgoimInterfaceTests +module DualQuadraturesTests using Test using Gridap diff --git a/test/AlgoimUtilsTests/TwoPhaseStokesAlgoimTests.jl b/test/AlgoimUtilsTests/TwoPhaseStokesAlgoimTests.jl index b0f4da8..2fff55f 100644 --- a/test/AlgoimUtilsTests/TwoPhaseStokesAlgoimTests.jl +++ b/test/AlgoimUtilsTests/TwoPhaseStokesAlgoimTests.jl @@ -8,7 +8,7 @@ using GridapPardiso using SparseMatricesCSR using Test -function main(n::Int,w::Float64,νˡ::Float64,νˢ::Float64,γ::Float64) +function steady_state(n::Int,w::Float64,νˡ::Float64,νˢ::Float64,γ::Float64) pmin = Point(-1.5,-1.0) pmax = Point( 4.5, 1.0) @@ -147,14 +147,219 @@ function main(n::Int,w::Float64,νˡ::Float64,νˢ::Float64,γ::Float64) end -n = 80 -w = 1.0e-1 +function dynamic(n::Int,w::Float64,νˡ::Float64,νˢ::Float64,γ::Float64,Δt₀::Float64) + + pmin = Point(-1.5,-2.5) + pmax = Point( 3.5, 2.5) + partition = (n,n) + + model = CartesianDiscreteModel(pmin,pmax,partition) + Ω = Triangulation(model) + dp = pmax - pmin + h = dp[2]/n + + order = 2 + degree = order == 1 ? 3 : 2*order + + R = 0.12 + φ(x) = 1.0 - ( ( x[1]*x[1] + x[2]*x[2] ) / R^2 ) + reffeᵠ = ReferenceFE(lagrangian,Float64,order+1) + Vbg = TestFESpace(Ω,reffeᵠ) + + # Buffer of active model and integration objects + buffer = Ref{Any}((Ωˡ=nothing,Ωˢ=nothing, + dΩˡ=nothing,dΩˢ=nothing, + dΓ=nothing,n_Γ=nothing, + aggsˡ=nothing,aggsˢ=nothing, + φ₋=nothing,cp₋=nothing,t=nothing)) + + function update_buffer!(t,dt,v₋₂) + + if buffer[].t == t + return true + else + + if buffer[].φ₋ === nothing + __φ₋ = AlgoimCallLevelSetFunction(φ,∇(φ)) + _φ₋ = compute_distance_fe_function(model,Vbg,__φ₋,order+1,cppdegree=3) + else + cp₋₂ = buffer[].cp₋ + φ₋₂ = buffer[].φ₋ + ϕ₋ = compute_displacement(cp₋₂,φ₋₂,v₋₂,dt,Ω) + ϕ₋ = get_free_dof_values(φ₋₂.φ) - ϕ₋ + _φ₋ = FEFunction(Vbg,ϕ₋) + end + + φ₋ = AlgoimCallLevelSetFunction(_φ₋,∇(_φ₋)) + cp₋ = compute_closest_point_projections(Vbg,φ₋,order+1,cppdegree=3) + + lquad = Quadrature(algoim,φ₋,degree,phase=IN) + Ωˡ,dΩˡ,cell_to_is_liquid = TriangulationAndMeasure(Ω,lquad) + + squad = Quadrature(algoim,φ₋,degree,phase=OUT) + Ωˢ,dΩˢ,cell_to_is_solid = TriangulationAndMeasure(Ω,squad) + + iquad = Quadrature(algoim,φ₋,degree) + _,dΓ,cell_to_is_cut = TriangulationAndMeasure(Ω,iquad) + n_Γ = normal(φ₋,Ω) # Exterior to liquid + + aggsˡ = aggregate(Ω,cell_to_is_liquid,cell_to_is_cut,IN) + aggsˢ = aggregate(Ω,cell_to_is_solid,cell_to_is_cut,OUT) + + buffer[] = (Ωˡ=Ωˡ,Ωˢ=Ωˢ,dΩˡ=dΩˡ,dΩˢ=dΩˢ, + dΓ=dΓ,n_Γ=n_Γ,aggsˡ=aggsˡ,aggsˢ=aggsˢ, + φ₋=φ₋,cp₋=cp₋,t=t) + return true + + end + + end + + reffeᵘ = ReferenceFE(lagrangian,VectorValue{2,Float64},order) + reffeˢ = ReferenceFE(lagrangian,VectorValue{2,Float64},order,space=:S) + reffeᵖ = ReferenceFE(lagrangian,Float64,order-1,space=:P) + + labels = get_face_labeling(model) + add_tag_from_tags!(labels,"inlet",[7]) + + function update_all!(t::Real,dt::Real,disp) + + update_buffer!(t,dt,disp) + + Ωˡ = buffer[].Ωˡ + Ωˢ = buffer[].Ωˢ + + dΩˡ = buffer[].dΩˡ + dΩˢ = buffer[].dΩˢ + + dΓ = buffer[].dΓ + n_Γ = buffer[].n_Γ + + aggsˡ = buffer[].aggsˡ + aggsˢ = buffer[].aggsˢ + + Vˡstd = TestFESpace(Ωˡ,reffeᵘ,dirichlet_tags=["inlet"]) + Vˡser = TestFESpace(Ωˡ,reffeˢ,conformity=:L2) + Qˡstd = TestFESpace(Ωˡ,reffeᵖ) + + Vˢstd = TestFESpace(Ωˢ,reffeᵘ) + Vˢser = TestFESpace(Ωˢ,reffeˢ,conformity=:L2) + Qˢstd = TestFESpace(Ωˢ,reffeᵖ) + + Vˡ = AgFEMSpace(Vˡstd,aggsˡ,Vˡser) + Qˡ = AgFEMSpace(Qˡstd,aggsˡ) + Vˢ = AgFEMSpace(Vˢstd,aggsˢ,Vˢser) + Qˢ = AgFEMSpace(Qˢstd,aggsˢ) + K = ConstantFESpace(model) + + uᵢ(x) = VectorValue(w*(1.0-x[2]*x[2]),0.0) + Uˡ = TrialFESpace(Vˡ,[uᵢ]) + Pˡ = TrialFESpace(Qˡ) + Uˢ = TrialFESpace(Vˢ) + Pˢ = TrialFESpace(Qˢ) + L = TrialFESpace(K) + + Y = MultiFieldFESpace([Vˡ,Qˡ,K,Vˢ,Qˢ,K]) + X = MultiFieldFESpace([Uˡ,Pˡ,L,Uˢ,Pˢ,L]) + + X,Y,dΩˡ,dΩˢ,Ωˡ,Ωˢ,dΓ,n_Γ + + end + + t₀ = 0.0 + Δt = Δt₀ + + X,Y,dΩˡ,dΩˢ,Ωˡ,Ωˢ,dΓ,n_Γ = update_all!(t₀,Δt,VectorValue(0.0,0.0)) + + wˡ = CellField(νˢ/(νˡ+νˢ),Ω) + wˢ = CellField(νˡ/(νˡ+νˢ),Ω) + νᵞ = 2*νˡ*νˢ/(νˡ+νˢ) + + jumpᵘ(uˡ,uˢ) = uˡ-uˢ + σ(ε,ν) = 2*ν*ε + τ(ε) = one(ε) + meanᵗ(uˡ,pˡ,νˡ,uˢ,pˢ,νˢ) = + wˡ*(σ(ε(uˡ),νˡ))-pˡ*(τ∘(ε(uˡ))) + + wˢ*(σ(ε(uˢ),νˢ))-pˢ*(τ∘(ε(uˢ))) + + a((uˡ,pˡ,lˡ,uˢ,pˢ,lˢ),(vˡ,qˡ,kˡ,vˢ,qˢ,kˢ)) = + ∫( 2*νˡ*(ε(uˡ)⊙ε(vˡ)) - qˡ*(∇⋅uˡ) - (∇⋅vˡ)*pˡ )dΩˡ + + ∫( 2*νˢ*(ε(uˢ)⊙ε(vˢ)) - qˢ*(∇⋅uˢ) - (∇⋅vˢ)*pˢ )dΩˢ + + ∫( pˡ*kˡ )dΩˡ + ∫( qˡ*lˡ )dΩˡ + + ∫( pˢ*kˢ )dΩˢ + ∫( qˢ*lˢ )dΩˢ + + ∫( (γ*νᵞ/h)*(jumpᵘ(uˡ,uˢ)⋅jumpᵘ(vˡ,vˢ)) - + (jumpᵘ(vˡ,vˢ)⋅(n_Γ⋅meanᵗ(uˡ,pˡ,νˡ,uˢ,pˢ,νˢ))) - + ((n_Γ⋅meanᵗ(vˡ,qˡ,νˡ,vˢ,qˢ,νˢ))⋅jumpᵘ(uˡ,uˢ)) )dΓ + + l((vˡ,qˡ,kˡ,vˢ,qˢ,kˢ)) = 0.0 + + assem = SparseMatrixAssembler(SymSparseMatrixCSR{1,Float64,Int},Vector{Float64},X,Y) + op = AffineFEOperator(a,l,X,Y,assem) + + A = get_matrix(op) + b = get_vector(op) + x = similar(b) + msglvl = 0 + # Customizing solver for real symmetric indefinite matrices + # See https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/sparse-solver-routines/onemkl-pardiso-parallel-direct-sparse-solver-iface/pardiso-iparm-parameter.html + iparm = new_iparm() + iparm[0+1] = 1 # Use default values (0 = enabled) + iparm[1+1] = 3 # Fill-in reducing ordering for the input matrix + iparm[9+1] = 8 # Pivoting perturbation + iparm[10+1] = 1 # Scaling vectors + iparm[12+1] = 1 # Improved accuracy using (non-) symmetric weighted matching + # Note that in the analysis phase (phase=11) you must provide the numerical values + # of the matrix A in array a in case of scaling and symmetric weighted matching. + iparm[17+1] = -1 # Report the number of non-zero elements in the factors. + iparm[18+1] = -1 # Report number of floating point operations (in 10^6 floating point operations) that are necessary to factor the matrix A. + iparm[20+1] = 1 # Pivoting for symmetric indefinite matrices. + ps = PardisoSolver(GridapPardiso.MTYPE_REAL_SYMMETRIC_INDEFINITE, iparm, msglvl) + ss = symbolic_setup(ps, A) + ns = numerical_setup(ss, A) + solve!(x, ns, b) + xh = FEFunction(X,x) + uhl, phl, _, uhs, phs, _ = xh + + _A = get_matrix(op) + _b = get_vector(op) + _x = get_free_dof_values(xh) + _r = _A*_x - _b + nr = norm(_r) + nb = norm(_b) + nx = norm(_x) + # @show nr, nr/nb, nr/nx + tol_warn = 1.0e-10 + if nr > tol_warn && nr/nb > tol_warn && nr/nx > tol_warn + @warn "Solver not accurate" + end + + # colors = color_aggregates(aggsˡ,model) + # writevtk(Ω,"res_bg_l",celldata=["aggregate"=>aggsˡ,"color"=>colors]) + # writevtk(Ω,"res_bg_s",celldata=["aggregate"=>aggsˢ,"color"=>colors]) + writevtk(Ωˡ,"res_l",cellfields=["uhl"=>uhl,"phl"=>phl]) + writevtk(dΩˢ,"res_s",cellfields=["uhs"=>uhs,"phs"=>phs]) + writevtk(dΓ,"res_gam",cellfields=["uhl"=>uhl,"uhs"=>uhs],qhulltype=convexhull) + # nh = interpolate_everywhere(n_Γ,Vˡstd) + # σn = νˡ*(ε(uhl)⋅nh) # - phl*nh + # writevtk(dΓ,"res_gam",cellfields=["uhl"=>uhl,"phl"=>phl,"sn"=>σn],qhulltype=convexhull) + + X,Y,dΩˡ,dΩˢ,Ωˡ,Ωˢ,dΓ,n_Γ = update_all!(t₀+Δt,Δt,uhs) + + writevtk(Ωˡ,"res_l_2") + writevtk(dΩˢ,"res_s_2") + writevtk(dΓ,"res_gam_2",qhulltype=convexhull) + +end + +n = 80 +w = 1.0e-1 νˡ = 1.0e-1 νˢ = 1.0e+2 -γ = 1.0e-1 # Scale with νˢ +γ = 1.0e-1 # Scale with νˢ +Δt = 1.0 @info "Values: n = $n, w = $w, νˡ = $νˡ, νˢ = $νˢ, γ = $γ" -main(n,w,νˡ,νˢ,γ) +dynamic(n,w,νˡ,νˢ,γ,Δt) end # module \ No newline at end of file