Misc bugfixes

ext/GridapPETScExt/HipmairXuSolvers.jl

@@ -12,9 +12,10 @@ function ADS_Solver(model,tags,order;rtol=1e-8,maxits=300)
 end
 
 function get_ads_spaces(model,order,tags)
+  D = num_cell_dims(model)
   reffe_H1_sc = ReferenceFE(lagrangian,Float64,order)
   V_H1_sc = FESpace(model,reffe_H1_sc;dirichlet_tags=tags)
-  
+
   reffe_H1 = ReferenceFE(lagrangian,VectorValue{D,Float64},order)
   V_H1 = FESpace(model,reffe_H1;dirichlet_tags=tags)
 

src/BlockSolvers/BlockDiagonalSolvers.jl

@@ -118,7 +118,7 @@ end
 
 function Gridap.Algebra.symbolic_setup(solver::BlockDiagonalSolver,mat::AbstractBlockMatrix,x::AbstractBlockVector)
   mat_blocks   = diag(blocks(mat))
-  block_ss     = map((b,m) -> block_symbolic_setup(b,m,x),solver.blocks,solver.solvers,mat_blocks)
+  block_ss     = map((b,s,m) -> block_symbolic_setup(b,s,m,x),solver.blocks,solver.solvers,mat_blocks)
   return BlockDiagonalSolverSS(solver,block_ss)
 end
 

src/BlockSolvers/BlockFEOperators.jl

@@ -56,7 +56,7 @@ function BlockFEOperator(
   global_op  = FEOperator(global_res,global_jac,trial,test,assem)
 
   block_ops = map(FEOperator,res,jac,blocks(trial),blocks(test),blocks(assem))
-  return BlockFEOperator{NB,SB,P}(global_op,block_ids,block_ops,nonlinear)
+  return BlockFEOperator{NB,SB,P}(global_op,ids,block_ops,nonlinear)
 end
 
 # BlockArrays API
@@ -94,7 +94,7 @@ function liform_from_blocks(ids, ranges, liforms)
     c = DomainContribution()
     for (I,lf) in zip(ids, liforms)
       vk = v[ranges[I]]
-      c += lf(uk,vk)
+      c += lf(vk)
     end
     return c
   end

src/LinearSolvers/RichardsonLinearSolvers.jl

@@ -34,11 +34,9 @@ function get_solver_caches(solver::RichardsonLinearSolver, A::AbstractMatrix)
   ω = solver.ω
   z = allocate_in_domain(A)
   r = allocate_in_domain(A)
-  α = allocate_in_domain(A)
   fill!(z,0.0)
   fill!(r,0.0)
-  fill!(α,1.0)
-  return ω, z, r, α
+  return ω, z, r
 end
 
 mutable struct RichardsonLinearNumericalSetup<:Gridap.Algebra.NumericalSetup
@@ -80,26 +78,24 @@ end
 
 function Gridap.Algebra.solve!(x::AbstractVector, ns:: RichardsonLinearNumericalSetup, b::AbstractVector)
   solver,A,Pl,caches = ns.solver,ns.A,ns.Pl_ns,ns.caches
-  ω, z, r, α = caches
+  ω, z, r = caches
   log = solver.log
-  # Relaxation parameters 
-  α .*= ω
   # residual
   r .= b
   mul!(r, A, x, -1, 1)
   done = init!(log,norm(r))
-  if !isa(ns.Pl_ns,Nothing) # Case when a preconditioner is applied
+  if !isnothing(ns.Pl_ns) # Case when a preconditioner is applied
     while !done
-      solve!(z, Pl, r) # Apply preconditioner r = PZ 
-      x .+= α.* z
+      solve!(z, Pl, r) # Apply preconditioner r = PZ
+      x .+= ω .* z
       r .= b
       mul!(r, A, x, -1, 1)
       done = update!(log,norm(r))
     end
     finalize!(log,norm(r))
   else    # Case when no preconditioner is applied
     while !done
-      x .+= α.* r
+      x .+= ω .* r
       r .= b
       mul!(r, A, x, -1, 1)
       done = update!(log,norm(r))

src/LinearSolvers/SchurComplementSolvers.jl

@@ -63,7 +63,7 @@ function Gridap.Algebra.solve!(x::AbstractBlockVector,ns::SchurComplementNumeric
 
   # Solve Schur complement
   solve!(x_u,A,y_u)                      # x_u = A^-1 y_u
-  copy!(bp,y_p); mul!(bp,C,du,-1.0,1.0)  # bp  = C*(A^-1 y_u) - y_p
+  copy!(bp,y_p); mul!(bp,C,x_u,-1.0,1.0)  # bp  = y_p - C*(A^-1 y_u)
   solve!(x_p,S,bp)                       # x_p = S^-1 bp
 
   mul!(bu,B,x_p)         # bu  = B*x_p

src/MultilevelTools/GridTransferOperators.jl

@@ -24,13 +24,13 @@ function RedistributionOperator(
     dof_ids_to = partition(get_free_dof_ids(space_to))
     to_cell_to_to_lid = _get_unified_cell_dof_ids(space_to)
   else
-    dof_ids_to, to_cell_to_to_lid = nothing, nothing, nothing
+    dof_ids_to, to_cell_to_to_lid = nothing, nothing
   end
   if !isnothing(space_from)
     dof_ids_from = partition(get_free_dof_ids(space_from))
     from_cell_to_from_lid = _get_unified_cell_dof_ids(space_from)
   else
-    dof_ids_from, from_cell_to_from_lid = nothing, nothing, nothing
+    dof_ids_from, from_cell_to_from_lid = nothing, nothing
   end
   indices_from, indices_to = GridapDistributed.redistribute_indices(
     dof_ids_from, from_cell_to_from_lid, to_cell_to_to_lid, model_to, glue; reverse,

src/MultilevelTools/HierarchicalArrays.jl

@@ -122,11 +122,9 @@ function Base.map!(f,a::HierarchicalArray,args::Vararg{HierarchicalArray,N}) whe
   @assert matching_level_parts(a,args...)
   ranks  = get_level_parts(a)
   arrays = map(a -> a.array, args)
-  map(ranks, a.array, arrays...) do ranks, ai, arrays...
+  map(eachindex(a.array), ranks, arrays...) do i, ranks, arrays...
     if i_am_in(ranks)
-      ai = f(arrays...)
-    else
-      nothing
+      a.array[i] = f(arrays...)
     end
   end
   return a

src/MultilevelTools/ModelHierarchies.jl

@@ -133,7 +133,7 @@ function CartesianModelHierarchy(
       model_ref  = Gridap.Adaptivity.refine(cmodel,level_nrefs[lev])
       ref_glue   = Gridap.Adaptivity.get_adaptivity_glue(model_ref)
     else
-      model_ref, ref_glue = nothing, nothing, nothing
+      model_ref, ref_glue = nothing, nothing
     end
     # Redistribution (if needed)
     if i_am_in(fparts) && (cparts !== fparts)

src/NonlinearSolvers/ContinuationFEOperators.jl

@@ -119,22 +119,42 @@ end
 
 function Algebra.allocate_residual(op::ContinuationFEOperator,u)
   switch!(op.switch, u, nothing)
-  ifelse(!has_switched(op), allocate_residual(op.op1, u), allocate_residual(op.op2, u))
+  if !has_switched(op)
+    allocate_residual(op.op1, u)
+  else
+    allocate_residual(op.op2, u)
+  end
 end
 
 function Algebra.residual!(b::AbstractVector,op::ContinuationFEOperator,u)
   switch!(op.switch, u, b)
-  ifelse(!has_switched(op), residual!(b, op.op1, u), residual!(b, op.op2, u))
+  if !has_switched(op)
+    residual!(b, op.op1, u)
+  else
+    residual!(b, op.op2, u)
+  end
 end
 
 function Algebra.allocate_jacobian(op::ContinuationFEOperator,u)
-  ifelse(!has_switched(op), allocate_jacobian(op.op1, u), allocate_jacobian(op.op2, u))
+  if !has_switched(op)
+    allocate_jacobian(op.op1, u)
+  else
+    allocate_jacobian(op.op2, u)
+  end
 end
 
 function Algebra.jacobian!(A::AbstractMatrix,op::ContinuationFEOperator,u)
   if op.reuse_caches
-    ifelse(!has_switched(op), jacobian!(A, op.op1, u), jacobian!(A, op.op2, u))
+    if !has_switched(op)
+      jacobian!(A, op.op1, u)
+    else
+      jacobian!(A, op.op2, u)
+    end
   else
-    ifelse(!has_switched(op), jacobian(op.op1, u), jacobian(op.op2, u))
+    if !has_switched(op)
+      jacobian(op.op1, u)
+    else
+      jacobian(op.op2, u)
+    end
   end
 end

src/PatchBasedSmoothers/PatchSolvers.jl

@@ -133,7 +133,7 @@ function Algebra.numerical_setup(ss::PatchSS{<:PatchSolverFromWeakform},mat::PSp
 end
 
 function Algebra.numerical_setup!(ns::PatchNS,mat::AbstractMatrix,vec::AbstractVector)
-  solver = nc.solver
+  solver = ns.solver
   @assert solver.is_nonlinear
   assem, trial, test = solver.assem, solver.trial, solver.test
 
@@ -197,7 +197,7 @@ function patch_solver_caches(
   patch_factorizations, caches = patch_solver_caches(
     patch_cols, patch_mats; collect_factorizations
   )
-  return patch_cols, patch_rows, patch_factorizations, caches
+  return patch_rows, patch_cols, patch_factorizations, caches
 end
 
 function patch_solver_caches(

src/SolverInterfaces/NullSpaces.jl

@@ -55,7 +55,7 @@ function is_orthogonal(N::NullSpace, v::AbstractVector; tol = 1.e-12)
 end
 
 function is_orthogonal(N::NullSpace, A::AbstractMatrix; tol = 1.e-12)
-  @assert length(v) == size(N,2)
+  @assert size(A,2) == size(N,2)
   v = allocate_in_range(A)
   for w in N.V
     mul!(v, A, w)

src/SolverInterfaces/PAExtras.jl

@@ -144,7 +144,7 @@ end
 function own_local_values_view(::Type{<:SparseMatrixCSR{Bi}},A::PSparseMatrix, new_cols) where Bi
   rows, cols = axes(A)
   map(partition(A),partition(rows),partition(cols),new_cols) do A, rows, cols, new_cols
-    m, rowptr, colval, nzval = A.n, A.rowptr, A.colval, A.nzval
+    m, rowptr, colval, nzval = A.m, A.rowptr, A.colval, A.nzval
 
     n_new = length(new_cols)
     old_to_new_cols = GridapDistributed.find_local_to_local_map(cols, new_cols)