Merge branch 'develop' of github.com:gridapapps/GridapMHD.jl into dev…

…elop

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"

src/Applications/expansion.jl

@@ -115,12 +115,12 @@ function _expansion(;
     error("Unknown formulation")
   end
 
-  params[:fespaces] = Dict(
+  params[:fespaces] = Dict{Symbol,Any}(
     :order_u => order
   )
 
-  params[:fluid] = Dict(
-    :domain => nothing, # whole domain
+  params[:fluid] = Dict{Symbol,Any}(
+    :domain => "fluid",
     :α => α,
     :β => β,
     :γ => γ,
@@ -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"],
@@ -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

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

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])
@@ -344,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
@@ -366,19 +361,24 @@ _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[:Ωf] = Ωs[1]
+    params[:Ω]  = Ω[1]
+    params[:Ωf] = Ωf[1]
     params[:Ωs] = Ωs[1]
   end
 end

src/Solvers/badia2024.jl

@@ -4,17 +4,19 @@ 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)
   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)
@@ -34,8 +36,10 @@ 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)
+  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

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

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,
@@ -211,6 +213,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
 
@@ -586,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

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)
 
@@ -22,6 +25,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
@@ -37,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
@@ -157,7 +163,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 +179,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
@@ -297,13 +303,13 @@ 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)*Π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Ω) = ∫( ζ*(∇⋅u)*(∇⋅v_u) )*dΩ
+a_al_j_j(j,v_j,ζ,dΩ) = ∫( ζ*(∇⋅j)*(∇⋅v_j) )*dΩ
 
 # Solid equations
 

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"
+)
+
+

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