diff --git a/.github/workflows/test-PR-struphy.yml b/.github/workflows/test-PR-struphy.yml
index 221135a89..97adbaeb9 100644
--- a/.github/workflows/test-PR-struphy.yml
+++ b/.github/workflows/test-PR-struphy.yml
@@ -63,13 +63,6 @@ jobs:
         run: |
           ls .testmon* || echo "No .testmondata"
 
-      # # temporary - remove after feectools merge
-      # - name: Install Struphy in container
-      #   uses: ./.github/actions/install/struphy-in-container
-      #   with:
-      #     compile_language: ${{ matrix.compile_language }}
-      #     struphy_branch: 146-install-feectools
-
       - name: Install Psydac in container
         uses: ./.github/actions/install/psydac-in-container
         with:
@@ -84,17 +77,35 @@ jobs:
       - name: Run struphy unit tests
         env:
           TESTMON_DATAFILE: ${{ github.workspace }}/.testmondata-struphy-${{ matrix.compile_language }}
+          TEST_FILES: "bsplines/tests/ \
+                      console/tests/ \
+                      feec/tests/ \
+                      fields_background/tests/ \
+                      geometry/tests/ \
+                      initial/tests/ \
+                      kinetic_background/tests/ \
+                      linear_algebra/tests/ \
+                      ode/tests/ \
+                      pic/tests/ \
+                      polar/tests/ \
+                      post_processing/tests/ \
+                      propagators/tests/"
         run: |
           source /struphy_${{ matrix.compile_language }}_/env_${{ matrix.compile_language }}_/bin/activate
-          struphy test Maxwell
-          struphy test unit
+          PYTHON=$(which python)
+          STRUPHY_PATH=$($PYTHON -c "import struphy; print(struphy.__path__[0])")
+          echo "Struphy is installed at: $STRUPHY_PATH"
+          echo "STRUPHY_PATH=${STRUPHY_PATH}" >> $GITHUB_ENV
+          mpirun -n 1 pytest -m single --testmon-forceselect -s --with-mpi --model-name Maxwell $STRUPHY_PATH
+          cd $STRUPHY_PATH
+          pytest --testmon $TEST_FILES
       
       - name: Run struphy model tests
         env:
           TESTMON_DATAFILE: ${{ github.workspace }}/.testmondata-struphy-${{ matrix.compile_language }}
         run: |
           source /struphy_${{ matrix.compile_language }}_/env_${{ matrix.compile_language }}_/bin/activate
-          struphy test models 
+          pytest -m models --testmon-forceselect  $STRUPHY_PATH/models/tests/default_params/  
 
       - name: Install mpi4py
         run: |
@@ -106,5 +117,5 @@ jobs:
           TESTMON_DATAFILE: ${{ github.workspace }}/.testmondata-struphy-${{ matrix.compile_language }}-mpi
         run: |
           source /struphy_${{ matrix.compile_language }}_/env_${{ matrix.compile_language }}_/bin/activate
-          struphy test Maxwell
-          struphy test models --mpi 2
\ No newline at end of file
+          mpirun -n 1 pytest -m single --testmon-forceselect -s --with-mpi --model-name Maxwell $STRUPHY_PATH
+          mpirun -oversubscribe -n 2 pytest -m models --testmon-forceselect --with-mpi $STRUPHY_PATH/models/tests/default_params/
\ No newline at end of file
diff --git a/.github/workflows/test-struphy.yml b/.github/workflows/test-struphy.yml
deleted file mode 100644
index 10608aa09..000000000
--- a/.github/workflows/test-struphy.yml
+++ /dev/null
@@ -1,112 +0,0 @@
-name: Test struphy
-
-on:
-  push:
-    branches: [devel-tiny]
-  workflow_dispatch:
-
-concurrency:
-  group: ${{ github.workflow }}-${{ github.ref }}
-  cancel-in-progress: true
-
-jobs:
-  test_struphy:
-    runs-on: ${{ matrix.os }}
-    strategy:
-      fail-fast: false
-      matrix:
-        os: [ ubuntu-latest ]
-        python-version: [ '3.12' ]
-        compile_language: ['c', 'fortran']
-
-    name: Test Struphy on ${{ matrix.os }} / Python ${{ matrix.python-version }} / ${{ matrix.compile_language }}
-
-    steps:
-      - name: Checkout repo
-        uses: actions/checkout@v4
-
-      - name: Set up environment variables
-        run: |
-          echo "PSYDAC_DIR=$GITHUB_WORKSPACE" >> $GITHUB_ENV
-          echo "STRUPHY_DIR=$GITHUB_WORKSPACE/struphy" >> $GITHUB_ENV
-
-      - name: Check .testmondata 1
-        run: |
-          ls .testmon* || echo "No .testmondata"
-
-      - name: Setup cache for testmon
-        uses: actions/cache@v4
-        with:
-          path: |
-            .testmondata-struphy-${{ matrix.compile_language }}
-            .testmondata-struphy-${{ matrix.compile_language }}-mpi
-          key: testmon-struphy-${{ matrix.compile_language }}-${{ github.ref_name }}-${{ github.event.number }}-${{ github.run_number }}
-          restore-keys: |
-            testmon-struphy-${{ matrix.compile_language }}-${{ github.ref_name }}-${{ github.event.number }}-
-            testmon-struphy-${{ matrix.compile_language }}-${{ github.ref_name }}
-            testmon-struphy-${{ matrix.compile_language }}-devel
-            testmon-struphy-${{ matrix.compile_language }}-
-            testmon-struphy-
-
-      - name: Check .testmondata 2
-        run: |
-          ls .testmon* || echo "No .testmondata"
-
-      - name: Install prerequisites
-        uses: ./.github/actions/install/psydac-req
-
-      - name: Install project without mpi
-        run: |
-          python -m pip install ".[test]" --no-cache-dir
-          python -m pip freeze
-
-      - name: Clone struphy from Github
-        run: |
-          git clone https://github.com/struphy-hub/struphy.git $STRUPHY_DIR
-      
-      - name: Install struphy without mpi #TODO: Set branch to devel
-        working-directory: ${{ env.STRUPHY_DIR }}
-        run: |
-          pip show feectools
-          pip uninstall feectools -y
-          python -m pip install ".[phys]" --no-cache-dir
-          echo "feectools location after installing struphy"
-          pip show feectools
-          
-      - name: Re-install feectools from the current branch
-        run: |
-          psydac-accelerate --cleanup --yes
-          pip show feectools
-          pip uninstall feectools -y
-          python -m pip install ${{ env.PSYDAC_DIR }} --force-reinstall --no-cache-dir
-          echo "feectools location after installing feectools from the current branch (should be the same as in the previous step)"
-          pip show feectools
-          
-      - name: Compile struphy
-        run: |
-          pip show feectools
-          struphy -h
-          struphy compile -y --language ${{ matrix.compile_language }}
-
-      - name: Run struphy unit tests
-        env:
-          TESTMON_DATAFILE: ${{ github.workspace }}/.testmondata-struphy-${{ matrix.compile_language }}
-        run: |
-            struphy test unit
-      
-      - name: Run struphy model tests
-        env:
-          TESTMON_DATAFILE: ${{ github.workspace }}/.testmondata-struphy-${{ matrix.compile_language }}
-        run: |
-            struphy test models 
-
-      - name: Install mpi4py
-        run: |
-            pip install -U mpi4py
-
-      - name: Run struphy model tests with mpi
-        env:
-          TESTMON_DATAFILE: ${{ github.workspace }}/.testmondata-struphy-${{ matrix.compile_language }}-mpi
-        run: |
-            struphy test Maxwell
-            struphy test models --mpi 2
diff --git a/feectools/linalg/basic.py b/feectools/linalg/basic.py
index 0497e7092..872a7027c 100644
--- a/feectools/linalg/basic.py
+++ b/feectools/linalg/basic.py
@@ -756,7 +756,7 @@ def set_scalar(self, c):
         self._scalar = c
 
     def toarray(self):
-        return self._scalar * self._operator.toarray() 
+        return self._scalar * self._operator.toarray
 
     def tosparse(self):
         return self._scalar * self._operator.tosparse().tocsr()
diff --git a/feectools/linalg/solvers.py b/feectools/linalg/solvers.py
index 31f084d6f..dbbac8e58 100644
--- a/feectools/linalg/solvers.py
+++ b/feectools/linalg/solvers.py
@@ -4,12 +4,14 @@
 
 """
 import numpy as np
-from math import sqrt
+from math import sqrt, inf
 
 from feectools.utilities.utils  import is_real
 from feectools.linalg.utilities import _sym_ortho
 from feectools.linalg.basic     import (Vector, LinearOperator,
         InverseLinearOperator, IdentityOperator, ScaledLinearOperator)
+from feectools.linalg.block import BlockVector, BlockVectorSpace
+
 
 __all__ = (
     'inverse',
@@ -20,7 +22,7 @@
     'PBiConjugateGradientStabilized',
     'MinimumResidual',
     'LSMR',
-    'GMRES'
+    'GMRES',
 )
 
 #===============================================================================
@@ -64,6 +66,8 @@ def inverse(A, solver, **kwargs):
         'minres'   : MinimumResidual,
         'lsmr'     : LSMR,
         'gmres'    : GMRES,
+        'uzawa'    : UzawaSolver,
+        'schur'    : SchurComplementSolver 
     }
 
     # Check solver input
@@ -1912,3 +1916,219 @@ def apply_givens_rotation(self, k, sn, cn):
     def dot(self, b, out=None):
         return self.solve(b, out=out)
 
+
+class UzawaSolver(InverseLinearOperator):
+    def __init__(self, A, *, A11, A22, B1, B2,
+                x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle=False,
+                inner_tol=1e-7, inner_maxiter=1000):
+
+
+        self._options = {
+            "x0": x0, "tol": tol, "maxiter": maxiter,
+            "verbose": verbose, "recycle": recycle,
+        }
+        super().__init__(A, **self._options)
+
+        self._inner_tol = inner_tol if inner_tol is not None else tol
+        self._inner_maxiter = inner_maxiter
+
+        self._A11 = A11
+        self._A22 = A22
+        self._B1  = B1
+        self._B2  = B2
+
+        self._A11inv = self._inner_solve(A11, tol=self._inner_tol, maxiter=self._inner_maxiter)
+        self._A22inv = self._inner_solve(A22, tol=self._inner_tol, maxiter=self._inner_maxiter)
+
+        self._info = None
+
+    @staticmethod
+    def _inner_solve(A, tol=1e-10, maxiter=200):
+        """Iterative GMRES inverse — works for non-symmetric positive definite operators."""
+        from feectools.linalg.solvers import inverse
+        return inverse(A, "gmres", tol=tol, maxiter=maxiter, recycle=True)
+
+    def update_A11(self, A11):
+        """Recompute A11 inverse when dt changes."""
+        self._A11 = A11
+        self._A11inv = self._inner_solve(A11, tol=self._inner_tol, maxiter=self._inner_maxiter)
+
+    def solve(self, b, out=None):
+
+        A11    = self._A11
+        A22    = self._A22
+        B1     = self._B1
+        B2     = self._B2
+        A11inv = self._A11inv
+        A22inv = self._A22inv
+
+        tol     = self._options["tol"]
+        maxiter = self._options["maxiter"]
+        verbose = self._options["verbose"]
+        recycle = self._options["recycle"]
+
+        F    = b[0]
+        g    = b[1]
+        f_u  = F[0]
+        f_ue = F[1]
+
+        # initial guess
+        x0 = self._options["x0"]
+        if x0 is not None:
+            u  = x0[0][0].copy()
+            ue = x0[0][1].copy()
+            p  = x0[1].copy()
+        else:
+            u  = A11.domain.zeros()
+            ue = A22.domain.zeros()
+            p  = B1.codomain.zeros()
+
+
+        if verbose:
+            print("Uzawa solver:")
+            print("+---------+---------------------+")
+            print("+ Iter. # | L2-norm of residual |")
+            print("+---------+---------------------+")
+            template = "| {:7d} | {:19.2e} |"
+
+        for iteration in range(1, maxiter + 1):
+
+            # solve A11 * u = f_u - B1^T * p
+            rhs_u = f_u - B1.T.dot(p)
+            u = A11inv.dot(rhs_u)
+
+            # solve A22 * ue = f_ue - B2^T * p
+            rhs_ue = f_ue - B2.T.dot(p)
+            ue = A22inv.dot(rhs_ue)
+
+            # constraint residual: R = B1*u + B2*ue - g
+            R = B1.dot(u) + B2.dot(ue) - g
+            residual_norm = sqrt(R.inner(R).real)
+
+            if verbose:
+                print(template.format(iteration, residual_norm))
+
+            if residual_norm < tol:
+                break
+
+            # pressure update: steepest descent step size
+            S_R   = B1.dot(A11inv.dot(B1.T.dot(R))) + B2.dot(A22inv.dot(B2.T.dot(R)))
+            alpha = R.inner(R).real / R.inner(S_R).real
+            p    += alpha * R
+
+        if verbose:
+            print("+---------+---------------------+")
+
+        self._info = {
+            'niter'   : iteration,
+            'success' : residual_norm < tol,
+            'res_norm': residual_norm,
+        }
+
+        if recycle:
+            block_u = BlockVector(self.domain.spaces[0], blocks=[u, ue])
+            self._options["x0"] = BlockVector(self.domain, blocks=[block_u, p])
+
+        if out is not None:
+            out[0][0] = u
+            out[0][1] = ue
+            out[1]    = p
+            return out
+
+        block_u = BlockVector(BlockVectorSpace(A11.domain, A22.domain), blocks=[u, ue])
+        return BlockVector(self.domain, blocks=[block_u, p])
+
+    def dot(self, b, out=None):
+        return self.solve(b, out=out)
+
+
+class SchurComplementSolver(InverseLinearOperator):
+    def __init__(self, A, *, A11, A22, B1, B2,
+                 x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle=False,
+                 inner_tol=1e-7, inner_maxiter=1000):
+
+        self._options = {
+            "x0": x0, "tol": tol, "maxiter": maxiter,
+            "verbose": verbose, "recycle": recycle,
+        }
+        super().__init__(A, **self._options)
+
+        self._inner_tol     = inner_tol if inner_tol is not None else tol
+        self._inner_maxiter = inner_maxiter
+
+        self._A11 = A11
+        self._A22 = A22
+        self._B1  = B1
+        self._B2  = B2
+
+        self._A11inv = self._make_inner_solver(A11, self._inner_tol, self._inner_maxiter)
+        self._A22inv = self._make_inner_solver(A22, self._inner_tol, self._inner_maxiter)
+
+        self._Sinv = self._build_schur_solver()
+
+        self._info = None
+
+    def _make_inner_solver(self, A, tol, maxiter):
+        from feectools.linalg.solvers import inverse
+        return inverse(A, "gmres", tol=tol, maxiter=maxiter, recycle=True)
+
+    def _build_schur_operator(self):
+        return (self._B1 @ self._A11inv @ self._B1.T
+              + self._B2 @ self._A22inv @ self._B2.T)
+
+    def _build_schur_solver(self):
+        from feectools.linalg.solvers import inverse
+        S = self._build_schur_operator()
+        return inverse(S, "gmres",
+                       tol=self._options["tol"],
+                       maxiter=self._options["maxiter"],
+                       verbose=self._options["verbose"],
+                       recycle=self._options["recycle"])
+
+    def update_A11(self, A11):
+        """Recompute A11 inverse and update Schur operator when dt changes."""
+        self._A11    = A11
+        self._A11inv = self._make_inner_solver(A11, self._inner_tol, self._inner_maxiter)
+        # update operator in-place so Sinv keeps its recycled state
+        self._Sinv.linop = self._build_schur_operator()
+
+    def solve(self, b, out=None):
+
+        A11inv = self._A11inv
+        A22inv = self._A22inv
+        B1     = self._B1
+        B2     = self._B2
+        Sinv   = self._Sinv
+
+        F    = b[0]
+        g    = b[1]
+        f_u  = F[0]
+        f_ue = F[1]
+
+        # Schur complement RHS: B1*u_f + B2*ue_f - g
+        rhs_p = B1.dot(A11inv.dot(f_u)) + B2.dot(A22inv.dot(f_ue)) - g
+
+        # solve S * p = rhs_p
+        p = Sinv.dot(rhs_p)
+
+        # recover velocities
+        u  = A11inv.dot(f_u  - B1.T.dot(p))
+        ue = A22inv.dot(f_ue - B2.T.dot(p))
+
+        self._info = Sinv.get_info()
+
+        if self._options["recycle"]:
+            block_u = BlockVector(self.domain.spaces[0], blocks=[u, ue])
+            self._options["x0"] = BlockVector(self.domain, blocks=[block_u, p])
+
+        if out is not None:
+            out[0][0] = u
+            out[0][1] = ue
+            out[1]    = p
+            return out
+
+        block_u = BlockVector(self.domain.spaces[0], blocks=[u, ue])
+        return BlockVector(self.domain, blocks=[block_u, p])
+
+    def dot(self, b, out=None):
+        return self.solve(b, out=out)
\ No newline at end of file
diff --git a/pyproject.toml b/pyproject.toml
index 80cef7ef3..dafc3f469 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name            = "feectools"
-version         = "0.1.2"
+version         = "0.1.3"
 description     = "Slimmed-down fork of Psydac (https://github.com/pyccel/psydac) with less functionality and fewer dependencies."
 readme          = "README.md"
 requires-python = ">= 3.10"