diff --git a/.vscode/settings.json b/.vscode/settings.json
new file mode 100644
index 0000000..9e26dfe
--- /dev/null
+++ b/.vscode/settings.json
@@ -0,0 +1 @@
+{}
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
index 1540531..6ca6f6a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -20,7 +20,7 @@ Gridap = "0.17"
 GridapDistributed = "0.2"
 MPI = "0.14, 0.15, 0.16, 0.17, 0.18, 0.19"
 PETSc_jll = "3.13"
-PartitionedArrays = "0.2.4"
+PartitionedArrays = "0.2.13"
 SparseMatricesCSR = "0.6.6"
 julia = "1.3"
 
diff --git a/src/Environment.jl b/src/Environment.jl
index f642c75..882efc0 100644
--- a/src/Environment.jl
+++ b/src/Environment.jl
@@ -3,6 +3,12 @@ function Init(;args=String[],file="",help="",finalize_atexit=true)
   if !MPI.Initialized()
       MPI.Init()
   end
+  #To avoid multiple printing of the same line in parallel
+  # if MPI.Comm_rank(MPI.COMM_WORLD) != 0
+  #   redirect_stderr(devnull)
+  #   redirect_stdout(devnull)
+  # end 
+
   if finalize_atexit
     atexit(Finalize)
   end
diff --git a/src/GridapPETSc.jl b/src/GridapPETSc.jl
index defef13..d1e34c0 100644
--- a/src/GridapPETSc.jl
+++ b/src/GridapPETSc.jl
@@ -5,6 +5,12 @@ using Libdl
 using Gridap.Helpers
 using Gridap.Algebra
 using Gridap.Arrays
+using Gridap.FESpaces
+using Gridap.MultiField
+
+using GridapDistributed
+using GridapDistributed.MultiField
+
 using LinearAlgebra
 using SparseArrays
 using SparseMatricesCSR
@@ -69,19 +75,25 @@ end
 include("PETSC.jl")
 
 using GridapPETSc.PETSC: @check_error_code
-using GridapPETSc.PETSC: PetscBool, PetscInt, PetscScalar, Vec, Mat, KSP, PC, SNES
+using GridapPETSc.PETSC: PetscBool, PetscInt, PetscScalar, Vec, Mat, KSP, PC, SNES, IS
 #export PETSC
 export @check_error_code
-export PetscBool, PetscInt, PetscScalar, Vec, Mat, KSP, PC
+export PetscBool, PetscInt, PetscScalar, Vec, Mat, KSP, PC, IS
 
 include("Environment.jl")
 
 export PETScVector
 export PETScMatrix
 export petsc_sparse
+export PETScIS
 include("PETScArrays.jl")
+include("PETScIndexes.jl")
 include("PartitionedArrays.jl")
 
+
+export PETScFieldSplit
+include("PETScFieldSplits.jl")
+
 export PETScLinearSolver
 include("PETScLinearSolvers.jl")
 
diff --git a/src/PETSC.jl b/src/PETSC.jl
index a8623d0..a0409a0 100644
--- a/src/PETSC.jl
+++ b/src/PETSC.jl
@@ -648,6 +648,9 @@ Base.unsafe_convert(::Type{Ptr{Cvoid}},v::PC) = v.ptr
 @wrapper(:PCFactorSetUpMatSolverType,PetscErrorCode,(PC,),(pc,),"https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCFactorSetUpMatSolverType.html")
 @wrapper(:PCFactorGetMatrix,PetscErrorCode,(PC,Ptr{Mat}),(ksp,mat),"https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCFactorGetMatrix.html")
 
+#PCType() = PCType(Ptr{Cvoid}())
+@wrapper(:PCGetType,PetscErrorCode,(PC, Ptr{Ptr{Cstring}}),(pc,typ),"https://petsc.org/main/docs/manualpages/PC/PCGetType/")
+
 
 """
 Julia alias for the `SNES` C type.
@@ -703,4 +706,88 @@ const SNESPATCH            = "patch"
 @wrapper(:PetscObjectRegisterDestroy,PetscErrorCode,(Ptr{Cvoid},),(obj,),"https://petsc.org/release/docs/manualpages/Sys/PetscObjectRegisterDestroy.html")
 @wrapper(:PetscObjectRegisterDestroyAll,PetscErrorCode,(),(),"https://petsc.org/release/docs/manualpages/Sys/PetscObjectRegisterDestroyAll.html")
 
+"""
+ Julia alias for the `IS` C type.
+ See [PETSc manual](https://petsc.org/release/docs/manualpages/IS/IS/)).
+ """
+ struct IS
+   ptr::Ptr{Cvoid}
+ end
+ IS() = IS(Ptr{Cvoid}())
+ Base.convert(::Type{IS},p::Ptr{Cvoid}) = IS(p)
+ Base.unsafe_convert(::Type{Ptr{Cvoid}},v::IS) = v.ptr
+
+
+"""
+Julia alias for `ISType` C type.
+
+See [PETSc manual](https://petsc.org/main/docs/manualpages/IS/ISType/).
+"""
+const ISType = Cstring
+
+
+ """
+ Julia alias for `PetscCopyMode` C type.
+ See [PETSc manual](https://petsc.org/main/docs/manualpages/Sys/PetscCopyMode/).
+ """
+ @enum PetscCopyMode begin
+   PETSC_COPY_VALUES
+   PETSC_OWN_POINTER
+   PETSC_USE_POINTER
+ end
+@wrapper(:PetscObjectSetName, PetscErrorCode, (Ptr{Cvoid}, Cstring), (field, name), "https://petsc.org/main/docs/manualpages/Sys/PetscObjectSetName/") 
+@wrapper(:ISCreateGeneral,PetscErrorCode,(MPI.Comm, PetscInt, Ptr{PetscInt}, PetscCopyMode, Ptr{IS}),(comm, n, idx, mode, is), "https://petsc.org/main/docs/manualpages/IS/ISCreateGeneral/")
+@wrapper(:ISView,PetscErrorCode,(IS, PetscViewer),(is, viewer), "https://petsc.org/main/docs/manualpages/IS/ISView/")
+@wrapper(:ISCreateBlock,PetscErrorCode,(MPI.Comm, PetscInt, PetscInt, Ptr{PetscInt}, PetscCopyMode, Ptr{IS}),(comm, bs, n, idx, mode, is), "https://petsc.org/release/docs/manualpages/IS/ISCreateBlock/")
+@wrapper(:ISGetIndices,PetscErrorCode,(IS, Ptr{PetscInt}),(is, ptr), "https://petsc.org/release/docs/manualpages/IS/ISGetIndices/")
+@wrapper(:ISExpand,PetscErrorCode,(IS, IS, Ptr{IS}),(is1, is2, isout), "https://petsc.org/release/docs/manualpages/IS/ISGetIndices/")
+@wrapper(:ISGetSize,PetscErrorCode,(IS, Ptr{PetscInt}), (is, nsize), "https://petsc.org/main/docs/manualpages/IS/ISGetSize/")
+@wrapper(:ISDestroy, PetscErrorCode,(Ptr{IS},), (pis,), "https://petsc.org/main/docs/manualpages/IS/ISDestroy/")
+
+@wrapper(:PCFieldSplitSetIS,PetscErrorCode,(PC, Cstring, IS),(pc, Cfieldname, is), "https://petsc.org/release/docs/manualpages/PC/PCFieldSplitSetIS/")
+
+"""
+Julia alias for `PCCompositeType` C type.
+See [PETSc manual](https://petsc.org/release/docs/manualpages/PC/PCCompositeType/).
+"""
+@enum PCCompositeType begin
+  PC_COMPOSITE_ADDITIVE
+  PC_COMPOSITE_MULTIPLICATIVE
+  PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE
+  PC_COMPOSITE_SPECIAL
+  PC_COMPOSITE_SCHUR
+  PC_COMPOSITE_GKB
+end
+
+
+@wrapper(:PCFieldSplitSetType,PetscErrorCode,(PC, PCCompositeType),(pc, pcctype), "https://petsc.org/release/docs/manualpages/PC/PCFieldSplitSetType/")
+
+
+#PETSc Print
+@wrapper(:PetscPrintf,PetscErrorCode,(MPI.Comm, Cstring ),(comm, args ),"https://petsc.org/release/docs/manualpages/Sys/PetscPrintf/")
+@wrapper(:PetscSynchronizedPrintf,PetscErrorCode,(MPI.Comm, Cstring),(comm, args),"https://petsc.org/main/docs/manualpages/Sys/PetscSynchronizedPrintf/")
+
+
+#PETSc Sleep
+@wrapper(:PetscSleep,PetscErrorCode,(PetscReal,),(s,),"https://petsc.org/main/docs/manualpages/Sys/PetscSleep/")
+
+
+#PETSc Alloc
+@wrapper(:PetscMallocA,PetscErrorCode,(PetscInt, PetscBool, PetscInt, Ptr{Cstring}, Ptr{Cstring}, Csize_t, Ptr{Cvoid},), (n, clear,  lineno, fun, fname, bytes0, ptr0,), "https://petsc.org/release/docs/manualpages/Sys/PetscMallocA/#petscmalloca")
+#PETSC_EXTERN PetscErrorCode PetscMallocA(int,PetscBool,int,const char *,const char *,size_t,void *,...);
+
+"""
+    @PetscMalloc1
+
+See [PETSc manual](https://petsc.org/release/docs/manualpages/Sys/PetscMalloc1/).
+"""
+macro PetscMalloc1(m1,r1)
+  quote
+    func = Ptr{Nothing}()
+    linec = Ptr{Nothing}()
+    PetscMallocA(1,PETSC_FALSE, Cint(11), Ref(Cstring(linec)), Ref(Cstring(func)), (Csize_t)($m1)*sizeof($r1),($r1))
+  end
+end
+# PetscMalloc1(m1,r1) PetscMallocA(1,PETSC_FALSE,__LINE__,PETSC_FUNCTION_NAME,__FILE__,(size_t)(m1)*sizeof(**(r1)),(r1))
+
 end # module
diff --git a/src/PETScFieldSplits.jl b/src/PETScFieldSplits.jl
new file mode 100644
index 0000000..9cdf9cc
--- /dev/null
+++ b/src/PETScFieldSplits.jl
@@ -0,0 +1,102 @@
+# PETScFieldSplit 
+# MF (multifield) is a vector where each element is a SingleFieldFESpace or DistributedSingleFieldFESpace, eg: [U,P]
+# tags is a vector containing the names of the field (in order), eg: ["vel", "pres"]
+# show_idx is Bool vector, if true it prints the indexes
+struct PETScFieldSplit
+  MF
+  tags::Vector{String}
+  show_idx::Bool
+end
+
+# PETScFieldSplit constructor from Multifield
+function PETScFieldSplit(MF::Union{MultiFieldFESpace,GridapDistributed.DistributedMultiFieldFESpace}, tags::Vector{String}; show_idx=false)
+  n_fields = length(MF)
+  n_tags = length(tags)
+  @assert n_tags == n_fields #Verify that at each field a name is assigned
+  M = [MF[1]] #Vector of SingleFieldFESpace
+  for i = 2:1:n_fields
+    M = [M..., MF[i]]
+  end
+  PETScFieldSplit(M, tags, show_idx)
+end
+
+function PETScFieldSplit(MF::Union{Vector{SingleFieldFESpace}, Vector{<:GridapDistributed.DistributedSingleFieldFESpace}}, tags; show_idx=false)
+  n_fields = length(MF)
+  n_tags = length(tags)
+  @assert n_tags == n_fields
+  PETScFieldSplit(MF, tags, show_idx)
+end
+
+
+
+
+function field_dof_split(U)
+  X = MultiFieldFESpace(U)
+  field_dof_split(X, U)
+end
+
+# It allows to obtain the dof of each field over each processor - dont know if the best/most elegant way
+function field_dof_split(X, U)
+  @assert length(X) == length(U)
+  xrows = X.gids.partition
+  urows = Any[]
+  for Ui in X
+    urowsi = Ui.gids.partition
+    push!(urows, urowsi)
+  end
+
+  offset = 0
+  ulrows = Any[]
+  for urowsi in urows
+    ulrowsi = map_parts(urowsi, xrows) do urowsij, xrows
+      num_dof_Ui = length(urowsij.oid_to_lid)
+      uloc_idx = collect(1+offset:1:num_dof_Ui+offset)
+      uri = xrows.lid_to_gid[xrows.oid_to_lid][uloc_idx]
+      offset = num_dof_Ui + offset
+      return uri
+    end
+    #It works only in MPI, non-MPI do not have .part fields
+    # .-1 to match the numbering in PETSc that starts from 0
+    ur = ulrowsi.part .- 1
+    push!(ulrows, ur)
+  end
+  return ulrows
+end
+
+
+
+function set_fieldsplit(sol, SplitField::PETScFieldSplit)
+
+  #Get PC-string name
+  pc = Ref{GridapPETSc.PETSC.PC}()
+  pctype = Ref{Ptr{Cstring}}()
+
+  #Get the PC, if the problem is linear ksp or non linear snes
+  #Maybe not the most elegant solution
+  if typeof(sol) <: Base.RefValue{GridapPETSc.PETSC.KSP}
+    ksp = sol
+  elseif typeof(sol) <: Base.RefValue{GridapPETSc.PETSC.SNES}
+    snes = sol
+    ksp = Ref{GridapPETSc.PETSC.KSP}()
+    @check_error_code GridapPETSc.PETSC.SNESGetKSP(snes[], ksp)
+  end
+
+  @check_error_code GridapPETSc.PETSC.KSPGetPC(ksp[], pc)
+  @check_error_code GridapPETSc.PETSC.PCGetType(pc[], pctype)
+
+  pc_ptr_conv = reinterpret(Ptr{UInt8}, pctype[])
+  GC.@preserve pc_name = unsafe_string(pc_ptr_conv)
+
+  @assert pc_name == "fieldsplit" #Check that the preconditioner requires splitting fields
+  n_tags = length(SplitField.tags)
+  U_Parray = field_dof_split(SplitField.MF)
+  for i = 1:1:n_tags
+    ISU = PETScIS(U_Parray[i])
+    @check_error_code GridapPETSc.PETSC.PCFieldSplitSetIS(pc[], SplitField.tags[i], ISU.is[])
+    if SplitField.show_idx
+      @check_error_code GridapPETSc.PETSC.ISView(ISU.is[], GridapPETSc.PETSC.@PETSC_VIEWER_STDOUT_SELF)
+    end
+  end
+
+
+end
\ No newline at end of file
diff --git a/src/PETScIndexes.jl b/src/PETScIndexes.jl
new file mode 100644
index 0000000..61221b4
--- /dev/null
+++ b/src/PETScIndexes.jl
@@ -0,0 +1,81 @@
+# Index
+
+mutable struct PETScIS <: AbstractVector{PetscInt}
+  is::Base.RefValue{IS}
+  initialized::Bool
+  size::Tuple{Int}
+  comm::MPI.Comm
+  PETScIS(comm::MPI.Comm) = new(Ref{IS}(), false, (-1,), comm)
+end
+
+
+function Init(a::PETScIS)
+  n = Ref{PetscInt}()
+  @check_error_code PETSC.ISGetSize(a.is[], n)
+  a.size = (Int(n[]),)
+  @assert Threads.threadid() == 1
+  _NREFS[] += 1
+  a.initialized = true
+  finalizer(Finalize, a)
+end
+
+function Finalize(a::PETScIS)
+  if a.initialized && GridapPETSc.Initialized()
+    if a.comm == MPI.COMM_SELF
+      @check_error_code PETSC.ISDestroy(a.is)
+    else
+      @check_error_code PETSC.PetscObjectRegisterDestroy(a.is[])
+    end
+    a.initialized = false
+    @assert Threads.threadid() == 1
+    _NREFS[] -= 1
+  end
+  nothing
+end
+
+function Base.size(v::PETScIS)
+  @assert v.initialized
+  v.size
+end
+
+Base.@propagate_inbounds function Base.getindex(v::PETScIS, i1::Integer)
+  @boundscheck checkbounds(v, i1)
+  n = one(PetscInt)
+  i0 = Ref(i1 - n)
+  pi0 = reinterpret(Ptr{PetscInt}, pointer_from_objref(i0))
+  @check_error_code PETSC.ISGetIndices(v.is[], pi0)
+  return pi0
+end
+
+function PETScIS(idx::Vector{PetscInt}, bs=1)
+  comm = MPI.COMM_SELF
+  is = PETScIS(comm)
+  n = length(idx)
+  @check_error_code GridapPETSc.PETSC.ISCreateGeneral(comm, n, idx, GridapPETSc.PETSC.PETSC_COPY_VALUES, is.is)
+  Init(is)
+end
+
+function PETScIS(idx::AbstractVector, bs=1)
+  idx = PetscInt.(idx)
+  PETScIS(idx)
+end
+
+# #Block Implementation
+#   function PETScIS(array::Vector{PetscInt},n, bs=1)
+#     comm = MPI.COMM_SELF
+#     is = PETScIS(comm)
+#     n = PetscInt(n)
+#     bs = PetscInt(bs)
+#     @check_error_code GridapPETSc.PETSC.ISCreateBlock(comm, n, bs, array, GridapPETSc.PETSC.PETSC_COPY_VALUES, is.is)
+#     Init(is)
+#   end
+
+# Constructors
+
+function PETScIS(n::Integer)
+  println("PETScIS")
+  v = Ref{Ptr{PetscInt}}()
+  println("Construct\n")
+  @check_error_code PETSC.PetscMalloc1(n, v)
+  Init(v)
+end
\ No newline at end of file
diff --git a/src/PETScLinearSolvers.jl b/src/PETScLinearSolvers.jl
index ff993a2..c87d5c3 100644
--- a/src/PETScLinearSolvers.jl
+++ b/src/PETScLinearSolvers.jl
@@ -1,19 +1,28 @@
-
 struct PETScLinearSolver{F} <: LinearSolver
   setup::F
+  fieldsplit::Union{PETScFieldSplit, Nothing}
 end
 
 ksp_from_options(ksp) = @check_error_code PETSC.KSPSetFromOptions(ksp[])
 
+function PETScLinearSolver(SplitField::PETScFieldSplit)
+  PETScLinearSolver(ksp_from_options,SplitField)
+end
+
+function PETScLinearSolver(ksp_options)
+  PETScLinearSolver(ksp_options,nothing)
+end
+
 function PETScLinearSolver()
-  PETScLinearSolver(ksp_from_options)
+  PETScLinearSolver(ksp_from_options, nothing)
 end
 
+
 struct PETScLinearSolverSS{F} <: SymbolicSetup
   solver::PETScLinearSolver{F}
 end
 
-function Algebra.symbolic_setup(solver::PETScLinearSolver,mat::AbstractMatrix)
+function Algebra.symbolic_setup(solver::PETScLinearSolver, mat::AbstractMatrix)
   PETScLinearSolverSS(solver)
 end
 
@@ -22,9 +31,9 @@ mutable struct PETScLinearSolverNS{T} <: NumericalSetup
   B::PETScMatrix
   ksp::Ref{KSP}
   initialized::Bool
-  function PETScLinearSolverNS(A,B::PETScMatrix)
-    T=typeof(A)
-    new{T}(A,B,Ref{KSP}(),false)
+  function PETScLinearSolverNS(A, B::PETScMatrix)
+    T = typeof(A)
+    new{T}(A, B, Ref{KSP}(), false)
   end
 end
 
@@ -32,7 +41,7 @@ function Init(a::PETScLinearSolverNS)
   @assert Threads.threadid() == 1
   _NREFS[] += 1
   a.initialized = true
-  finalizer(Finalize,a)
+  finalizer(Finalize, a)
 end
 
 function Finalize(ns::PETScLinearSolverNS)
@@ -49,53 +58,58 @@ function Finalize(ns::PETScLinearSolverNS)
   nothing
 end
 
-function Algebra.numerical_setup(ss::PETScLinearSolverSS,A::AbstractMatrix)
-  B = convert(PETScMatrix,A)
-  ns = PETScLinearSolverNS(A,B)
-  @check_error_code PETSC.KSPCreate(B.comm,ns.ksp)
-  @check_error_code PETSC.KSPSetOperators(ns.ksp[],ns.B.mat[],ns.B.mat[])
+function Algebra.numerical_setup(ss::PETScLinearSolverSS, A::AbstractMatrix)
+  B = convert(PETScMatrix, A)
+  ns = PETScLinearSolverNS(A, B)
+  @check_error_code PETSC.KSPCreate(B.comm, ns.ksp)
+  @check_error_code PETSC.KSPSetOperators(ns.ksp[], ns.B.mat[], ns.B.mat[])
   ss.solver.setup(ns.ksp)
+
+  if typeof(ss.solver.fieldsplit) == PETScFieldSplit
+    set_fieldsplit(ns.ksp, ss.solver.fieldsplit)
+  end
+
   @check_error_code PETSC.KSPSetUp(ns.ksp[])
   Init(ns)
 end
 
-function Algebra.solve!(x::PETScVector,ns::PETScLinearSolverNS,b::AbstractVector)
+function Algebra.solve!(x::PETScVector, ns::PETScLinearSolverNS, b::AbstractVector)
   if (x.comm != MPI.COMM_SELF)
     gridap_petsc_gc() # Do garbage collection of PETSc objects
   end
 
-  B = convert(PETScVector,b)
-  @check_error_code PETSC.KSPSolve(ns.ksp[],B.vec[],x.vec[])
+  B = convert(PETScVector, b)
+  @check_error_code PETSC.KSPSolve(ns.ksp[], B.vec[], x.vec[])
   x
 end
 
-function Algebra.solve!(x::Vector{PetscScalar},ns::PETScLinearSolverNS,b::AbstractVector)
-  X = convert(PETScVector,x)
-  solve!(X,ns,b)
+function Algebra.solve!(x::Vector{PetscScalar}, ns::PETScLinearSolverNS, b::AbstractVector)
+  X = convert(PETScVector, x)
+  solve!(X, ns, b)
   x
 end
 
-function Algebra.solve!(x::AbstractVector,ns::PETScLinearSolverNS,b::AbstractVector)
-  X = convert(Vector{PetscScalar},x)
-  solve!(X,ns,b)
+function Algebra.solve!(x::AbstractVector, ns::PETScLinearSolverNS, b::AbstractVector)
+  X = convert(Vector{PetscScalar}, x)
+  solve!(X, ns, b)
   x .= X
   x
 end
 
-function Algebra.solve!(x::PVector,ns::PETScLinearSolverNS,b::PVector)
-  X = similar(b,(axes(ns.A)[2],))
-  B = similar(b,(axes(ns.A)[2],))
-  copy!(X,x)
-  copy!(B,b)
-  Y = convert(PETScVector,X)
-  solve!(Y,ns,B)
-  copy!(x,Y)
+function Algebra.solve!(x::PVector, ns::PETScLinearSolverNS, b::PVector)
+  X = similar(b, (axes(ns.A)[2],))
+  B = similar(b, (axes(ns.A)[2],))
+  copy!(X, x)
+  copy!(B, b)
+  Y = convert(PETScVector, X)
+  solve!(Y, ns, B)
+  copy!(x, Y)
 end
 
-function Algebra.numerical_setup!(ns::PETScLinearSolverNS,A::AbstractMatrix)
+function Algebra.numerical_setup!(ns::PETScLinearSolverNS, A::AbstractMatrix)
   ns.A = A
-  ns.B = convert(PETScMatrix,A)
-  @check_error_code PETSC.KSPSetOperators(ns.ksp[],ns.B.mat[],ns.B.mat[])
+  ns.B = convert(PETScMatrix, A)
+  @check_error_code PETSC.KSPSetOperators(ns.ksp[], ns.B.mat[], ns.B.mat[])
   @check_error_code PETSC.KSPSetUp(ns.ksp[])
   ns
 end
diff --git a/src/PETScNonlinearSolvers.jl b/src/PETScNonlinearSolvers.jl
index ecd630b..47a1ea1 100644
--- a/src/PETScNonlinearSolvers.jl
+++ b/src/PETScNonlinearSolvers.jl
@@ -1,6 +1,7 @@
 
 mutable struct PETScNonlinearSolver{F} <: NonlinearSolver
   setup::F
+  fieldsplit::Union{PETScFieldSplit, Nothing}
 end
 
 mutable struct PETScNonlinearSolverCache{A,B,C,D,E}
@@ -108,8 +109,18 @@ end
 
 snes_from_options(snes) = @check_error_code PETSC.SNESSetFromOptions(snes[])
 
+
+function PETScNonlinearSolver(SplitField::PETScFieldSplit)
+  PETScNonlinearSolver(snes_from_options,SplitField)
+end
+
+function PETScNonlinearSolver(snes_options)
+  PETScNonlinearSolver(snes_options,nothing)
+end
+
+
 function PETScNonlinearSolver()
-  PETScNonlinearSolver(snes_from_options)
+  PETScNonlinearSolver(snes_from_options,nothing)
 end
 
 function _set_petsc_residual_function!(nls::PETScNonlinearSolver, cache)
@@ -182,6 +193,10 @@ function Algebra.solve!(x::AbstractVector,nls::PETScNonlinearSolver,op::Nonlinea
 
   nls.setup(cache.snes)
 
+  if typeof(nls.fieldsplit) == PETScFieldSplit
+    set_fieldsplit(cache.snes, nls.fieldsplit)
+  end
+  
   @check_error_code PETSC.SNESSolve(cache.snes[],C_NULL,cache.x_petsc.vec[])
   copy!(x,cache.x_petsc)
   cache
diff --git a/test/IndexesTests.jl b/test/IndexesTests.jl
new file mode 100644
index 0000000..f3e29a6
--- /dev/null
+++ b/test/IndexesTests.jl
@@ -0,0 +1,33 @@
+using GridapPETSc
+using Test
+using SparseArrays
+using SparseMatricesCSR
+using GridapPETSc: PetscScalar, PetscInt
+using LinearAlgebra
+using MPI
+using PartitionedArrays
+
+
+function main(parts)
+options = "-info"
+  GridapPETSc.with(args=split(options)) do
+    backend = get_backend(parts)
+    if backend == MPIBackend()
+      comm = MPI.COMM_WORLD
+      procid = PartitionedArrays.get_part_id(comm)
+      nprocs = PartitionedArrays.num_parts(comm)
+    elseif backend == SequentialBackend()
+      procid = 1
+      nprocs = 1
+    end
+
+    println(procid)
+    array = ones(5)
+    array = array .+ procid
+    is =  PETScIS(array)
+    @check_error_code GridapPETSc.PETSC.ISView(is.is[], GridapPETSc.PETSC.@PETSC_VIEWER_STDOUT_WORLD)
+    @test is.size[1] == length(array)*nprocs
+  end
+end
+
+#mpiexecjl --project=. -n 4 julia test/IndexesTests.jl
\ No newline at end of file
diff --git a/test/mpi/IndexesTests.jl b/test/mpi/IndexesTests.jl
new file mode 100644
index 0000000..ebb49d1
--- /dev/null
+++ b/test/mpi/IndexesTests.jl
@@ -0,0 +1,3 @@
+include("../IndexesTests.jl")
+nparts = (2,2)
+prun(main,mpi,nparts)
diff --git a/test/sequential/PETScIndexesTests.jl b/test/sequential/PETScIndexesTests.jl
new file mode 100644
index 0000000..2f9af77
--- /dev/null
+++ b/test/sequential/PETScIndexesTests.jl
@@ -0,0 +1,25 @@
+module PETScIndexesTests
+
+using GridapPETSc
+using Test
+using SparseArrays
+using SparseMatricesCSR
+using GridapPETSc: PetscScalar, PetscInt
+using LinearAlgebra
+
+
+
+
+options = "-info"
+GridapPETSc.with(args=split(options)) do
+
+array = collect(1:1:100)
+is = PETScIS(array)
+#@check_error_code GridapPETSc.PETSC.ISView(is.is[], GridapPETSc.PETSC.@PETSC_VIEWER_STDOUT_WORLD)
+@check_error_code GridapPETSc.PETSC.ISView(is.is[], GridapPETSc.PETSC.@PETSC_VIEWER_STDOUT_SELF)
+@test is.size[1] == length(array)
+@test is.initialized == true
+end
+
+
+end # module
diff --git a/test/sequential/StokesTest.jl b/test/sequential/StokesTest.jl
new file mode 100644
index 0000000..c0db00b
--- /dev/null
+++ b/test/sequential/StokesTest.jl
@@ -0,0 +1,114 @@
+module StokesTest
+
+using Gridap
+using Test
+using GridapDistributed
+
+using PartitionedArrays
+using SparseArrays
+using GridapPETSc
+
+#You can provide the following options string (which will be automatically used when calling PetscLinearSolver())
+# or the PescLinearSolver(mykspsetup)
+#They are equivalent (with the exeption of -ksp_monitor)
+#You can use PescLinearSolver(mykspsetup) and at the same time "-ksp_monitor", they are merged
+
+function mykspsetup(ksp)
+    pc = Ref{GridapPETSc.PETSC.PC}()
+    umfpack = Ref{GridapPETSc.PETSC.Mat}()
+    @check_error_code GridapPETSc.PETSC.KSPSetType(ksp[], GridapPETSc.PETSC.KSPPREONLY)
+    @check_error_code GridapPETSc.PETSC.KSPGetPC(ksp[], pc)
+    @check_error_code GridapPETSc.PETSC.PCSetType(pc[], GridapPETSc.PETSC.PCLU)
+    @check_error_code GridapPETSc.PETSC.PCFactorSetMatSolverType(pc[], GridapPETSc.PETSC.MATSOLVERUMFPACK)
+    @check_error_code GridapPETSc.PETSC.PCFactorSetUpMatSolverType(pc[])
+    @check_error_code GridapPETSc.PETSC.PCFactorGetMatrix(pc[], umfpack)
+    @check_error_code GridapPETSc.PETSC.KSPSetFromOptions(ksp[])
+end
+
+
+function main(solver_options)
+    if solver_options == :petsc_linecommand
+        options = "-pc_type lu -ksp_type preonly -ksp_max_it 10 -ksp_monitor -pc_factor_mat_solver_type umfpack"
+    elseif solver_options == :petsc_mykspsetup
+        options = "-ksp_monitor"
+    elseif solver_options == :julia
+        options = " "
+    else
+        error("Solver $(solver_options) not valid. Use instead:\n\
+         :petsc_linecommand\n\
+         :petsc_mykspsetup\n\
+         :julia\n")
+    end
+    tt = 0.0
+    GridapPETSc.with(args=split(options)) do
+        n = 50
+        domain = (0, 1, 0, 1)
+        partition = (n, n)
+        model = CartesianDiscreteModel(domain, partition)
+
+        labels = get_face_labeling(model)
+        add_tag_from_tags!(labels, "diri1", [6,])
+        add_tag_from_tags!(labels, "diri0", [1, 2, 3, 4, 5, 7, 8])
+        add_tag_from_tags!(labels, "dirip", [1])
+
+        order = 2
+        reffeᵤ = ReferenceFE(lagrangian, VectorValue{2,Float64}, order)
+        reffeₚ = ReferenceFE(lagrangian, Float64, order - 1; space=:P)
+
+        V = TestFESpace(model, reffeᵤ, labels=labels, dirichlet_tags=["diri0", "diri1"], conformity=:H1)
+        Q = TestFESpace(model, reffeₚ, conformity=:L2, dirichlet_tags="dirip")
+        Y = MultiFieldFESpace([V, Q])
+
+        u0 = VectorValue(0, 0)
+        u1 = VectorValue(1, 0)
+        U = TrialFESpace(V, [u0, u1])
+        P = TrialFESpace(Q, 0.0)
+        X = MultiFieldFESpace([U, P])
+
+        degree = order
+        Ω = Triangulation(model)
+        dΩ = Measure(Ω, degree)
+
+        f = VectorValue(0.0, 0.0)
+        h = 1 / n
+        τ = (h .^ 2) ./ 4
+
+        a((u, p), (v, q)) = ∫(∇(v) ⊙ ∇(u) - (∇ ⋅ v) * p + q * (∇ ⋅ u))dΩ + ∫((τ ⋅ ∇(q))' ⋅ (∇(p)))dΩ
+        #The last term is added in order to have elements on the main diagonal, if not PETSc does not work properly
+
+        l((v, q)) = ∫(v ⋅ f)dΩ
+        res((u, p), (v, q)) = a((u, p), (v, q)) - l((v, q))
+        op = AffineFEOperator(a, l, X, Y)
+
+
+        if solver_options == :petsc_linecommand
+            #Solve using Petsc solver, with the options given in line command style
+            solver = PETScLinearSolver()
+            uh, ph = solve(solver, op)
+            tt = @elapsed uh, ph = solve(solver, op)
+        elseif solver_options == :petsc_mykspsetup
+            #Solve using Petsc solver, with the options given in  mykspsetup
+            solver = PETScLinearSolver(mykspsetup)
+            uh, ph = solve(solver, op)
+            tt = @elapsed uh, ph = solve(solver, op)
+        elseif solver_options == :julia
+            #Solve using default julia solver, which is a LU decomposition
+            uh, ph = solve(op)
+            tt = @elapsed uh, ph = solve(op)
+        end
+
+    end
+    return tt
+
+end
+#The PETSc and Julia in this case both use a LU decomposition, the time should be really close
+
+j = main(:julia)
+pl = main(:petsc_linecommand)
+pm = main(:petsc_mykspsetup)
+
+@test isapprox(j,pl; rtol = 0.1)
+@test isapprox(j,pm; rtol = 0.1)
+@test isapprox(pm,pm; rtol = 0.1)
+
+end #end module
diff --git a/test/sequential/runtests.jl b/test/sequential/runtests.jl
index a42e404..5d53ce0 100644
--- a/test/sequential/runtests.jl
+++ b/test/sequential/runtests.jl
@@ -7,6 +7,8 @@ using MPI
 
 @time @testset "PETScArrays" begin include("PETScArraysTests.jl") end
 
+@time @testset "PETScIndexes" begin include("PETScIndexesTests.jl") end
+
 @time @testset "PartitionedArrays (sequential)" begin include("PartitionedArraysTests.jl") end
 
 @time @testset "PETScLinearSolvers" begin include("PETScLinearSolversTests.jl") end
@@ -19,6 +21,8 @@ using MPI
 
 @time @testset "ElasticityDriver" begin include("ElasticityDriver.jl") end
 
+@time @testset "StokesTest" begin include("StokesTest.jl") end
+
 @time @testset "DarcyDriver" begin include("DarcyDriver.jl") end
 
 @time @testset "PLaplacianDriver" begin include("PLaplacianDriver.jl") end