From c3fcaa79391711e8e0c85c604e2a1fd24cc412c2 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Thu, 31 Oct 2024 09:41:45 +1100 Subject: [PATCH] Updated GMG drivers --- src/Applications/hunt.jl | 63 +++++++++++-------- src/Fixes.jl | 2 +- src/Main.jl | 50 ++++++++++----- src/Meshers/hunt_mesher.jl | 19 ++++++ src/Solvers/gmg.jl | 125 ++++++++++++++++++++++++++----------- src/parameters.jl | 3 +- src/utils.jl | 4 +- src/weakforms.jl | 3 +- test/dev/gmg.jl | 27 ++++++++ 9 files changed, 212 insertions(+), 84 deletions(-) create mode 100644 test/dev/gmg.jl diff --git a/src/Applications/hunt.jl b/src/Applications/hunt.jl index fc0b5ff..bcda236 100644 --- a/src/Applications/hunt.jl +++ b/src/Applications/hunt.jl @@ -54,6 +54,8 @@ function _hunt(; σw1=0.1, σw2=10.0, tw=0.0, + order = 2, + order_j = order, nsums = 10, vtk=true, title = "test", @@ -63,12 +65,14 @@ function _hunt(; jac_assemble = false, solve = true, solver = :julia, + formulation = :mhd, + rt_scaling = false, verbose = true, BL_adapted = true, kmap_x = 1, kmap_y = 1, ranks_per_level = nothing - ) +) info = Dict{Symbol,Any}() params = Dict{Symbol,Any}( @@ -105,12 +109,21 @@ function _hunt(; N = Ha^2/Re f̄ = (L/(ρ*u0^2))*VectorValue(f) B̄ = (1/B0)*VectorValue(B) - α = 1.0 - β = 1.0/Re - γ = N σ̄1 = σw1/σ σ̄2 = σw2/σ + if formulation == :cfd # Option 1 (CFD) + α = 1.0 + β = 1.0/Re + γ = N + elseif formulation == :mhd # Option 2 (MHD) is chosen in the experimental article + α = (1.0/N) + β = (1.0/Ha^2) + γ = 1.0 + else + error("Unknown formulation") + end + # DiscreteModel in terms of reduced quantities model = hunt_mesh(parts,params,nc,rank_partition,L,tw,Ha,kmap_x,kmap_y,BL_adapted,ranks_per_level) @@ -119,30 +132,28 @@ function _hunt(; writevtk(model,"hunt_model") end + params[:fluid] = Dict{Symbol,Any}( + :domain=>nothing, + :α=>α, + :β=>β, + :γ=>γ, + :f=>f̄, + :B=>B̄, + :ζ=>ζ, + ) + if tw > 0.0 σ_Ω = solid_conductivity(σ̄1,σ̄2,Ω,get_cell_gids(model),get_face_labeling(model)) params[:solid] = Dict(:domain=>"solid",:σ=>σ_Ω) - params[:fluid] = Dict( - :domain=>"fluid", - :α=>α, - :β=>β, - :γ=>γ, - :B=>B̄, - :f=>f̄, - :ζ=>ζ, - ) - else - params[:fluid] = Dict( - :domain=>nothing, - :α=>α, - :β=>β, - :γ=>γ, - :f=>f̄, - :B=>B̄, - :ζ=>ζ, - ) + params[:fluid][:domain] = "fluid" end + params[:fespaces] = Dict{Symbol,Any}( + :order_u => order, + :order_j => order_j, + :rt_scaling => rt_scaling ? 1.0/get_mesh_size(model) : nothing + ) + # Boundary conditions params[:bcs] = Dict( @@ -255,13 +266,13 @@ end function hunt_mesh( parts,params, nc::Tuple,np::Tuple,L::Real,tw::Real,Ha::Real,kmap_x::Number,kmap_y::Number,BL_adapted::Bool, - ranks_per_level) + ranks_per_level +) 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) params[:model] = model else # Multigrid - base_model = Meshers.hunt_generate_base_mesh(nc,L,tw,Ha,kmap_x,kmap_y,BL_adapted) - mh = Meshers.generate_mesh_hierarchy(parts,base_model,0,ranks_per_level) + 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/Fixes.jl b/src/Fixes.jl index 860ae4b..29539e6 100644 --- a/src/Fixes.jl +++ b/src/Fixes.jl @@ -8,4 +8,4 @@ function Geometry.get_polytopes(model::GridapDistributed.DistributedDiscreteModel) polys = map(get_polytopes,local_views(model)) return PartitionedArrays.getany(polys) -end +end \ No newline at end of file diff --git a/src/Main.jl b/src/Main.jl index 6dd1ab3..dc828d5 100644 --- a/src/Main.jl +++ b/src/Main.jl @@ -269,7 +269,7 @@ function _fe_space(::Val{:j},params) uses_mg = space_uses_multigrid(params[:solver])[3] model = uses_mg ? params[:multigrid][:mh] : params[:model] - phi = rt_scaling(model,params[:fespaces]) + phi = rt_scaling(params[:model],params[:fespaces]) reffe_j = ReferenceFE(raviart_thomas,Float64,k-1;basis_type=:jacobi,phi=phi) params[:fespaces][:reffe_j] = reffe_j @@ -356,10 +356,11 @@ const TriangulationTypes = Union{Gridap.Triangulation,GridapDistributed.Distribu _interior(model,domain::DiscreteModelTypes) = Interior(domain) _interior(model,domain::TriangulationTypes) = domain -_interior(model,domain::Nothing) = Triangulation(model) # This should be removed, but Gridap needs fixes +_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) @@ -375,18 +376,25 @@ function _setup_trians!(params) solid = params[:solid] Ωs = !isnothing(solid) ? _interior(params[:model],solid[:domain]) : nothing - if !uses_multigrid(params[:solver]) + #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] + #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 @@ -441,21 +449,29 @@ function _allocate_solution(op::TransientFEOperator,args...) nothing end -function _get_cell_size(t::Triangulation) +# Mesh 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) - map(m->m^(1/d),meas) + return collect(Float64, meas .^ (1/d)) end function _get_cell_size(t::GridapDistributed.DistributedTriangulation) map(_get_cell_size,local_views(t)) end -_get_cell_size(m::DiscreteModelTypes) = _get_cell_size(Triangulation(m)) - -_get_mesh_size(m::DiscreteModel) = minimum(_get_cell_size(m)) +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)) +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/Meshers/hunt_mesher.jl b/src/Meshers/hunt_mesher.jl index 0338524..5afee45 100644 --- a/src/Meshers/hunt_mesher.jl +++ b/src/Meshers/hunt_mesher.jl @@ -53,6 +53,10 @@ end function hunt_add_tags!(model,L::Real,tw::Real) labels = get_face_labeling(model) + hunt_add_tags!(labels,model,L,tw) +end + +function hunt_add_tags!(labels,model,L::Real,tw::Real) if tw > 0.0 ## add solid tags # When model is part of a distributed model we are using # that the maximum entity (9) is the same in all parts, @@ -128,3 +132,18 @@ function hunt_generate_base_mesh( hunt_add_tags!(model,L,tw) return model end + +function hunt_generate_mesh_hierarchy( + parts,np_per_level,nc,L,tw,Ha,kmap_x,kmap_y,BL_adapted +) + Lt = L+tw + _nc = (nc[1],nc[2],3) + domain = (-1.0,1.0,-1.0,1.0,0.0,0.1) + CartesianModelHierarchy( + parts,np_per_level,domain,_nc; + nrefs = (2,2,1), + isperiodic = (false,false,true), + map = hunt_stretch_map(Lt,Ha,kmap_x,kmap_y,BL_adapted), + add_labels! = labels -> hunt_add_tags!(labels,nothing,L,tw) # TODO: will not work for solid + ) +end diff --git a/src/Solvers/gmg.jl b/src/Solvers/gmg.jl index af97167..76b75d1 100644 --- a/src/Solvers/gmg.jl +++ b/src/Solvers/gmg.jl @@ -4,33 +4,59 @@ function get_block_solver(::Val{:gmg},params) gmg_solver(Val(vars),params) end -function gmg_solver(::Val{(1,3)},params) +function gmg_solver(::Val{(:u,:j)},params) mh = params[:multigrid][:mh] - trials = MultiFieldFESpace(map(s -> s, params[:multigrid][:trials][[1,3]])) - tests = MultiFieldFESpace(map(s -> s, params[:multigrid][:tests][[1,3]])) - tests_u, tests_j = params[:multigrid][:tests][1], params[:multigrid][:tests][3] + trials_u = params[:multigrid][:trials][:u] + tests_u = params[:multigrid][:tests][:u] + trials_j = params[:multigrid][:trials][:j] + tests_j = params[:multigrid][:tests][:j] + trials = MultiFieldFESpace([trials_u, trials_j]) + tests = MultiFieldFESpace([tests_u, tests_j]) nlevs = num_levels(mh) k = params[:fespaces][:k] qdegree = map(lev -> 2*k+1,1:nlevs) + is_nonlinear = params[:fluid][:convection] + + reffe_p = params[:fespaces][:reffe_p] + Πp = MultilevelTools.LocalProjectionMap(divergence,reffe_p,2*k) function jacobian_uj(model) - _, dΩf, α, β, γ, σf, f, B, ζ = retrieve_fluid_params(model,params,k) - _, dΩs, σs = retrieve_solid_params(model,params,k) - _, params_thin_wall, _, _ = retrieve_bcs_params(model,params,k) - Πp = params[:fespaces][:Πp] + Ωf = Triangulation(model) + dΩf = Measure(Ωf,2*k) + _, _, α, β, γ, σf, f, B, ζ = retrieve_fluid_params(model,params,k) + _, params_tw, _, _, params_Λ = retrieve_bcs_params(model,params,k) function a(x,dx,y) u, j = x du, dj = dx v_u, v_j = y r = a_mhd_u_u(du,v_u,β,dΩf) + a_mhd_u_j(dj,v_u,γ,B,dΩf) + a_mhd_j_u(du,v_j,σf,B,dΩf) + a_mhd_j_j(dj,v_j,dΩf) - r = r + dc_mhd_u_u(u,du,v_u,α,dΩf) - for p in params_thin_wall + if is_nonlinear + r = r + dc_mhd_u_u(u,du,v_u,α,dΩf) + end + for p in params_tw r = r + a_thin_wall_j_j(dj,v_j,p...) end - if !isnothing(dΩs) - r = r + a_solid_j_j(dj,v_j,dΩs) + if abs(ζ) > eps(typeof(ζ)) + r = r + a_al_u_u(du,v_u,ζ,Πp,dΩf) + a_al_j_j(dj,v_j,ζ,dΩf) + end + return r + end + return a + end + + function biform_uj(model) + Ωf = Triangulation(model) + dΩf = Measure(Ωf,2*k) + _, _, α, β, γ, σf, f, B, ζ = retrieve_fluid_params(model,params,k) + _, params_tw, _, _, params_Λ = retrieve_bcs_params(model,params,k) + function a(dx,y) + du, dj = dx + v_u, v_j = y + r = a_mhd_u_u(du,v_u,β,dΩf) + a_mhd_u_j(dj,v_u,γ,B,dΩf) + a_mhd_j_u(du,v_j,σf,B,dΩf) + a_mhd_j_j(dj,v_j,dΩf) + for p in params_tw + r = r + a_thin_wall_j_j(dj,v_j,p...) end if abs(ζ) > eps(typeof(ζ)) r = r + a_al_u_u(du,v_u,ζ,Πp,dΩf) + a_al_j_j(dj,v_j,ζ,dΩf) @@ -41,10 +67,14 @@ function gmg_solver(::Val{(1,3)},params) end function jacobian_u(model) - _, dΩf, α, β, γ, σf, f, B, ζ = retrieve_fluid_params(model,params,k) - Πp = params[:fespaces][:Πp] + _, _, α, β, γ, σf, f, B, ζ = retrieve_fluid_params(model,params,k) + Ωf = Triangulation(model) + dΩf = Measure(Ωf,2*k) function a(u,du,v_u) - r = a_mhd_u_u(du,v_u,β,dΩf) + dc_mhd_u_u(u,du,v_u,α,dΩf) + r = a_mhd_u_u(du,v_u,β,dΩf) + if is_nonlinear + dc_mhd_u_u(u,du,v_u,α,dΩf) + end if abs(ζ) > eps(typeof(ζ)) r = r + a_al_u_u(du,v_u,ζ,Πp,dΩf) end @@ -53,22 +83,36 @@ function gmg_solver(::Val{(1,3)},params) return a end - function projection_rhs(model) - _, dΩf, α, β, γ, σf, f, B, ζ = retrieve_fluid_params(model,params,k) - Πp = params[:fespaces][:Πp] - #l(du,v_u) = a_al_u_u(du,v_u,ζ,Πp,dΩf) - l((du,dj),(v_u,v_j)) = a_al_u_u(du,v_u,ζ,Πp,dΩf) + a_al_j_j(dj,v_j,ζ,dΩf) + function rhs(model) + _, _, α, β, γ, σf, f, B, ζ = retrieve_fluid_params(model,params,k) + Ωf = Triangulation(model) + dΩf = Measure(Ωf,2*k) + l((u,j),(v_u,v_j)) = a_al_u_u(u,v_u,ζ,Πp,dΩf) + a_al_j_j(j,v_j,ζ,dΩf) + return l + end + function rhs_u(model) + _, _, α, β, γ, σf, f, B, ζ = retrieve_fluid_params(model,params,k) + Ωf = Triangulation(model) + dΩf = Measure(Ωf,2*k) + l(u,v_u) = a_al_u_u(u,v_u,ζ,Πp,dΩf) return l end weakforms = map(mhl -> jacobian_uj(GridapSolvers.get_model(mhl)),mh) - smoothers = gmg_patch_smoothers(mh,tests,jacobian_uj) - #restrictions = setup_restriction_operators(tests,qdegree;mode=:residual,solver=LUSolver()) - #prolongations_u = gmg_patch_prolongations(tests_u,jacobian_u,projection_rhs) - #prolongations_j = setup_prolongation_operators(tests_j,qdegree) - #prolongations = MultiFieldTransferOperator(tests,[prolongations_u,prolongations_j]) - prolongations = gmg_patch_prolongations(tests,jacobian_uj,projection_rhs) - restrictions = gmg_patch_restrictions(tests,prolongations,projection_rhs,qdegree;solver=LUSolver()) + smoothers = gmg_patch_smoothers(mh,trials,jacobian_uj;is_nonlinear,w=-0.2) + + mf_prolongation = false + if mf_prolongation + prolongations_u = gmg_patch_prolongations(tests_u,jacobian_u,rhs_u;is_nonlinear) + prolongations_j = setup_prolongation_operators(tests_j,qdegree) + prolongations = MultiFieldTransferOperator(tests,[prolongations_u,prolongations_j];op_type=:prolongation) + restrictions_u = gmg_patch_restrictions(tests_u,prolongations_u,rhs_u,qdegree;solver=LUSolver()) + restrictions_j = setup_restriction_operators(tests_j,qdegree) + restrictions = MultiFieldTransferOperator(tests,[restrictions_u,restrictions_j];op_type=:restriction) + else + prolongations = gmg_patch_prolongations(tests,biform_uj,biform_uj;is_nonlinear=false) + restrictions = gmg_patch_restrictions(tests,prolongations,biform_uj,qdegree;solver=LUSolver()) + end return gmg_solver(mh,trials,tests,weakforms,restrictions,prolongations,smoothers) end @@ -84,35 +128,46 @@ function gmg_solver(mh,trials,tests,weakforms,restrictions,prolongations,smoothe pre_smoothers=smoothers, post_smoothers=smoothers, coarsest_solver=coarsest_solver, maxiter=3, verbose=i_am_main(ranks), mode=:preconditioner, is_nonlinear=true ) - gmg.log.depth += 4 + gmg.log.depth = 6 solver = FGMRESSolver(10,gmg;m_add=5,maxiter=30,rtol=1.0e-6,verbose=i_am_main(ranks),name="UJ Block - FGMRES+GMG") - solver.log.depth += 3 # For printing purposes + solver.log.depth = 4 return solver end -function gmg_patch_smoothers(mh,tests,weakform) +function gmg_patch_smoothers( + mh,tests,weakform; + niter = 5, + w = 0.2, + is_nonlinear = true, patch_decompositions = PatchDecomposition(mh) +) spaces = view(map(GridapSolvers.get_fe_space,tests),1:num_levels(tests)-1) patch_spaces = PatchFESpace(tests,patch_decompositions) smoothers = map(patch_decompositions,patch_spaces,spaces) do PD, Ph, Vh - psolver = PatchBasedLinearSolver(weakform(PD),Ph,Vh,is_nonlinear=true) - RichardsonSmoother(psolver,5,0.2) + psolver = PatchBasedLinearSolver(weakform(PD),Ph,Vh;is_nonlinear) + if w < 0 + solver = GMRESSolver(niter;Pr=psolver,maxiter=niter) + patch_smoother = RichardsonSmoother(solver,1,1.0) + else + patch_smoother = RichardsonSmoother(psolver,niter,w) + end end return smoothers end -function gmg_patch_prolongations(sh,lhs,rhs) +function gmg_patch_prolongations(sh,lhs,rhs;is_nonlinear=true) map(view(linear_indices(sh),1:num_levels(sh)-1)) do lev cparts = get_level_parts(sh,lev+1) if i_am_in(cparts) model = get_model_before_redist(sh,lev) - PD = PatchDecomposition(model) + PD = GridapSolvers.PatchBasedSmoothers.CoarsePatchDecomposition(model) lhs_i = lhs(PD) rhs_i = rhs(PD) else PD, lhs_i, rhs_i = nothing, nothing, nothing end - PatchProlongationOperator(lev,sh,PD,lhs_i,rhs_i;is_nonlinear=true) + # TODO: We should give it both tests and trials in the nonlinear case + PatchProlongationOperator(lev,sh,PD,lhs_i,rhs_i;is_nonlinear) end end diff --git a/src/parameters.jl b/src/parameters.jl index c3a9be2..31ca46d 100644 --- a/src/parameters.jl +++ b/src/parameters.jl @@ -222,7 +222,7 @@ uses_petsc(solver::Symbol) = uses_petsc(Val(solver)) uses_petsc(::Val{:julia}) = false uses_petsc(::Val{:petsc}) = true uses_petsc(::Val{:li2019}) = true -uses_petsc(::Val{:badia2024}) = true +uses_petsc(::Val{:badia2024}) = false uses_multigrid(solver::Dict) = any(space_uses_multigrid(solver)) space_uses_multigrid(solver::Dict) = space_uses_multigrid(Val(solver[:solver]),solver) @@ -340,6 +340,7 @@ function params_multigrid(params::Dict{Symbol,Any}) # Init some variables multigrid[:trials] = Dict{Symbol,Any}() multigrid[:tests] = Dict{Symbol,Any}() + multigrid[:variables] = getindex([:u,:p,:j,:φ],findall(space_uses_multigrid(solver))) return multigrid end diff --git a/src/utils.jl b/src/utils.jl index 051ff18..34a087b 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -48,7 +48,7 @@ end # Function that given a filter set a name to those elements that satisfy the # conditions of the filter -function add_entity!(model,in,name) +function add_entity!(model,is_in,name) labels = get_face_labeling(model) node_coordinates = get_node_coordinates(model) entity = num_entities(labels) + 1 @@ -56,7 +56,7 @@ function add_entity!(model,in,name) facets = get_face_nodes(model,d) for (i,facet) in enumerate(facets) coord = sum(node_coordinates[facet])/length(facet) - if in(coord) + if is_in(coord) labels.d_to_dface_to_entity[d+1][i] = entity end end diff --git a/src/weakforms.jl b/src/weakforms.jl index 596ede0..3a66f53 100644 --- a/src/weakforms.jl +++ b/src/weakforms.jl @@ -237,8 +237,7 @@ function retrieve_bcs_params(model,params,k) for i in 1:length(params[:bcs][:stabilization]) Λ = _skeleton(model,params[:bcs][:stabilization][i][:domain]) dΛ = Measure(Λ,2*k) - _h = _get_cell_size(Λ) - h = CellField(_h,Λ) + h = get_cell_size(Λ) μ = params[:bcs][:stabilization][i][:μ] push!(params_Λ,(μ,h,dΛ)) end diff --git a/test/dev/gmg.jl b/test/dev/gmg.jl new file mode 100644 index 0000000..cd8f5d5 --- /dev/null +++ b/test/dev/gmg.jl @@ -0,0 +1,27 @@ + +using Gridap, GridapDistributed, PartitionedArrays, GridapMHD +using SparseArrays + +solver = Dict( + :solver => :badia2024, + :matrix_type => SparseMatrixCSC{Float64,Int64}, + :vector_type => Vector{Float64}, + :petsc_options => "", + :block_solvers => [:gmg,:cg_jacobi,:cg_jacobi], +) + +hunt( + nc=(4,4), + L=1.0, + B=(0.,10.,0.), + ζ=10.0, + order = 2, + order_j = 3, + vtk=false, + title="hunt", + solver=solver, + ranks_per_level=fill((1,1,1),2), + formulation = :mhd, + rt_scaling = true, + BL_adapted = false, +)