Bump to Gridap v0.17

Project.toml

@@ -14,7 +14,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 [compat]
 FastClosures = "0.2, 0.3"
 ForwardDiff = "0.9.0, 0.10.0"
-Gridap = "0.15.5"
+Gridap = "^0.17.8"
 NLPModels = "0.15, 0.16, 0.17, 0.18"
 julia = "^1.6.0"
 

src/GridapPDENLPModel.jl

@@ -16,10 +16,20 @@ The following arguments are accepted:
 - nvar_con : number of dofs in the control functions
 - nparam: : number of real unknowns
 - nnzh_obj : number of nonzeros elements in the objective hessian
-- Hrows : store the structure for the hessian of the lagrangian
-- Hcols : store the structure for the hessian of the lagrangian
-- Jrows : store the structure for the hessian of the jacobian
-- Jcols : store the structure for the hessian of the jacobian
+- Hobjkrows : store the structure for the hessian of the objective w.r.t. k
+- Hobjkcols : store the structure for the hessian of the objective w.r.t. k
+- Hobjrows : store the structure for the hessian of the objective
+- Hobjcols : store the structure for the hessian of the objective
+- Hkrows : store the structure for the hessian of the constraints w.r.t. k
+- Hkcols : store the structure for the hessian of the constraints w.r.t. k
+- Hrows : store the structure for the hessian of the constraints
+- Hcols : store the structure for the hessian of the constraints
+- Jkrows : store the structure for the hessian of the jacobian w.r.t. k
+- Jkcols : store the structure for the hessian of the jacobian w.r.t. k
+- Jyrows : store the structure for the hessian of the jacobian w.r.t. y
+- Jycols : store the structure for the hessian of the jacobian w.r.t. y
+- Jurows : store the structure for the hessian of the jacobian w.r.t. u
+- Jucols : store the structure for the hessian of the jacobian w.r.t. u
 """
 struct PDENLPMeta{NRJ <: AbstractEnergyTerm, Op <: Union{FEOperator, Nothing}}
   # For the objective function
@@ -42,10 +52,20 @@ struct PDENLPMeta{NRJ <: AbstractEnergyTerm, Op <: Union{FEOperator, Nothing}}
   nnzh_obj::Int
 
   # store the structure for hessian and jacobian matrix
+  Hobjkrows::AbstractVector{Int}
+  Hobjkcols::AbstractVector{Int}
+  Hobjrows::AbstractVector{Int}
+  Hobjcols::AbstractVector{Int}
+  Hkrows::AbstractVector{Int}
+  Hkcols::AbstractVector{Int}
   Hrows::AbstractVector{Int}
   Hcols::AbstractVector{Int}
-  Jrows::AbstractVector{Int}
-  Jcols::AbstractVector{Int}
+  Jkrows::AbstractVector{Int}
+  Jkcols::AbstractVector{Int}
+  Jyrows::AbstractVector{Int}
+  Jycols::AbstractVector{Int}
+  Jurows::AbstractVector{Int}
+  Jucols::AbstractVector{Int}
 end
 
 struct PDEWorkspace{T, S, M}
@@ -187,9 +207,35 @@ function hess_coord!(
     luh = _obj_integral(nlp.pdemeta.tnrj, κ, yu)
     lag_hess = Gridap.FESpaces._hessian(x -> _obj_integral(nlp.pdemeta.tnrj, κ, x), yu, luh)
 
-    matdata = Gridap.FESpaces.collect_cell_matrix(lag_hess)
+    matdata = Gridap.FESpaces.collect_cell_matrix(nlp.pdemeta.Y, nlp.pdemeta.X, lag_hess)
     assem = SparseMatrixAssembler(nlp.pdemeta.Y, nlp.pdemeta.X)
-    nini = fill_hess_coo_numeric!(vals, assem, matdata, n = nini)
+    ###############################################################
+    # TO BE IMPROVED
+
+    m1 = Gridap.FESpaces.nz_counter(
+      Gridap.FESpaces.get_matrix_builder(assem),
+      (Gridap.FESpaces.get_rows(assem), Gridap.FESpaces.get_cols(assem)),
+    ) # Gridap.Algebra.CounterCS
+    Gridap.FESpaces.symbolic_loop_matrix!(m1, assem, matdata)
+    m2 = Gridap.FESpaces.nz_allocation(m1) # Gridap.Algebra.InserterCSC
+    #@show size(Gridap.FESpaces.create_from_nz(m2))
+    Gridap.FESpaces.numeric_loop_matrix!(m2, assem, matdata)
+    m3 = sparse(LowerTriangular(Gridap.FESpaces.create_from_nz(m2)))
+    _, _, v = findnz(m3)
+    vals[(nini + 1):(nini + length(v))] .= v
+    nini += length(v)
+
+    ##############################################################
+    #=
+    nnzh = length(nlp.pdemeta.Hobjrows)
+    @show nlp.meta.nvar, nlp.pdemeta.nparam
+    ms = sparse(nlp.pdemeta.Hobjrows, nlp.pdemeta.Hobjcols, zeros(T, nnzh), nlp.meta.nvar - nlp.pdemeta.nparam, nlp.meta.nvar - nlp.pdemeta.nparam) # SparseArrays.sparse!
+    Gridap.FESpaces.numeric_loop_matrix!(ms, assem, matdata)
+    ms_lower = sparse(LowerTriangular(Gridap.FESpaces.create_from_nz(ms)))
+    vals[(nini + 1):(nini + nnzh)] .= ms_lower.nzval
+    nini += nnzh
+    =#
+
   end
 
   if nini < nlp.meta.nnzh
@@ -224,7 +270,8 @@ function hprod!(
 end
 
 function hess_structure(nlp::GridapPDENLPModel)
-  return (nlp.pdemeta.Hrows, nlp.pdemeta.Hcols)
+  np = nlp.pdemeta.nparam
+  return (vcat(nlp.pdemeta.Hobjkrows, nlp.pdemeta.Hobjrows .+ np, nlp.pdemeta.Hkrows, nlp.pdemeta.Hrows .+ np), vcat(nlp.pdemeta.Hobjkcols, nlp.pdemeta.Hobjcols .+ np, nlp.pdemeta.Hkcols, nlp.pdemeta.Hcols .+ np))
 end
 
 function hess_structure!(
@@ -233,8 +280,9 @@ function hess_structure!(
   cols::AbstractVector{T},
 ) where {T <: Integer}
   @lencheck nlp.meta.nnzh rows cols
-  rows .= T.(nlp.pdemeta.Hrows)
-  cols .= T.(nlp.pdemeta.Hcols)
+  np = nlp.pdemeta.nparam
+  rows .= T.(vcat(nlp.pdemeta.Hobjkrows, nlp.pdemeta.Hobjrows .+ np, nlp.pdemeta.Hkrows, nlp.pdemeta.Hrows .+ np))
+  cols .= T.(vcat(nlp.pdemeta.Hobjkcols, nlp.pdemeta.Hobjcols .+ np, nlp.pdemeta.Hkcols, nlp.pdemeta.Hcols .+ np))
   return (rows, cols)
 end
 
@@ -301,10 +349,38 @@ function hess_coord!(
     end
     luh = split_res(xh, λf)
 
-    lag_hess = Gridap.FESpaces._hessian(x -> split_res(x, λf), xh, luh)
-    matdata = Gridap.FESpaces.collect_cell_matrix(lag_hess)
+    # lag_hess = Gridap.FESpaces.jacobian(Gridap.FESpaces._gradient(x -> split_res(x, λf), xh, luh), xh) # 
+    # lag_hess = Gridap.FESpaces._hessian(x -> split_res(x, λf), xh, luh)
+    lag_hess = _hessianv1(x -> split_res(x, λf), xh, luh)
+    matdata = Gridap.FESpaces.collect_cell_matrix(nlp.pdemeta.Y, nlp.pdemeta.X, lag_hess)
     assem = SparseMatrixAssembler(nlp.pdemeta.Y, nlp.pdemeta.X)
-    fill_hess_coo_numeric!(vals, assem, matdata, n = nini)
+    ##############################################################
+    # TO BE IMPROVED
+
+    m1 = Gridap.FESpaces.nz_counter(
+      Gridap.FESpaces.get_matrix_builder(assem),
+      (Gridap.FESpaces.get_rows(assem), Gridap.FESpaces.get_cols(assem)),
+    ) # Gridap.Algebra.CounterCS
+    Gridap.FESpaces.symbolic_loop_matrix!(m1, assem, matdata)
+    m2 = Gridap.FESpaces.nz_allocation(m1) # Gridap.Algebra.InserterCSC
+    Gridap.FESpaces.numeric_loop_matrix!(m2, assem, matdata)
+    m3 = sparse(LowerTriangular(Gridap.FESpaces.create_from_nz(m2)))
+    _, _, v = findnz(m3)
+    vals[(nini + 1):(nini + length(v))] .= v
+    nini += length(v)
+
+    # nini = fill_hess_coo_numeric!(vals, assem, matdata, n = nini)
+    ##############################################################
+    #=
+    # https://github.com/gridap/Gridap.jl/blob/master/src/FESpaces/SparseMatrixAssemblers.jl
+    nnzh = length(nlp.pdemeta.Hrows)
+    @show nlp.meta.nvar, nlp.pdemeta.nparam
+    ms = sparse(nlp.pdemeta.Hrows, nlp.pdemeta.Hcols, zeros(T, nnzh), nlp.meta.nvar - nlp.pdemeta.nparam, nlp.meta.nvar - nlp.pdemeta.nparam) # SparseArrays.sparse!
+    Gridap.FESpaces.numeric_loop_matrix!(ms, assem, matdata)
+    ms_lower = sparse(LowerTriangular(Gridap.FESpaces.create_from_nz(ms)))
+    vals[(nini + 1):(nini + nnzh)] .= ms_lower.nzval
+    nini += nnzh
+    =#
   end
 
   return vals
@@ -380,7 +456,9 @@ function jtprod!(nlp::GridapPDENLPModel, x::AbstractVector, v::AbstractVector, J
 end
 
 function jac_structure(nlp::GridapPDENLPModel)
-  return (nlp.pdemeta.Jrows, nlp.pdemeta.Jcols)
+  nparam = nlp.pdemeta.nparam
+  ny = num_free_dofs(nlp.pdemeta.Ypde)
+  return (vcat(nlp.pdemeta.Jkrows, nlp.pdemeta.Jyrows, nlp.pdemeta.Jurows), vcat(nlp.pdemeta.Jkcols, nlp.pdemeta.Jycols .+ nparam, nlp.pdemeta.Jucols .+ nparam .+ ny))
 end
 
 function jac_structure!(
@@ -389,8 +467,10 @@ function jac_structure!(
   cols::AbstractVector{T},
 ) where {T <: Integer}
   @lencheck nlp.meta.nnzj rows cols
-  rows .= T.(nlp.pdemeta.Jrows)
-  cols .= T.(nlp.pdemeta.Jcols)
+  nparam = nlp.pdemeta.nparam
+  ny = num_free_dofs(nlp.pdemeta.Ypde)
+  rows .= T.(vcat(nlp.pdemeta.Jkrows, nlp.pdemeta.Jyrows, nlp.pdemeta.Jurows))
+  cols .= T.(vcat(nlp.pdemeta.Jkcols, nlp.pdemeta.Jycols .+ nparam, nlp.pdemeta.Jucols .+ nparam .+ ny))
   return rows, cols
 end
 
@@ -407,6 +487,10 @@ function jac_coord!(nlp::GridapPDENLPModel, x::AbstractVector, vals::AbstractVec
     nlp.pdemeta.Xpde,
     nlp.pdemeta.Ycon,
     x,
+    nlp.pdemeta.Jyrows,
+    nlp.pdemeta.Jycols,
+    nlp.pdemeta.Jurows,
+    nlp.pdemeta.Jucols,
     vals,
     nlp.workspace.c,
     nlp.workspace.Jk,

src/PDENLPModels.jl

@@ -1,7 +1,6 @@
 module PDENLPModels
 
 #This package contains a list of NLPModels for PDE optimization.
-using LinearAlgebra: length
 using ForwardDiff, LinearAlgebra, SparseArrays, FastClosures
 
 #JSO packages
@@ -39,6 +38,9 @@ function Gridap.FESpaces.scatter_free_and_dirichlet_values(
   lazy_map(Broadcasting(Gridap.Arrays.PosNegReindex(free_values, dirichlet_values)), cell_dof_ids)
 end
 
+include("gridap_autodiff.jl")
+include("gridap_utils.jl")
+
 #Regroup the different types of FEFunction
 const FEFunctionType =
   Union{Gridap.FESpaces.SingleFieldFEFunction, Gridap.MultiField.MultiFieldFEFunction}

src/additional_constructors.jl

@@ -20,7 +20,8 @@ function GridapPDENLPModel(
 
   @assert nparam ≥ 0 throw(DimensionError("x0", nvar_pde, nvar))
 
-  rows, cols, nnzh = _compute_hess_structure(tnrj, Y, X, x0, nparam)
+  robj, rck, cobj, cck, no, nk = _compute_hess_structure(tnrj, Y, X, x0, nparam)
+  nnzh = nk + no
 
   if NRJ <: NoFETerm && typeof(lvar) <: AbstractVector && typeof(uvar) <: AbstractVector
     lv, uv = lvar, uvar
@@ -54,8 +55,18 @@ function GridapPDENLPModel(
     nvar_con,
     nparam,
     nnzh,
-    rows,
-    cols,
+    robj,
+    cobj,
+    rck, 
+    cck,
+    Int[],
+    Int[],
+    Int[],
+    Int[],
+    Int[],
+    Int[],
+    Int[],
+    Int[],
     Int[],
     Int[],
   )
@@ -169,23 +180,7 @@ function GridapPDENLPModel(
 
   @assert nparam >= 0 throw(DimensionError("x0", nvar_pde + nvar_con, nvar))
 
-  if !(typeof(Xcon) <: VoidFESpace) && !(typeof(Ycon) <: VoidFESpace)
-    _xpde = _fespace_to_multifieldfespace(Xpde)
-    _xcon = _fespace_to_multifieldfespace(Xcon)
-    #Handle the case where Ypde or Ycon are single field FE space(s).
-    _ypde = _fespace_to_multifieldfespace(Ypde)
-    _ycon = _fespace_to_multifieldfespace(Ycon)
-    #Build Y (resp. X) the trial (resp. test) space of the Multi Field function [y,u]
-    X = MultiFieldFESpace(vcat(_xpde.spaces, _xcon.spaces))
-    Y = MultiFieldFESpace(vcat(_ypde.spaces, _ycon.spaces))
-  elseif (typeof(Xcon) <: VoidFESpace) ⊻ (typeof(Ycon) <: VoidFESpace)
-    throw(ErrorException("Error: Xcon or Ycon are both nothing or must be specified."))
-  else
-    #_xpde = _fespace_to_multifieldfespace(Xpde)
-    X = Xpde #_xpde
-    #_ypde = _fespace_to_multifieldfespace(Ypde)
-    Y = Ypde #_ypde
-  end
+  Y, X = _to_multifieldfespace(Ypde, Xpde, Ycon, Xcon)
 
   if NRJ == NoFETerm && typeof(lvary) <: AbstractVector && typeof(uvary) <: AbstractVector
     lvar, uvar = vcat(lvary, lvaru, lvark), vcat(uvary, uvaru, uvark)
@@ -203,15 +198,17 @@ function GridapPDENLPModel(
   @lencheck nvar lvar uvar
   @lencheck ncon ucon y0
 
-  rows, cols, nnzh = _compute_hess_structure(tnrj, c, Y, Ypde, Ycon, X, Xpde, x0, nparam)
-  _, _, nnzh_obj = _compute_hess_structure(tnrj, Y, X, x0, nparam)
+  rk, ro, rck, rc, ck, co, cck, cc, nk, no, nck, nc = _compute_hess_structure(tnrj, c, Y, Ypde, Ycon, X, Xpde, x0, nparam)
+  nnzh_obj = nk + no
+  nnzh = nnzh_obj + nck + nc
+  #_, _, nnzh_obj = _compute_hess_structure(tnrj, Y, X, x0, nparam)
 
   if typeof(c) <: AffineFEOperator #Here we expect ncon = nvar_pde
     lin = 1:ncon
   end
   Jkrows, Jkcols, nnz_jac_k = jac_k_structure(nparam, ncon)
-  Jrows, Jcols, nini = _jacobian_struct(c, x0, Y, Xpde, Ypde, Ycon)
-  nnzj = nini + nnz_jac_k
+  Jyrows, Jurows, Jycols, Jucols, nini_y, nini_u = _jacobian_struct(c, x0, Y, Xpde, Ypde, Ycon)
+  nnzj = nini_y + nini_u + nnz_jac_k
 
   meta = NLPModelMeta{T, S}(
     nvar,
@@ -243,10 +240,20 @@ function GridapPDENLPModel(
     nvar_con,
     nparam,
     nnzh_obj,
-    rows,
-    cols,
-    vcat(Jkrows, Jrows),
-    vcat(Jkcols, Jcols .+ nparam),
+    rk,
+    ck,
+    ro,
+    co,
+    rck, 
+    cck, 
+    rc, 
+    cc,
+    Jkrows,
+    Jkcols,
+    Jyrows,
+    Jycols,
+    Jurows,
+    Jucols,
   )
 
   workspace = PDEWorkspace(T, S, nvar, ncon, nparam, nnzh, nnzj)

src/additional_obj_terms.jl

@@ -1,3 +1,11 @@
+trial_check(::TrialFESpace) = true
+trial_check(::Gridap.FESpaces.DirichletFESpace) = true
+trial_check(::Gridap.FESpaces.FESpaceWithConstantFixed) = true
+trial_check(::Gridap.FESpaces.ZeroMeanFESpace) = true
+trial_check(::Gridap.FESpaces.UnconstrainedFESpace) = true # When Ycon = TrialFESpace(Xcon), we have typeof(Xcon) == typeof(Ycon)
+trial_check(::VoidFESpace) = true
+trial_check(Y::MultiFieldFESpace) = all(trial_check, Y)
+
 abstract type AbstractEnergyTerm end
 
 """
@@ -128,6 +136,12 @@ struct EnergyFETerm <: AbstractEnergyTerm
   trian::Triangulation
   Ypde::FESpace
   Ycon::FESpace
+
+  function EnergyFETerm(f, trian, Ypde, Ycon)
+    @assert trial_check(Ypde)
+    @assert trial_check(Ycon)
+    return new(f, trian, Ypde, Ycon)
+  end
 end
 
 function EnergyFETerm(f::Function, trian::Triangulation, Ypde)
@@ -153,14 +167,16 @@ function _compute_gradient!(
   X::FESpace,
 )
   @lencheck 0 κ
-  cell_yu = Gridap.FESpaces.get_cell_dof_values(yu)
-  cell_id_yu = Gridap.Arrays.IdentityVector(length(cell_yu))
-  cell_r_yu = get_array(gradient(x -> _obj_integral(tnrj, κ, x), yu))
-  #Put the result in the format expected by Gridap.FESpaces.assemble_vector!
-  vecdata_yu = [[cell_r_yu], [cell_id_yu]] #TODO would replace by Tuple work?
+  fsplit(x) = if typeof(tnrj.Ycon) <: VoidFESpace
+    return tnrj.f(x)
+  else
+    y, u = _split_FEFunction(x, tnrj.Ypde, tnrj.Ycon)
+    return tnrj.f(y, u)
+  end
+  vec = Gridap.FESpaces.collect_cell_vector(X, gradient(fsplit, yu))
   #Assemble the gradient in the "good" space
   assem = Gridap.FESpaces.SparseMatrixAssembler(Y, X)
-  Gridap.FESpaces.assemble_vector!(g, assem, vecdata_yu)
+  Gridap.FESpaces.assemble_vector!(g, assem, vec)
   return g
 end
 
@@ -216,6 +232,8 @@ struct MixedEnergyFETerm <: AbstractEnergyTerm
     Ycon::FESpace,
   )
     @assert n > 0
+    @assert trial_check(Ypde)
+    @assert trial_check(Ycon)
     return new(f, trian, n, inde, Ypde, Ycon)
   end
 end
@@ -255,18 +273,19 @@ function _compute_gradient!(
   nyu = num_free_dofs(Y)
   @lencheck tnrj.nparam + nyu g
 
-  cell_yu = Gridap.FESpaces.get_cell_dof_values(yu)
-  cell_id_yu = Gridap.Arrays.IdentityVector(length(cell_yu))
-
-  cell_r_yu = get_array(gradient(x -> _obj_integral(tnrj, κ, x), yu))
-  #Put the result in the format expected by Gridap.FESpaces.assemble_vector
-  vecdata_yu = [[cell_r_yu], [cell_id_yu]] #TODO would replace by Tuple work?
+  fsplit(x) = if typeof(tnrj.Ycon) <: VoidFESpace
+    return tnrj.f(κ, x)
+  else
+    y, u = _split_FEFunction(x, tnrj.Ypde, tnrj.Ycon)
+    return tnrj.f(κ, y, u)
+  end
+  vec = Gridap.FESpaces.collect_cell_vector(X, gradient(fsplit, yu))
   #Assemble the gradient in the "good" space
   assem = Gridap.FESpaces.SparseMatrixAssembler(Y, X)
   Gridap.FESpaces.assemble_vector!(
     view(g, (tnrj.nparam + 1):(tnrj.nparam + nyu)),
     assem,
-    vecdata_yu,
+    vec,
   )
 
   _compute_gradient_k!(view(g, 1:(tnrj.nparam)), tnrj, κ, yu)