Minor

src/MultilevelTools/DistributedGridTransferOperators.jl

@@ -88,7 +88,7 @@ function _get_projection_cache(lev::Int,sh::FESpaceHierarchy,qdegree::Int,mode::
 
     model_H = get_model(mh,lev+1)
     UH   = get_fe_space(sh,lev+1)
-    VH   = get_test_space(UH)
+    VH   = get_fe_space(sh,lev+1)
     ΩH   = Triangulation(model_H)
     dΩH  = Measure(ΩH,qdegree)
     dΩhH = Measure(ΩH,Ωh,qdegree)
@@ -106,7 +106,7 @@ function _get_projection_cache(lev::Int,sh::FESpaceHierarchy,qdegree::Int,mode::
     u,v    = get_trial_fe_basis(UH), get_fe_basis(VH)
     data   = collect_cell_matrix_and_vector(UH,VH,aH(u,v),lH(v,u00),u_dir)
     AH,bH0 = assemble_matrix_and_vector(assem,data)
-    xH     = pfill(0.0,partition(axes(AH,2)))
+    xH     = allocate_in_domain(AH)
     bH     = copy(bH0)
 
     cache_refine = model_h, Uh, fv_h, dv_h, VH, AH, lH, xH, bH, bH0, assem

src/MultilevelTools/FESpaceHierarchies.jl

@@ -161,11 +161,19 @@ end
 
 # Computing system matrices
 
-function compute_hierarchy_matrices(trials::FESpaceHierarchy,a::Function,l::Function,qdegree::Integer)
-  return compute_hierarchy_matrices(trials,a,l,Fill(qdegree,num_levels(trials)))
-end
-
-function compute_hierarchy_matrices(trials::FESpaceHierarchy,a::Function,l::Function,qdegree::AbstractArray{<:Integer})
+function compute_hierarchy_matrices(trials::FESpaceHierarchy,
+                                    tests::FESpaceHierarchy,
+                                    a::Function,
+                                    l::Function,
+                                    qdegree::Integer)
+  return compute_hierarchy_matrices(trials,tests,a,l,Fill(qdegree,num_levels(trials)))
+end
+
+function compute_hierarchy_matrices(trials::FESpaceHierarchy,
+                                    tests::FESpaceHierarchy,
+                                    a::Function,
+                                    l::Function,
+                                    qdegree::AbstractArray{<:Integer})
   nlevs = num_levels(trials)
   mh    = trials.mh
 
@@ -179,7 +187,7 @@ function compute_hierarchy_matrices(trials::FESpaceHierarchy,a::Function,l::Func
     if i_am_in(parts)
       model = get_model(mh,lev)
       U = get_fe_space(trials,lev)
-      V = get_test_space(U)
+      V = get_fe_space(tests,lev)
       Ω = Triangulation(model)
       dΩ = Measure(Ω,qdegree[lev])
       ai(u,v) = a(u,v,dΩ)
@@ -195,14 +203,3 @@ function compute_hierarchy_matrices(trials::FESpaceHierarchy,a::Function,l::Func
   end
   return mats, A, b
 end
-
-function get_test_space(U::GridapDistributed.DistributedSingleFieldFESpace)
-  spaces = map(local_views(U)) do U
-    if isa(U,Gridap.FESpaces.UnconstrainedFESpace)
-      U
-    else
-      U.space
-    end
-  end
-  return GridapDistributed.DistributedSingleFieldFESpace(spaces,U.gids,U.vector_type)
-end

src/PatchBasedSmoothers/seq/PatchBasedLinearSolvers.jl

@@ -1,9 +1,10 @@
 
-struct PatchBasedLinearSolver{A,B} <: Gridap.Algebra.LinearSolver
+struct PatchBasedLinearSolver{A,B,C,D} <: Gridap.Algebra.LinearSolver
   bilinear_form  :: Function
-  Ph             :: A
-  Vh             :: B
-  local_solver   :: Gridap.Algebra.LinearSolver
+  patch_space    :: A
+  space          :: B
+  measure        :: C
+  local_solver   :: D    
 end
 
 struct PatchBasedSymbolicSetup <: Gridap.Algebra.SymbolicSetup
@@ -16,17 +17,19 @@ end
 
 struct PatchBasedSmootherNumericalSetup{A,B,C} <: Gridap.Algebra.NumericalSetup
   solver   :: PatchBasedLinearSolver
-  Ap_ns    :: A
+  local_ns :: A
   weights  :: B
   caches   :: C
 end
 
 function Gridap.Algebra.numerical_setup(ss::PatchBasedSymbolicSetup,A::AbstractMatrix)
-  Ph, Vh, solver = ss.solver.Ph, ss.solver.Vh, ss.solver
+  solver = ss.solver
+  Ph, Vh, dΩ, solver = solver.patch_space, solver.space, solver.measure
   weights = compute_weight_operators(Ph,Vh)
 
   assembler = SparseMatrixAssembler(Ph,Ph)
-  Ap        = assemble_matrix(solver.bilinear_form,assembler,Ph,Ph)
+  ap(u,v)   = solver.bilinear_form(u,v,dΩ)
+  Ap        = assemble_matrix(ap,assembler,Ph,Ph)
   Ap_ns     = numerical_setup(symbolic_setup(solver.local_solver,Ap),Ap)
 
   # Caches
@@ -38,17 +41,15 @@ function Gridap.Algebra.numerical_setup(ss::PatchBasedSymbolicSetup,A::AbstractM
 end
 
 function Gridap.Algebra.numerical_setup(ss::PatchBasedSymbolicSetup,A::PSparseMatrix)
-  Ph, Vh, solver = ss.solver.Ph, ss.solver.Vh, ss.solver
-  weights = compute_weight_operators(Ph,Vh)
-
-  # Patch system solver
-  # Only local systems need to be solved
-  u = get_trial_fe_basis(Ph)
-  v = get_fe_basis(Ph)
-  matdata = collect_cell_matrix(Ph,Ph,solver.bilinear_form(u,v))
-  Ap_ns = map(local_views(Ph),matdata) do Ph, matdata
-    assemb = SparseMatrixAssembler(Ph,Ph)
-    Ap     = assemble_matrix(assemb,matdata)
+  solver = ss.solver
+  Ph, Vh, dΩ = solver.patch_space, solver.space, solver.measure
+  #weights = compute_weight_operators(Ph,Vh)
+
+  # Patch system solver (only local systems need to be solved)
+  Ap_ns = map(local_views(Ph),local_views(dΩ)) do Ph, dΩ
+    assembler = SparseMatrixAssembler(Ph,Ph)
+    ap(u,v)   = solver.bilinear_form(u,v,dΩ)
+    Ap        = assemble_matrix(ap,assembler,Ph,Ph)
     return numerical_setup(symbolic_setup(solver.local_solver,Ap),Ap)
   end
 
@@ -59,15 +60,15 @@ function Gridap.Algebra.numerical_setup(ss::PatchBasedSymbolicSetup,A::PSparseMa
   x      = pfill(0.0,partition(Vh.gids))
   caches = (rp,dxp,r,x)
 
-  return PatchBasedSmootherNumericalSetup(solver,Ap_ns,weights,caches)
+  return PatchBasedSmootherNumericalSetup(solver,Ap_ns,nothing,caches)
 end
 
 function Gridap.Algebra.numerical_setup!(ns::PatchBasedSmootherNumericalSetup, A::AbstractMatrix)
   @notimplemented
 end
 
 function Gridap.Algebra.solve!(x::AbstractVector,ns::PatchBasedSmootherNumericalSetup,r::AbstractVector)
-  Ap_ns, weights, caches = ns.Ap_ns, ns.weights, ns.caches
+  Ap_ns, weights, caches = ns.local_ns, ns.weights, ns.caches
 
   Ph = ns.solver.Ph
   w, w_sums = weights
@@ -81,9 +82,9 @@ function Gridap.Algebra.solve!(x::AbstractVector,ns::PatchBasedSmootherNumerical
 end
 
 function Gridap.Algebra.solve!(x_mat::PVector,ns::PatchBasedSmootherNumericalSetup,r_mat::PVector)
-  Ap_ns, weights, caches = ns.Ap_ns, ns.weights, ns.caches
+  Ap_ns, weights, caches = ns.local_ns, ns.weights, ns.caches
 
-  Ph = ns.solver.Ph
+  Ph = ns.solver.patch_space
   w, w_sums = weights
   rp, dxp, r, x = caches
 

src/PatchBasedSmoothers/seq/PatchMultiFieldFESpaces.jl

@@ -5,6 +5,8 @@ struct PatchDistributedMultiFieldFESpace{A,B}
   gids   :: B
 end
 
+GridapDistributed.local_views(a::PatchDistributedMultiFieldFESpace) = a.spaces
+
 ## PatchFESpace from MultiFieldFESpace
 
 function PatchFESpace(space::Gridap.MultiField.MultiFieldFESpace,

test/_dev/GMG/GMG_Multifield.jl

@@ -48,10 +48,7 @@ function compute_matrices(trials,tests,a::Function,l::Function,qdegree)
   return mats, A, b
 end
 
-function get_patch_smoothers(tests,patch_spaces,patch_decompositions,biform,qdegree)
-  tests_u, tests_j = tests
-  patch_spaces_u, patch_spaces_j = patch_spaces
-
+function get_patch_smoothers(tests,patch_spaces,patch_decompositions,biform,qdegree)  
   mh = tests.mh
   nlevs = num_levels(mh)
   smoothers = Vector{RichardsonSmoother}(undef,nlevs-1)
@@ -63,9 +60,8 @@ function get_patch_smoothers(tests,patch_spaces,patch_decompositions,biform,qdeg
       Vh = get_fe_space(tests,lev)
       Ω  = Triangulation(PD)
       dΩ = Measure(Ω,qdegree)
-      a(u,v) = biform(u,v,dΩ)
-      local_solver = PETScLinearSolver(set_ksp_options) # IS_ConjugateGradientSolver(;reltol=1.e-6)
-      patch_smoother = PatchBasedLinearSolver(a,Ph,Vh,local_solver)
+      local_solver = LUSolver() # IS_ConjugateGradientSolver(;reltol=1.e-6)
+      patch_smoother = PatchBasedLinearSolver(biform,Ph,Vh,dΩ,local_solver)
       smoothers[lev] = RichardsonSmoother(patch_smoother,10,0.2)
     end
   end
@@ -74,7 +70,7 @@ end
 
 np       = 1    # Number of processors
 D        = 3    # Problem dimension
-n_refs_c = 2    # Number of refinements for the coarse model
+n_refs_c = 1    # Number of refinements for the coarse model
 n_levels = 2    # Number of refinement levels
 order    = 1    # FE order
 
@@ -116,32 +112,34 @@ smatrices, A, b = compute_matrices(trials,tests,biform,liform,qdegree);
 pbs = GridapSolvers.PatchBasedSmoothers.PatchBoundaryExclude()
 patch_decompositions = PatchDecomposition(mh;patch_boundary_style=pbs)
 patch_spaces = PatchFESpace(tests,patch_decompositions);
-
 smoothers = get_patch_smoothers(tests,patch_spaces,patch_decompositions,biform,qdegree)
 
+smoother_ns = numerical_setup(symbolic_setup(smoothers[1],A),A)
+
 restrictions, prolongations = setup_transfer_operators(trials,qdegree;mode=:residual);
 
-GridapPETSc.with() do
-  gmg = GMGLinearSolver(mh,
-                        smatrices,
-                        prolongations,
-                        restrictions,
-                        pre_smoothers=smoothers,
-                        post_smoothers=smoothers,
-                        coarsest_solver=PETScLinearSolver(set_ksp_options),
-                        maxiter=1,
-                        rtol=1.0e-10,
-                        verbose=false,
-                        mode=:preconditioner)
-
-  solver = CGSolver(gmg;maxiter=100,atol=1e-10,rtol=1.e-6,verbose=i_am_main(ranks))
-  ns = numerical_setup(symbolic_setup(solver,A),A)
-
-  x = pfill(0.0,partition(axes(A,2)))
-  solve!(x,ns,b)
-  @time begin
-    fill!(x,0.0)
-    solve!(x,ns,b)
-  end
-  println("n_dofs = ", length(x))
-end
+
+#GridapPETSc.with() do
+#  gmg = GMGLinearSolver(mh,
+#                        smatrices,
+#                        prolongations,
+#                        restrictions,
+#                        pre_smoothers=smoothers,
+#                        post_smoothers=smoothers,
+#                        coarsest_solver=PETScLinearSolver(set_ksp_options),
+#                        maxiter=1,
+#                        rtol=1.0e-10,
+#                        verbose=false,
+#                        mode=:preconditioner)
+#
+#  solver = CGSolver(gmg;maxiter=100,atol=1e-10,rtol=1.e-6,verbose=i_am_main(ranks))
+#  ns = numerical_setup(symbolic_setup(solver,A),A)
+#
+#  x = pfill(0.0,partition(axes(A,2)))
+#  solve!(x,ns,b)
+#  @time begin
+#    fill!(x,0.0)
+#    solve!(x,ns,b)
+#  end
+#  println("n_dofs = ", length(x))
+#end