From d83bf0ed2ec545072b6e2fba70facf6404c1a7a5 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Thu, 17 Oct 2024 23:10:26 +1100 Subject: [PATCH 1/3] Minor --- src/Solvers/petsc.jl | 6 +++++- test/dev/bugfix.jl | 51 +++++++++++++++++++++++++++++++++++--------- 2 files changed, 46 insertions(+), 11 deletions(-) diff --git a/src/Solvers/petsc.jl b/src/Solvers/petsc.jl index c6b4fb1..827141b 100644 --- a/src/Solvers/petsc.jl +++ b/src/Solvers/petsc.jl @@ -6,6 +6,10 @@ get_block_solver(::Val{:petsc_gmres_schwarz},params) = PETScLinearSolver(petsc_g get_block_solver(::Val{:petsc_gmres_amg},params) = PETScLinearSolver(petsc_gmres_amg_setup) get_block_solver(::Val{:petsc_from_options},params) = PETScLinearSolver() +# MUMPS + +MUMPS_MAX_MEM_INCREASE = 40 # TODO: Create preferences? + function petsc_mumps_setup(ksp) pc = Ref{GridapPETSc.PETSC.PC}() mumpsmat = Ref{GridapPETSc.PETSC.Mat}() @@ -17,9 +21,9 @@ function petsc_mumps_setup(ksp) @check_error_code GridapPETSc.PETSC.PCFactorGetMatrix(pc[],mumpsmat) # @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[], 4, 1) @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[], 7, 0) + @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[], 14, MUMPS_MAX_MEM_INCREASE) @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[], 28, 2) @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[], 29, 1) - # @check_error_code GridapPETSc.PETSC.MatMumpsSetCntl(mumpsmat[], 3, 1.0e-6) @check_error_code GridapPETSc.PETSC.KSPView(ksp[],C_NULL) end diff --git a/test/dev/bugfix.jl b/test/dev/bugfix.jl index 81c54da..2059627 100644 --- a/test/dev/bugfix.jl +++ b/test/dev/bugfix.jl @@ -11,7 +11,7 @@ using GridapGmsh using GridapMHD np = (2,2) -ranks = with_mpi() do distribute +ranks = with_debug() do distribute distribute(LinearIndices((prod(np),))) end @@ -20,12 +20,18 @@ mesh = Dict{Symbol,Any}( :base_mesh => "solid", ) +mesh = Dict{Symbol,Any}( + :mesher => :p4est_SG, + :base_mesh => "solid", + :num_refs => 0, +) + params = Dict{Symbol,Any}( :debug=>false, :solve=>true, :res_assemble=>false, :jac_assemble=>false, - :solver=> GridapMHD.default_solver_params(Val(:julia)) + :solver=> GridapMHD.default_solver_params(Val(:badia2024)) ) params[:fluid] = Dict{Symbol,Any}( :domain => "fluid", @@ -34,13 +40,14 @@ params[:fluid] = Dict{Symbol,Any}( :γ => 0.0, :f => VectorValue(0.0,0.0,0.0), :B => VectorValue(0.0,1.0,0.0), - :ζ => 0.0, + :ζ => 1.0, :convection => true, ) +u_in = GridapMHD.u_inlet(:parabolic,1.0,4.0,1.0) params[:bcs] = Dict{Symbol,Any}( :u => Dict( :tags => ["inlet", "wall"], - :values => [VectorValue(0.0, 0.0, 0.0), VectorValue(0.0, 0.0, 0.0)] + :values => [u_in, VectorValue(0.0, 0.0, 0.0)] ), :j => Dict( :tags => ["inlet", "outlet"], @@ -48,6 +55,10 @@ params[:bcs] = Dict{Symbol,Any}( ), ) params[:solid] = Dict(:domain => "wall",:σ => 1.0) +params[:fespaces] = Dict{Symbol,Any}( + :order_u => 2, + :order_j => 3, +) model = GridapMHD.expansion_mesh(mesh,ranks,params) @@ -63,19 +74,39 @@ mfs = GridapMHD._multi_field_style(params) V = MultiFieldFESpace([V_u,V_p,V_j,V_φ];style=mfs) U = MultiFieldFESpace([U_u,U_p,U_j,U_φ];style=mfs) -res, jac = GridapMHD.weak_form(params) +op = GridapMHD._fe_operator(U,V,params) +xh = GridapMHD._allocate_solution(op,params) +r = residual(op,xh) +j = jacobian(op,xh) -xh = zero(U); u = get_trial_fe_basis(U); v = get_fe_basis(V); contr = jac(xh,u,v) map(local_views(contr)) do a for strian in Gridap.CellData.get_domains(a) - scell_mat = Gridap.CellData.get_contribution(a,strian) - for mat in scell_mat - println(mat) - end + cell_mat = Gridap.CellData.get_contribution(a,strian) + #for mat in cell_mat + # println(mat) + #end + println(first(cell_mat)) end end + +res, jac = GridapMHD.weak_form(params) + +# res, jac = GridapMHD.weak_form(params) +# xh = zero(U); +# u = get_trial_fe_basis(U); +# v = get_fe_basis(V); +# +# contr = jac(xh,u,v) +# map(local_views(contr)) do a +# for strian in Gridap.CellData.get_domains(a) +# scell_mat = Gridap.CellData.get_contribution(a,strian) +# for mat in scell_mat +# println(mat) +# end +# end +# end #cell_mat = collect_cell_matrix(U, V, contr) From fb9b2866daa824b3ac9466c4985944e26cd63495 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Wed, 23 Oct 2024 08:48:45 +1100 Subject: [PATCH 2/3] Minor --- src/weakforms.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/weakforms.jl b/src/weakforms.jl index 5a41ba3..88e8536 100644 --- a/src/weakforms.jl +++ b/src/weakforms.jl @@ -308,7 +308,7 @@ function a_al(x,dy,ζ,Πp,dΩf,dΩ) v_u, v_p, v_j, v_φ = dy a_al_u_u(u,v_u,ζ,Πp,dΩf) + a_al_j_j(j,v_j,ζ,dΩ) end -a_al_u_u(u,v_u,ζ,Πp,dΩ) = ∫( ζ*(∇⋅u)*(∇⋅v_u) )*dΩ +a_al_u_u(u,v_u,ζ,Πp,dΩ) = ∫( ζ*Πp(u)*(∇⋅v_u) )*dΩ a_al_j_j(j,v_j,ζ,dΩ) = ∫( ζ*(∇⋅j)*(∇⋅v_j) )*dΩ # Solid equations From 00c03be24f149c8fd96cd67ce1f5720de3f86231 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Wed, 23 Oct 2024 14:57:39 +1100 Subject: [PATCH 3/3] Added RT scaling --- src/Applications/expansion.jl | 2 ++ src/GridapMHD.jl | 1 + src/Main.jl | 12 +++++++++++- src/parameters.jl | 14 ++++++++++++++ 4 files changed, 28 insertions(+), 1 deletion(-) diff --git a/src/Applications/expansion.jl b/src/Applications/expansion.jl index 4ecbf71..f45e2bc 100644 --- a/src/Applications/expansion.jl +++ b/src/Applications/expansion.jl @@ -63,6 +63,7 @@ function _expansion(; solid_coupling = :none, niter = nothing, convection = :true, + rt_scaling = false, savelines = false, petsc_options = "", ) @@ -119,6 +120,7 @@ function _expansion(; params[:fespaces] = Dict{Symbol,Any}( :order_u => order, :order_j => order_j, + :rt_scaling => rt_scaling ? 1.0/_get_mesh_size(model) : nothing ) params[:fluid] = Dict{Symbol,Any}( diff --git a/src/GridapMHD.jl b/src/GridapMHD.jl index 3c7c7d6..3151370 100644 --- a/src/GridapMHD.jl +++ b/src/GridapMHD.jl @@ -14,6 +14,7 @@ using DrWatson using Gridap using Gridap.Helpers, Gridap.Algebra, Gridap.CellData, Gridap.ReferenceFEs using Gridap.Geometry, Gridap.FESpaces, Gridap.MultiField, Gridap.ODEs +using Gridap.Fields using PartitionedArrays using PartitionedArrays: getany diff --git a/src/Main.jl b/src/Main.jl index 2ecc25c..6dd1ab3 100644 --- a/src/Main.jl +++ b/src/Main.jl @@ -269,7 +269,8 @@ function _fe_space(::Val{:j},params) uses_mg = space_uses_multigrid(params[:solver])[3] model = uses_mg ? params[:multigrid][:mh] : params[:model] - reffe_j = ReferenceFE(raviart_thomas,Float64,k-1) + phi = rt_scaling(model,params[:fespaces]) + reffe_j = ReferenceFE(raviart_thomas,Float64,k-1;basis_type=:jacobi,phi=phi) params[:fespaces][:reffe_j] = reffe_j j_bc = params[:bcs][:j][:values] @@ -449,3 +450,12 @@ 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)) + +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/parameters.jl b/src/parameters.jl index 000e0b9..c3a9be2 100644 --- a/src/parameters.jl +++ b/src/parameters.jl @@ -269,6 +269,7 @@ function params_fespaces(params::Dict{Symbol,Any}) :p_space => false, :p_constraint => false, :φ_constraint => false, + :rt_scaling => false, ) optional = Dict( :order_u => 2, @@ -276,6 +277,7 @@ function params_fespaces(params::Dict{Symbol,Any}) :p_space => :P, :p_constraint => nothing, :φ_constraint => nothing, + :rt_scaling => nothing, ) fespaces = _add_optional(params[:fespaces],mandatory,optional,params,"[:fespaces]") fespaces[:p_conformity] = p_conformity(params[:model],fespaces) @@ -306,6 +308,18 @@ end p_conformity(model::DiscreteModel,feparams) = p_conformity(Interior(model),feparams) p_conformity(model::GridapDistributed.DistributedDiscreteModel,feparams) = p_conformity(Interior(model),feparams) +# Scaling for Raviart-Thomas basis functions +function rt_scaling(model,feparams) + Dc = num_cell_dims(model) + if isnothing(feparams[:rt_scaling]) + phi = GenericField(identity) + else + ξ = feparams[:rt_scaling] + phi = AffineMap(ξ*one(TensorValue{Dc,Dc,Float64}),zero(VectorValue{Dc,Float64})) + end + return phi +end + """ Valid keys for `params[:multigrid]` are the following: