Added Li2019 preconditioner to the library

Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.9.2"
 manifest_format = "2.0"
-project_hash = "fc82c8fcba89b0e1a03d17a5669f11d0aaedef3f"
+project_hash = "a707abf4993889e607550240d993374a45ea273e"
 
 [[deps.AbstractFFTs]]
 deps = ["LinearAlgebra"]
@@ -817,9 +817,9 @@ version = "0.4.19"
 
 [[deps.PrecompileTools]]
 deps = ["Preferences"]
-git-tree-sha1 = "9673d39decc5feece56ef3940e5dafba15ba0f81"
+git-tree-sha1 = "03b4c25b43cb84cee5c90aa9b5ea0a78fd848d2f"
 uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
-version = "1.1.2"
+version = "1.2.0"
 
 [[deps.Preferences]]
 deps = ["TOML"]

Project.toml

@@ -1,21 +1,8 @@
-authors = ["Jesús Bonilla", "Francesc Verdugo", "Large Scale Scientific Computing"]
 name = "GridapMHD"
 uuid = "afadc09e-922c-4d62-8330-7a6898d360cc"
+authors = ["Jesús Bonilla", "Francesc Verdugo", "Large Scale Scientific Computing"]
 version = "0.4.0"
 
-[compat]
-BSON = "0.3"
-FileIO = "1"
-Gridap = "0.17.18"
-GridapDistributed = "0.3"
-GridapGmsh = "0.7"
-GridapPETSc = "0.5"
-MPI = "0.20"
-PackageCompiler = "2"
-PartitionedArrays = "0.3"
-SparseMatricesCSR = "0.6.6"
-julia = "1.7"
-
 [deps]
 BSON = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
@@ -26,6 +13,7 @@ GridapDistributed = "f9701e48-63b3-45aa-9a63-9bc6c271f355"
 GridapGmsh = "3025c34a-b394-11e9-2a55-3fee550c04c8"
 GridapPETSc = "bcdc36c2-0c3e-11ea-095a-c9dadae499f1"
 GridapSolvers = "6d3209ee-5e3c-4db7-a716-942eb12ed534"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
 PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d"
 PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9"
@@ -34,6 +22,19 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1"
 
+[compat]
+BSON = "0.3"
+FileIO = "1"
+Gridap = "0.17.18"
+GridapDistributed = "0.3"
+GridapGmsh = "0.7"
+GridapPETSc = "0.5"
+MPI = "0.20"
+PackageCompiler = "2"
+PartitionedArrays = "0.3"
+SparseMatricesCSR = "0.6.6"
+julia = "1.7"
+
 [extras]
 MPIPreferences = "3da0fdf6-3ccc-4f1b-acd9-58baa6c99267"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

src/GridapMHD.jl

@@ -5,25 +5,33 @@ module GridapMHD
 # using MPI
 
 using Random
+using LinearAlgebra
 using SparseArrays
 using SparseMatricesCSR
-using PartitionedArrays
+using BlockArrays
+
+using FileIO
+using BSON
+
 using Gridap
 using Gridap.Helpers
 using Gridap.Algebra
 using Gridap.CellData
 using Gridap.ReferenceFEs
 using Gridap.Geometry
 using Gridap.FESpaces
+using Gridap.MultiField
+
 using GridapDistributed
-using GridapGmsh
 using GridapPETSc
 using GridapPETSc.PETSC
-using FileIO
-using BSON
 using GridapGmsh
+
 using GridapSolvers
+using GridapSolvers.LinearSolvers
+using GridapSolvers.LinearSolvers: allocate_col_vector, allocate_row_vector
 
+using PartitionedArrays
 using PartitionedArrays: getany
 
 include("Main.jl")
@@ -36,6 +44,6 @@ include("expansion.jl")
 
 include("cavity.jl")
 
-include("block_solvers.jl")
+include("block_solver_li2019.jl")
 
 end # module

src/Main.jl

@@ -110,7 +110,7 @@ function add_default_params(_params)
   params
 end
 
-default_ptimer(model) = PTimer(DebugArray(collect(LinearIndices(1))))
+default_ptimer(model) = PTimer(DebugArray(LinearIndices((1,))))
 default_ptimer(model::GridapDistributed.DistributedDiscreteModel) = PTimer(get_parts(model))
 
 """

src/block_solver_li2019.jl

@@ -0,0 +1,106 @@
+
+"""
+  This preconditioner is based on [(Li,2019)](https://doi.org/10.1137/19M1260372)
+"""
+struct LI2019_Solver <: Gridap.Algebra.LinearSolver
+  Dj_solver
+  Fk_solver
+  Δp_solver
+  Ip_solver
+  Iφ_solver
+  Dj
+  Δp
+  Ip
+  Ij
+  Iφ
+  params
+end
+
+struct LI2019_SS <: Gridap.Algebra.SymbolicSetup
+  solver::LI2019_Solver
+end
+
+function Gridap.Algebra.symbolic_setup(solver::LI2019_Solver, A::AbstractMatrix)
+  return LI2019_SS(solver)
+end
+
+mutable struct LI2019_NS <: Gridap.Algebra.NumericalSetup
+  solver
+  Dj_ns
+  Fk_ns
+  Δp_ns
+  Ip_ns
+  Iφ_ns
+  sysmat
+  caches
+end
+
+function allocate_caches(solver::LI2019_Solver,A::BlockMatrix)
+  du = allocate_col_vector(A[Block(1,1)])
+  dp = allocate_col_vector(A[Block(2,2)])
+  dj = allocate_col_vector(A[Block(3,3)])
+  dφ = allocate_col_vector(A[Block(4,4)])
+  return du, dp, dj, dφ
+end
+
+function Gridap.Algebra.numerical_setup(ss::LI2019_SS, A::BlockMatrix)
+  solver = ss.solver
+
+  Fu = A[Block(1,1)]; K = A[Block(3,1)]; Kᵗ = A[Block(1,3)]; κ = solver.params[:fluid][:γ]
+  Fk = Fu - (1.0/κ^2) * Kᵗ * K
+
+  Dj_ns = numerical_setup(symbolic_setup(solver.Dj_solver,solver.Dj),solver.Dj)
+  Fk_ns = numerical_setup(symbolic_setup(solver.Fk_solver,Fk),Fk)
+  Δp_ns = numerical_setup(symbolic_setup(solver.Δp_solver,solver.Δp),solver.Δp)
+  Ip_ns = numerical_setup(symbolic_setup(solver.Δp_solver,solver.Ip),solver.Ip)
+  Iφ_ns = numerical_setup(symbolic_setup(solver.Iφ_solver,solver.Iφ),solver.Iφ)
+  caches = allocate_caches(solver,A)
+  return LI2019_NS(ss.solver,Dj_ns,Fk_ns,Δp_ns,Ip_ns,Iφ_ns,A,caches)
+end
+
+function Gridap.Algebra.numerical_setup!(ns::LI2019_NS, A::BlockMatrix)
+  solver = ns.solver
+
+  #! Pattern of matrix changes, so we need to recompute everything.
+  # This will get fixed when we are using iterative solvers for Fk
+  Fu = A[Block(1,1)]; K = A[Block(3,1)]; Kᵗ = A[Block(1,3)]; κ = solver.params[:fluid][:γ]
+  Fk = Fu - (1.0/κ^2) * Kᵗ * K
+  # numerical_setup!(ns.Fk_ns,Fk)
+
+  ns.Fk_ns  = numerical_setup(symbolic_setup(solver.Fk_solver,Fk),Fk)
+  ns.sysmat = A
+
+  return ns
+end
+
+# Follows Algorithm 4.1 in (Li,2019)
+function Gridap.Algebra.solve!(x::BlockVector,ns::LI2019_NS,b::BlockVector)
+  sysmat, caches, params = ns.sysmat, ns.caches, ns.solver.params
+  fluid = params[:fluid]; ζ = params[:ζ]; iRe = fluid[:β]
+  κ = fluid[:γ]; α1 = ζ + iRe;
+
+  bu, bp, bj, bφ = blocks(b)
+  u , p , j , φ  = blocks(x)
+  du, dp, dj, dφ = caches
+
+  # Solve for p
+  solve!(p,ns.Δp_ns,bp)
+  solve!(dp,ns.Ip_ns,bp)
+  p .= -α1 .* dp .- p
+
+  #  Solve for φ
+  #dφ .= -bφ
+  solve!(φ,ns.Iφ_ns,bφ)
+
+  # Solve for u
+  copy!(du,bu); mul!(du,sysmat[Block(1,2)],p,-1.0,1.0) # du = bu - Aup * p
+  solve!(u,ns.Fk_ns,du)                                # u = Fu \ (bu - Aup * p)
+
+  # Solve for j
+  copy!(dj,bj)
+  mul!(dj,sysmat[Block(3,1)],u,-2.0,1.0) # dj = bj - 2.0 * Aju * u
+  mul!(dj,sysmat[Block(3,4)],φ,-2.0,1.0) # dj = bj - 2.0 * Aju * u - 2.0 * Ajφ * φ
+  solve!(j,ns.Dj_ns,dj)                  # j = Dj \ (bj - 2.0 * Aju * u - 2.0 * Ajφ * φ)
+
+  return x
+end

src/cavity.jl

@@ -8,18 +8,18 @@ function cavity(;
 
   if isa(backend,Nothing)
     @assert isa(np,Nothing)
-    info, t = _cavity(;title=title,path=path,mesh=mesh,kwargs...)
+    info, t = _cavity(;title=title,path=path,kwargs...)
   else
     @assert backend ∈ [:sequential,:mpi]
     @assert !isa(np,Nothing)
     np = isa(np,Int) ? (np,) : np
     if backend == :sequential
       info,t = with_debug() do distribute
-        _cavity(;distribute=distribute,rank_partition=np,title=_title,path=path,mesh=mesh,kwargs...)
+        _cavity(;distribute=distribute,rank_partition=np,title=title,path=path,kwargs...)
       end
     else
       info,t = with_mpi() do distribute
-        _cavity(;distribute=distribute,rank_partition=np,title=_title,path=path,mesh=mesh,kwargs...)
+        _cavity(;distribute=distribute,rank_partition=np,title=title,path=path,kwargs...)
       end
     end
   end
@@ -50,7 +50,9 @@ function _cavity(;
   B0=norm(B),
   vtk=true,
   title="Cavity",
+  path='.',
   solver=:julia,
+  verbose=true,
   )
 
   @assert length(nc) == 3
@@ -61,7 +63,7 @@ function _cavity(;
     rank_partition = (1,)
     distribute = DebugArray
   end
-  @assert length(rank_partition) == length(nc)
+  @assert is_serial || (length(rank_partition) == length(nc))
   parts = distribute(LinearIndices((prod(rank_partition),)))
 
   t = PTimer(parts,verbose=verbose)
@@ -103,6 +105,7 @@ function _cavity(;
       :solve => true,
       :res_assemble => false,
       :jac_assemble => false,
+      :check_valid => false,
       :model => model,
       :fluid => Dict(
           :domain => model,
@@ -115,7 +118,8 @@ function _cavity(;
       :bcs => Dict(
           :u => Dict(:tags => ["wall", "lid"], :values => [uw, ul]),
           :j => Dict(:tags => "insulating", :values => ji),
-      )
+      ),
+      :k => 2,
       :ζ => 0.0 # Augmented-Lagragian term 
   )
 
@@ -124,7 +128,7 @@ function _cavity(;
 
   tic!(t; barrier=true)
   # ReferenceFEs
-  k = 2
+  k = params[:k]
   T = Float64
   model = params[:model]
   D = num_cell_dims(model)
@@ -158,7 +162,7 @@ function _cavity(;
   op = FEOperator(res, jac, U, V, assem)
 
   if solver === :block_gmres
-    xh = block_gmres_solver(op,U,V)
+    xh = block_gmres_solver(op,U,V,Ω,params)
   else
     xh = zero(U)
     solver = NLSolver(show_trace=true, method=:newton)
@@ -173,18 +177,69 @@ function _cavity(;
     ph = (ρ * u0^2) * p̄h
     jh = (σ * u0 * B0) * j̄h
     φh = (u0 * B0 * L) * φ̄h
-    writevtk(Ω, title, order=2, cellfields=["uh" => uh, "ph" => ph, "jh" => jh, "phi" => φh])
+    writevtk(Ω, joinpath(path,title), order=2, cellfields=["uh" => uh, "ph" => ph, "jh" => jh, "phi" => φh])
     toc!(t, "vtk")
   end
 
   info = Dict{Symbol,Any}()
-  info[:ncells] = num_cells(model)
+  info[:ncells]  = num_cells(model)
   info[:ndofs_u] = length(get_free_dof_values(ūh))
   info[:ndofs_p] = length(get_free_dof_values(p̄h))
   info[:ndofs_j] = length(get_free_dof_values(j̄h))
   info[:ndofs_φ] = length(get_free_dof_values(φ̄h))
-  info[:ndofs] = length(get_free_dof_values(xh))
-  info[:Re] = Re
-  info[:Ha] = Ha
-  save("$title.bson", info)
-end
+  info[:ndofs]   = sum([info[:ndofs_u], info[:ndofs_p], info[:ndofs_j], info[:ndofs_φ]])
+  info[:Re]      = Re
+  info[:Ha]      = Ha
+
+  return info, t
+end
+
+
+function block_gmres_solver(op,U,V,Ω,params)
+  U_u, U_p, U_j, U_φ = U
+  V_u, V_p, V_j, V_φ = V
+
+  al_op = Gridap.FESpaces.get_algebraic_operator(op)
+
+  k  = params[:k]
+  γ  = params[:fluid][:γ]
+  dΩ = Measure(Ω,2*k)
+  Dj = assemble_matrix((j,v_j) -> ∫(γ*j⋅v_j + γ*(∇⋅j)⋅(∇⋅v_j))*dΩ ,U_j,V_j)
+  Ij = assemble_matrix((j,v_j) -> ∫(j⋅v_j)*dΩ ,U_j,V_j)
+  Δp = assemble_matrix((p,v_p) -> ∫(∇(p)⋅∇(v_p))*dΩ ,U_p,V_p)
+  Ip = assemble_matrix((p,v_p) -> ∫(p*v_p)*dΩ,V_p,V_p)
+  Iφ = assemble_matrix((φ,v_φ) -> ∫(-γ*φ*v_φ)*dΩ ,U_φ,V_φ)
+
+  Dj_solver = LUSolver()
+  Fk_solver = LUSolver()
+  Δp_solver = LUSolver()
+  Ip_solver = LUSolver()
+  Iφ_solver = LUSolver()
+
+  block_solvers = [Dj_solver,Fk_solver,Δp_solver,Ip_solver,Iφ_solver]
+  block_mats = [Dj,Δp,Ip,Ij,Iφ]
+  P = LI2019_Solver(block_solvers...,block_mats...,params)
+
+  sysmat_solver = GMRESSolver(300,P,1e-8)
+
+  # Gridap's Newton-Raphson solver
+  xh = zero(U)
+  sysmat = jacobian(op,xh)
+  sysmat_ns = numerical_setup(symbolic_setup(sysmat_solver,sysmat),sysmat)
+
+  x  = allocate_col_vector(sysmat)
+  dx = allocate_col_vector(sysmat)
+  b  = allocate_col_vector(sysmat)
+
+  A = allocate_jacobian(al_op,x)
+  nlsolver = NewtonRaphsonSolver(sysmat_solver,1e-5,10)
+  nlsolver_cache = Gridap.Algebra.NewtonRaphsonCache(A,b,dx,sysmat_ns)
+  solve!(x,nlsolver,al_op,nlsolver_cache)
+
+  ūh = FEFunction(U_u,x.blocks[1])
+  p̄h = FEFunction(U_p,x.blocks[2])
+  j̄h = FEFunction(U_j,x.blocks[3])
+  φ̄h = FEFunction(U_φ,x.blocks[4])
+  return [ūh, p̄h, j̄h, φ̄h]
+end
+