From 71d2885530f55f4e3f7c26934728ce8c4aa77290 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Thu, 26 Sep 2024 11:23:22 +1000 Subject: [PATCH 1/6] More fixes --- run.sh | 16 ---------------- src/Applications/expansion.jl | 10 +++++----- src/GridapMHD.jl | 3 ++- src/Main.jl | 6 +++++- src/Solvers/badia2024.jl | 1 + src/parameters.jl | 1 + src/weakforms.jl | 11 +++++++---- test/dev/buhler2006.jl | 31 +++++++++++++++++++++++++++++++ 8 files changed, 52 insertions(+), 27 deletions(-) delete mode 100755 run.sh create mode 100644 test/dev/buhler2006.jl diff --git a/run.sh b/run.sh deleted file mode 100755 index 8ec41a5..0000000 --- a/run.sh +++ /dev/null @@ -1,16 +0,0 @@ -mpiexec -n 4 julia --project=. -e' - using GridapMHD: hunt; - using SparseMatricesCSR; - using GridapPETSc; - - solver = Dict( - :solver => :badia2024, - :matrix_type => SparseMatrixCSR{0,PetscScalar,PetscInt}, - :vector_type => Vector{PetscScalar}, - :block_solvers => [:petsc_mumps,:cg_jacobi,:cg_jacobi], - :petsc_options => "-ksp_error_if_not_converged true -ksp_converged_reason", - :niter => 10, - ) - - hunt(backend=:mpi,np=(2,2),solver=solver,title="hunt",ζ=10.0) - ' diff --git a/src/Applications/expansion.jl b/src/Applications/expansion.jl index 9157687..1c3e3ad 100644 --- a/src/Applications/expansion.jl +++ b/src/Applications/expansion.jl @@ -115,11 +115,11 @@ function _expansion(; error("Unknown formulation") end - params[:fespaces] = Dict( + params[:fespaces] = Dict{Symbol,Any}( :order_u => order ) - params[:fluid] = Dict( + params[:fluid] = Dict{Symbol,Any}( :domain => nothing, # whole domain :α => α, :β => β, @@ -141,16 +141,16 @@ function _expansion(; ) if solid_coupling == :none - params[:bcs][:j] = Dict( + params[:bcs][:j] = Dict{Symbol,Any}( :tags => ["wall", "inlet", "outlet"], :values=>[VectorValue(0.0,0.0,0.0), VectorValue(0.0,0.0,0.0), VectorValue(0.0,0.0,0.0)], ) elseif solid_coupling == :thin_wall - params[:bcs][:j] = Dict( + params[:bcs][:j] = Dict{Symbol,Any}( :tags => ["inlet", "outlet"], :values=>[VectorValue(0.0,0.0,0.0), VectorValue(0.0,0.0,0.0)] ) - params[:bcs][:thin_wall] = Dict( + params[:bcs][:thin_wall] = Dict{Symbol,Any}( :τ=>τ, :cw=>cw, :domain => ["wall"], diff --git a/src/GridapMHD.jl b/src/GridapMHD.jl index dfa74ab..3c7c7d6 100644 --- a/src/GridapMHD.jl +++ b/src/GridapMHD.jl @@ -24,7 +24,8 @@ using GridapDistributed: i_am_main using GridapGmsh using GridapPETSc -using GridapSolvers, GridapSolvers.MultilevelTools, GridapSolvers.PatchBasedSmoothers +using GridapSolvers +using GridapSolvers.SolverInterfaces, GridapSolvers.MultilevelTools, GridapSolvers.PatchBasedSmoothers using GridapSolvers.LinearSolvers, GridapSolvers.NonlinearSolvers, GridapSolvers.BlockSolvers # Mesh generation diff --git a/src/Main.jl b/src/Main.jl index eeb60d2..5ccbb48 100644 --- a/src/Main.jl +++ b/src/Main.jl @@ -227,6 +227,7 @@ function _fe_space(::Val{:u},params) T = VectorValue{num_cell_dims(model),Float64} reffe_u = ReferenceFE(lagrangian,T,k) + params[:fespaces][:reffe_u] = reffe_u u_bc = params[:bcs][:u][:values] V_u = TestFESpace(Ωf,reffe_u;dirichlet_tags=params[:bcs][:u][:tags]) @@ -250,6 +251,7 @@ function _fe_space(::Val{:p},params) reffe_p = ReferenceFE(lagrangian,Float64,k-1;space=params[:fespaces][:p_space]) conformity = p_conformity(Ωf,params[:fespaces]) constraint = params[:fespaces][:p_constraint] + params[:fespaces][:reffe_p] = reffe_p V_p = TestFESpace(Ωf,reffe_p;conformity,constraint) U_p = _trial_fe_space(V_p,nothing,params[:transient]) @@ -263,6 +265,7 @@ function _fe_space(::Val{:j},params) model = uses_mg ? params[:multigrid][:mh] : params[:model] reffe_j = ReferenceFE(raviart_thomas,Float64,k-1) + params[:fespaces][:reffe_j] = reffe_j j_bc = params[:bcs][:j][:values] V_j = TestFESpace(model,reffe_j;dirichlet_tags=params[:bcs][:j][:tags]) @@ -285,6 +288,7 @@ function _fe_space(::Val{:φ},params) reffe_φ = ReferenceFE(lagrangian,Float64,k-1) conformity = :L2 constraint = params[:fespaces][:φ_constraint] + params[:fespaces][:reffe_φ] = reffe_φ V_φ = TestFESpace(model,reffe_φ;conformity,constraint) U_φ = _trial_fe_space(V_φ,nothing,params[:transient]) @@ -378,7 +382,7 @@ function _setup_trians!(params) else params[:multigrid][:Ωf] = Ωf params[:multigrid][:Ωs] = Ωs - params[:Ωf] = Ωs[1] + params[:Ωf] = Ωf[1] params[:Ωs] = Ωs[1] end end diff --git a/src/Solvers/badia2024.jl b/src/Solvers/badia2024.jl index e41c150..2665ca9 100644 --- a/src/Solvers/badia2024.jl +++ b/src/Solvers/badia2024.jl @@ -34,6 +34,7 @@ function Badia2024Solver(op::FEOperator,params) m = params[:solver][:niter] l_solver = FGMRESSolver(m,P;rtol=l_rtol,atol=1e-14,verbose=verbose,name="Global System - FGMRES + Badia2024") + SolverInterfaces.set_depth!(l_solver,2) # Nonlinear Solver nl_solver = GridapSolvers.NewtonSolver(l_solver,maxiter=1,atol=1e-14,rtol=nl_rtol,verbose=verbose) diff --git a/src/parameters.jl b/src/parameters.jl index 391f81e..3b003e3 100644 --- a/src/parameters.jl +++ b/src/parameters.jl @@ -211,6 +211,7 @@ function default_solver_params(::Val{:badia2024}) :block_solvers => [:petsc_mumps,:petsc_cg_jacobi,:petsc_cg_jacobi], :niter => 80, :rtol => 1e-5, + :initial_values => nothing, ) end diff --git a/src/weakforms.jl b/src/weakforms.jl index d5bc842..3d17069 100644 --- a/src/weakforms.jl +++ b/src/weakforms.jl @@ -22,6 +22,9 @@ function _weak_form(params,k) bcs_params = retrieve_bcs_params(params,k) params_φ, params_thin_wall, params_f, params_B, params_Λ = bcs_params + reffe_p = params[:fespaces][:reffe_p] + Πp = MultilevelTools.LocalProjectionMap(divergence,reffe_p,2*k) + function a(x,dy) r = a_mhd(x,dy,β,γ,B,σf,dΩf) for p in params_thin_wall @@ -157,7 +160,7 @@ retrieve_fluid_params(params,k) = retrieve_fluid_params(params[:model],params,k) function retrieve_fluid_params(model,params,k) fluid = params[:fluid] - Ωf = _interior(model,fluid[:domain]) + Ωf = params[:Ωf] dΩf = Measure(Ωf,2*k) α, β, γ, σf = fluid[:α], fluid[:β], fluid[:γ], fluid[:σ] @@ -173,7 +176,7 @@ retrieve_solid_params(params,k) = retrieve_solid_params(params[:model],params,k) function retrieve_solid_params(model,params,k) solid = params[:solid] if solid !== nothing - Ωs = _interior(model,solid[:domain]) + Ωs = params[:Ωs] dΩs = Measure(Ωs,2*k) σs = solid[:σ] return Ωs, dΩs, σs @@ -302,8 +305,8 @@ function a_al(x,dy,ζ,Πp,dΩ) v_u, v_p, v_j, v_φ = dy a_al_u_u(u,v_u,ζ,Πp,dΩ) + a_al_j_j(j,v_j,ζ,dΩ) end -a_al_u_u(u,v_u,ζ,Πp,dΩ) = ∫( ζ*Πp(∇⋅u)*Πp(∇⋅v_u) ) * dΩ -a_al_j_j(j,v_j,ζ,dΩ) = ∫( ζ*(∇⋅j)*(∇⋅v_j) ) * 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 diff --git a/test/dev/buhler2006.jl b/test/dev/buhler2006.jl new file mode 100644 index 0000000..39d7739 --- /dev/null +++ b/test/dev/buhler2006.jl @@ -0,0 +1,31 @@ +using GridapMHD: expansion; +using SparseMatricesCSR; +using GridapPETSc; + +solver = Dict( + :solver => :badia2024, + :matrix_type => SparseMatrixCSR{0,PetscScalar,PetscInt}, + :vector_type => Vector{PetscScalar}, + :block_solvers => [:petsc_mumps,:cg_jacobi,:cg_jacobi], + :petsc_options => "-ksp_error_if_not_converged true -ksp_converged_reason", +) + +expansion( + backend = :mpi, + np = 1, + mesh = "/scratch/bt62/jm3247/mhd-validation/experiments/buhler2006/meshes/expansion_Ha_1000_N_1000000_np_16.msh", + Ha = 1000, + N = 1000000, + ζ = 10.0, + cw = 0.028, + Z = 4.0, + b = 1.0, + inlet = :parabolic, + solid_coupling = :solid, + solver = solver, + savelines = true, + title = "expansion_Ha_1000_N_1000000_np_16", + path = "/scratch/bt62/jm3247/mhd-validation/experiments/buhler2006/data" +) + + From aea88ca8a3632f41294fa885c9913cbe4da9ba67 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Thu, 26 Sep 2024 11:33:54 +1000 Subject: [PATCH 2/6] Minor --- src/Solvers/badia2024.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Solvers/badia2024.jl b/src/Solvers/badia2024.jl index 2665ca9..e3fcde3 100644 --- a/src/Solvers/badia2024.jl +++ b/src/Solvers/badia2024.jl @@ -34,7 +34,8 @@ function Badia2024Solver(op::FEOperator,params) m = params[:solver][:niter] l_solver = FGMRESSolver(m,P;rtol=l_rtol,atol=1e-14,verbose=verbose,name="Global System - FGMRES + Badia2024") - SolverInterfaces.set_depth!(l_solver,2) + #SolverInterfaces.set_depth!(l_solver,2) + l_solver.log.depth = 2 # Nonlinear Solver nl_solver = GridapSolvers.NewtonSolver(l_solver,maxiter=1,atol=1e-14,rtol=nl_rtol,verbose=verbose) From 53ed8146a0f8e07e031171ce1574f103d494f034 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Mon, 30 Sep 2024 09:31:24 +1000 Subject: [PATCH 3/6] Minor --- src/Applications/expansion.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/Applications/expansion.jl b/src/Applications/expansion.jl index 1c3e3ad..34f34e1 100644 --- a/src/Applications/expansion.jl +++ b/src/Applications/expansion.jl @@ -120,7 +120,7 @@ function _expansion(; ) params[:fluid] = Dict{Symbol,Any}( - :domain => nothing, # whole domain + :domain => "fluid", :α => α, :β => β, :γ => γ, @@ -162,7 +162,6 @@ function _expansion(; :values=>[VectorValue(0.0,0.0,0.0), VectorValue(0.0,0.0,0.0)] ) params[:solid] = Dict(:domain => "wall",:σ => 1.0) - params[:fluid][:domain] = ["fluid"] end if μ > 0 From 154e94bf2677a09bf6b1f4554c1b9da6c3c0c645 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Mon, 30 Sep 2024 11:40:02 +1000 Subject: [PATCH 4/6] Fixed AL when solid --- src/Main.jl | 16 ++++++---------- src/Solvers/badia2024.jl | 10 ++++++---- src/weakforms.jl | 9 ++++++--- 3 files changed, 18 insertions(+), 17 deletions(-) diff --git a/src/Main.jl b/src/Main.jl index 5ccbb48..46154e2 100644 --- a/src/Main.jl +++ b/src/Main.jl @@ -348,15 +348,6 @@ end const DiscreteModelTypes = Union{Gridap.DiscreteModel,GridapDistributed.DistributedDiscreteModel} const TriangulationTypes = Union{Gridap.Triangulation,GridapDistributed.DistributedTriangulation} -function _fluid_mesh(model,domain::DiscreteModelTypes) - msg = "params[:fluid][:domain] is a discrete model, but params[:fluid][:domain]===params[:model] is not true." - @assert model === domain msg - return domain -end -_fluid_mesh(model,domain::TriangulationTypes) = domain -_fluid_mesh(model,domain::Nothing) = model # This should be removed, but Gridap needs fixes -_fluid_mesh(model,domain) = Interior(model,tags=domain) - _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 @@ -370,18 +361,23 @@ _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 = _fluid_mesh(params[:model],fluid[:domain]) + Ω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 diff --git a/src/Solvers/badia2024.jl b/src/Solvers/badia2024.jl index e3fcde3..14b7c58 100644 --- a/src/Solvers/badia2024.jl +++ b/src/Solvers/badia2024.jl @@ -4,11 +4,13 @@ function Badia2024Solver(op::FEOperator,params) # Preconditioner model = params[:model] k = params[:fespaces][:k] - Ωf = _interior(model,params[:fluid][:domain]) - dΩ = Measure(Ωf,2*k) + Ω = params[:Ω] + dΩ = Measure(Ω,2*k) + Ωf = params[:Ωf] + dΩf = Measure(Ωf,2*k) α_p = -1.0/(params[:fluid][:β] + params[:fluid][:ζ]) α_φ = -1.0/(1.0 + params[:fluid][:ζ]) - a_Ip(p,v_p) = ∫(α_p*p*v_p)*dΩ + a_Ip(p,v_p) = ∫(α_p*p*v_p)*dΩf a_Iφ(φ,v_φ) = ∫(α_φ*φ*v_φ)*dΩ U_u, U_p, U_j, U_φ = get_trial(op) @@ -38,6 +40,6 @@ function Badia2024Solver(op::FEOperator,params) l_solver.log.depth = 2 # Nonlinear Solver - nl_solver = GridapSolvers.NewtonSolver(l_solver,maxiter=1,atol=1e-14,rtol=nl_rtol,verbose=verbose) + nl_solver = GridapSolvers.NewtonSolver(l_solver,maxiter=10,atol=1e-14,rtol=nl_rtol,verbose=verbose) return nl_solver end diff --git a/src/weakforms.jl b/src/weakforms.jl index 3d17069..88e8536 100644 --- a/src/weakforms.jl +++ b/src/weakforms.jl @@ -13,6 +13,9 @@ function weak_form(params,k) end function _weak_form(params,k) + Ω = params[:Ω] + dΩ = Measure(Ω,2*k) + fluid = params[:fluid] Ωf, dΩf, α, β, γ, σf, f, B, ζ, g = retrieve_fluid_params(params,k) @@ -40,7 +43,7 @@ function _weak_form(params,k) r = r + a_solid(x,dy,σs,dΩs) end if abs(ζ) > eps(typeof(ζ)) - r = r + a_al(x,dy,ζ,Πp,dΩf) + r = r + a_al(x,dy,ζ,Πp,dΩf,dΩ) end r end @@ -300,10 +303,10 @@ dc_mhd_u_u(u,du,v_u,α,dΩ) = ∫( α*v_u⋅( (conv∘(u,∇(du))) + (conv∘(du # Augmented lagrangian -function a_al(x,dy,ζ,Πp,dΩ) +function a_al(x,dy,ζ,Πp,dΩf,dΩ) u, p, j, φ = x v_u, v_p, v_j, v_φ = dy - a_al_u_u(u,v_u,ζ,Πp,dΩ) + a_al_j_j(j,v_j,ζ,dΩ) + 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Ω) = ∫( ζ*Πp(u)*(∇⋅v_u) )*dΩ a_al_j_j(j,v_j,ζ,dΩ) = ∫( ζ*(∇⋅j)*(∇⋅v_j) )*dΩ From 2ee431679083420581161bfa5330f2359394436f Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Mon, 30 Sep 2024 15:36:24 +1000 Subject: [PATCH 5/6] Minor --- src/parameters.jl | 2 ++ src/weakforms.jl | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/parameters.jl b/src/parameters.jl index 3b003e3..20cd51b 100644 --- a/src/parameters.jl +++ b/src/parameters.jl @@ -162,6 +162,8 @@ function params_solver(params::Dict{Symbol,Any}) return solver end +default_solver_params(s::Symbol) = default_solver_params(Val(s)) + function default_solver_params(::Val{:julia}) return Dict( :solver => :julia, diff --git a/src/weakforms.jl b/src/weakforms.jl index 88e8536..5a41ba3 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Ω) = ∫( ζ*Πp(u)*(∇⋅v_u) )*dΩ +a_al_u_u(u,v_u,ζ,Πp,dΩ) = ∫( ζ*(∇⋅u)*(∇⋅v_u) )*dΩ a_al_j_j(j,v_j,ζ,dΩ) = ∫( ζ*(∇⋅j)*(∇⋅v_j) )*dΩ # Solid equations From 0c46fd31eacb9efb4f1bc1fa75e00414cc9851a3 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Tue, 1 Oct 2024 17:30:05 +1000 Subject: [PATCH 6/6] More fixes --- Project.toml | 1 - src/Solvers/badia2024.jl | 2 +- src/Solvers/petsc.jl | 6 ++-- src/parameters.jl | 2 +- test/seq/expansion_badia2024_tests.jl | 49 +++++++++++++-------------- 5 files changed, 29 insertions(+), 31 deletions(-) diff --git a/Project.toml b/Project.toml index 18089f2..b000522 100644 --- a/Project.toml +++ b/Project.toml @@ -30,7 +30,6 @@ gmsh_jll = "630162c2-fc9b-58b3-9910-8442a8a132e6" BSON = "0.3" FileIO = "1" Gridap = "0.18" -GridapDistributed = "0.4.0" GridapGmsh = "0.7.2" GridapP4est = "0.3.7" GridapPETSc = "0.5.0" diff --git a/src/Solvers/badia2024.jl b/src/Solvers/badia2024.jl index 14b7c58..6ddc0ec 100644 --- a/src/Solvers/badia2024.jl +++ b/src/Solvers/badia2024.jl @@ -16,7 +16,7 @@ function Badia2024Solver(op::FEOperator,params) U_u, U_p, U_j, U_φ = get_trial(op) V_u, V_p, V_j, V_φ = get_test(op) - diag_solvers = map(s -> get_block_solver(Val(s),params),params[:solver][:block_solvers]) + diag_solvers = map(s -> get_block_solver(Val(s),params), params[:solver][:block_solvers]) uj_block = NonlinearSystemBlock(1) p_block = BiformBlock(a_Ip,U_p,V_p) diff --git a/src/Solvers/petsc.jl b/src/Solvers/petsc.jl index 4bd12af..c6b4fb1 100644 --- a/src/Solvers/petsc.jl +++ b/src/Solvers/petsc.jl @@ -17,10 +17,10 @@ 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[], 28, 2) - # @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[], 29, 2) + @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) + @check_error_code GridapPETSc.PETSC.KSPView(ksp[],C_NULL) end # Standalone AMG solver, 5 iterations diff --git a/src/parameters.jl b/src/parameters.jl index 20cd51b..000e0b9 100644 --- a/src/parameters.jl +++ b/src/parameters.jl @@ -589,5 +589,5 @@ function params_bcs_stabilization(params::Dict{Symbol,Any}) :μ=>true, ) optional = Dict(:domain=>params[:fluid][:domain]) - _check_mandatory_and_add_optional_weak(params[:bcs][:stabilization],mandatory,optional,params,"[:bcs][:thin_wall]") + _check_mandatory_and_add_optional_weak(params[:bcs][:stabilization],mandatory,optional,params,"[:bcs][:stabilization]") end diff --git a/test/seq/expansion_badia2024_tests.jl b/test/seq/expansion_badia2024_tests.jl index d5f3bf3..98e4b60 100644 --- a/test/seq/expansion_badia2024_tests.jl +++ b/test/seq/expansion_badia2024_tests.jl @@ -2,38 +2,37 @@ module ExpansionBadia2024TestsSequential using GridapMHD: expansion using SparseArrays +using SparseMatricesCSR +using GridapPETSc -cw = 0.028 Re = 1.0 -Ha = 100.0 +Ha = 10.0 N = Ha^2/Re -mesh = Dict( - :mesher => :p4est_MG, - :base_mesh => "coarse", - :num_refs_coarse => 0, - :ranks_per_level => [1,1], -) +# solver = Dict( +# :solver => :badia2024, +# :matrix_type => SparseMatrixCSC{Float64,Int64}, +# :vector_type => Vector{Float64}, +# :block_solvers => [:julia,:cg_jacobi,:cg_jacobi], +# :petsc_options => "-ksp_error_if_not_converged true -ksp_converged_reason" +# ) + +# solver = Dict( +# :solver => :badia2024, +# :matrix_type => SparseMatrixCSR{0,PetscScalar,PetscInt}, +# :vector_type => Vector{PetscScalar}, +# :block_solvers => [:petsc_from_options,:cg_jacobi,:cg_jacobi], +# :petsc_options => "-ksp_error_if_not_converged true -ksp_converged_reason -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps", +# ) + solver = Dict( :solver => :badia2024, - :matrix_type => SparseMatrixCSC{Float64,Int64}, - :vector_type => Vector{Float64}, - :block_solvers => [:gmg,:cg_jacobi,:cg_jacobi], - :petsc_options => "-ksp_error_if_not_converged true -ksp_converged_reason" + :matrix_type => SparseMatrixCSR{0,PetscScalar,PetscInt}, + :vector_type => Vector{PetscScalar}, + :block_solvers => [:petsc_mumps,:cg_jacobi,:cg_jacobi], + :petsc_options => "-ksp_error_if_not_converged true -ksp_converged_reason", ) -expansion(np=1,backend=:mpi,mesh=mesh,solver=solver,order=2,ζ=100.0,N=N,Ha=Ha,cw=cw,title="ExpansionGMG",formulation=:mhd) -""" -# ReferenceSolution -expansion(np=1,backend=:mpi,solver=:julia,order=2,ζ=0.0,N=N,Ha=Ha,cw=cw,title="ExpansionRef") +expansion(np=4,backend=:mpi,mesh="710",solver=solver,order=2,ζ=10.0,N=N,Ha=Ha,title="Expansion") -expansion(np=1,backend=:mpi,solver=:julia,order=2,ζ=0.0,N=N,Ha=Ha,cw=cw,title="ExpansionRefBis",mesh=mesh) - -meshbis = Dict( - :mesher => :p4est_SG, - :base_mesh => "coarse", - :num_refs => 1 -) -expansion(np=1,backend=:mpi,solver=:julia,order=2,ζ=0.0,N=N,Ha=Ha,cw=cw,title="ExpansionRefBis",mesh=meshbis) -""" end # module \ No newline at end of file