From 0a1fd0eaa9af6c1a6be22f68e30477664f953b10 Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Thu, 29 Jan 2026 13:30:25 +0100 Subject: [PATCH 1/9] remove struphy test from PR workflow --- .github/workflows/test-PR-struphy.yml | 34 +++++++++++++++++---------- 1 file changed, 22 insertions(+), 12 deletions(-) diff --git a/.github/workflows/test-PR-struphy.yml b/.github/workflows/test-PR-struphy.yml index 221135a89..146d96a29 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,34 @@ 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 + 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 +116,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 From 6893e03c551aab4cb5e2d58168cafe0d065cc6bc Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Thu, 29 Jan 2026 13:31:04 +0100 Subject: [PATCH 2/9] delete struphy workflow on push --- .github/workflows/test-struphy.yml | 112 ----------------------------- 1 file changed, 112 deletions(-) delete mode 100644 .github/workflows/test-struphy.yml 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 From 8dd343323721569653b651154027a2f67c2b4d73 Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Thu, 29 Jan 2026 13:41:13 +0100 Subject: [PATCH 3/9] change struphy path for unit tests --- .github/workflows/test-PR-struphy.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/test-PR-struphy.yml b/.github/workflows/test-PR-struphy.yml index 146d96a29..97adbaeb9 100644 --- a/.github/workflows/test-PR-struphy.yml +++ b/.github/workflows/test-PR-struphy.yml @@ -97,6 +97,7 @@ jobs: 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 From 278908bdb513a402ae4121843f4887467b4a61b2 Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Thu, 29 Jan 2026 13:42:29 +0100 Subject: [PATCH 4/9] bump version number --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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" From 1981de121fe6949b4a0797a3d61c8089c25c0b9f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?D=C3=A1vid=20Sz=C3=A9kedi?= Date: Thu, 12 Mar 2026 22:21:50 +0000 Subject: [PATCH 5/9] Initial implementation of Uzawa algorithm. --- feectools/linalg/solvers.py | 131 +++++++++++++++++++++++++++++++++++- 1 file changed, 129 insertions(+), 2 deletions(-) diff --git a/feectools/linalg/solvers.py b/feectools/linalg/solvers.py index 31f084d6f..9428e85fc 100644 --- a/feectools/linalg/solvers.py +++ b/feectools/linalg/solvers.py @@ -4,7 +4,7 @@ """ 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 @@ -20,7 +20,7 @@ 'PBiConjugateGradientStabilized', 'MinimumResidual', 'LSMR', - 'GMRES' + 'GMRES', ) #=============================================================================== @@ -64,6 +64,7 @@ def inverse(A, solver, **kwargs): 'minres' : MinimumResidual, 'lsmr' : LSMR, 'gmres' : GMRES, + 'uzawa' : UzawaSolver } # Check solver input @@ -1912,3 +1913,129 @@ 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_solver='gmres', inner_tol=1e-10): # TODO specify these from outside + + self._options = { + "x0": x0, "tol": tol, "maxiter": maxiter, + "verbose": verbose, "recycle": recycle, + } + super().__init__(A, **self._options) + + self._A11 = A11 + self._A22 = A22 + self._B1 = B1 + self._B2 = B2 + + # inner solves for A11^{-1} and A22^{-1} + self._A11inv = inverse(A11, inner_solver, tol=inner_tol, maxiter=maxiter, verbose=False) + self._A22inv = inverse(A22, inner_solver, tol=inner_tol, maxiter=maxiter, verbose=False) + + # pre-allocate temporaries + self._tmps_u = A11.domain.zeros() + self._tmps_ue = A22.domain.zeros() + self._tmps_p = B1.codomain.zeros() + + self._info = None + + def solve(self, b, out=None): + """ + Uzawa iteration on the saddle-point system. + + b is a BlockVector [F, g] where F = [f_u, f_ue] and g is the + constraint RHS. + """ + + 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) + rhs_u -= A11.dot(u) + u += A11inv.dot(rhs_u) + + # solve A22 * ue = f_ue - B2^T * p + rhs_ue = f_ue - B2.T.dot(p) + rhs_ue -= A22.dot(ue) + 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 on Schur complement) + 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: + # store solution as next initial guess + from feectools.linalg.block import BlockVector, BlockVectorSpace + block_u = BlockVector(BlockVectorSpace(A11.domain, A22.domain), blocks=[u, ue]) + self._options["x0"] = BlockVector(self.domain, blocks=[block_u, p]) + + # pack solution back into BlockVector matching the full system + 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]) + sol = BlockVector(self.domain, blocks=[block_u, p]) + return sol + + def dot(self, b, out=None): + return self.solve(b, out=out) \ No newline at end of file From d2a48ef19d31ae22ba238369cd5a53c679c6f1d8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?D=C3=A1vid=20Sz=C3=A9kedi?= Date: Mon, 16 Mar 2026 14:07:25 +0000 Subject: [PATCH 6/9] added todos --- feectools/linalg/solvers.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/feectools/linalg/solvers.py b/feectools/linalg/solvers.py index 9428e85fc..b0c36e814 100644 --- a/feectools/linalg/solvers.py +++ b/feectools/linalg/solvers.py @@ -1931,7 +1931,7 @@ def __init__(self, A, *, A11, A22, B1, B2, self._B2 = B2 # inner solves for A11^{-1} and A22^{-1} - self._A11inv = inverse(A11, inner_solver, tol=inner_tol, maxiter=maxiter, verbose=False) + self._A11inv = inverse(A11, inner_solver, tol=inner_tol, maxiter=maxiter, verbose=False) # TODO option to use direct solve self._A22inv = inverse(A22, inner_solver, tol=inner_tol, maxiter=maxiter, verbose=False) # pre-allocate temporaries @@ -2006,7 +2006,7 @@ def solve(self, b, out=None): if residual_norm < tol: break - # pressure update (steepest descent on Schur complement) + # pressure update (steepest descent on Schur complement) TODO preconditioning 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 @@ -2026,7 +2026,6 @@ def solve(self, b, out=None): block_u = BlockVector(BlockVectorSpace(A11.domain, A22.domain), blocks=[u, ue]) self._options["x0"] = BlockVector(self.domain, blocks=[block_u, p]) - # pack solution back into BlockVector matching the full system if out is not None: out[0][0] = u out[0][1] = ue From ae6859fcb16c765f7bb7cabb82eb6268d1967014 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?D=C3=A1vid=20Sz=C3=A9kedi?= Date: Sun, 19 Apr 2026 17:32:17 +0000 Subject: [PATCH 7/9] Fixed a small bug. --- feectools/linalg/basic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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() From a47ee5249f5e8885acac920e5461cd3cb6ebc669 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?D=C3=A1vid=20Sz=C3=A9kedi?= Date: Tue, 7 Jul 2026 16:22:15 +0000 Subject: [PATCH 8/9] Changed Uzawa to use direct inner solves. --- feectools/linalg/solvers.py | 64 +++++++++++++++++-------------------- 1 file changed, 30 insertions(+), 34 deletions(-) diff --git a/feectools/linalg/solvers.py b/feectools/linalg/solvers.py index b0c36e814..648cacb85 100644 --- a/feectools/linalg/solvers.py +++ b/feectools/linalg/solvers.py @@ -10,6 +10,8 @@ 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', @@ -1916,8 +1918,7 @@ def dot(self, b, out=None): class UzawaSolver(InverseLinearOperator): def __init__(self, A, *, A11, A22, B1, B2, - x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle=False, - inner_solver='gmres', inner_tol=1e-10): # TODO specify these from outside + x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle=False): self._options = { "x0": x0, "tol": tol, "maxiter": maxiter, @@ -1930,29 +1931,29 @@ def __init__(self, A, *, A11, A22, B1, B2, self._B1 = B1 self._B2 = B2 - # inner solves for A11^{-1} and A22^{-1} - self._A11inv = inverse(A11, inner_solver, tol=inner_tol, maxiter=maxiter, verbose=False) # TODO option to use direct solve - self._A22inv = inverse(A22, inner_solver, tol=inner_tol, maxiter=maxiter, verbose=False) - - # pre-allocate temporaries - self._tmps_u = A11.domain.zeros() - self._tmps_ue = A22.domain.zeros() - self._tmps_p = B1.codomain.zeros() + self._A11inv = self._factorize(A11) + self._A22inv = self._factorize(A22) self._info = None - def solve(self, b, out=None): - """ - Uzawa iteration on the saddle-point system. + @staticmethod + def _factorize(A): + """Compute sparse LU factorization of A.""" + from feectools.linalg.direct_solvers import SparseSolver + spmat = A.tosparse() + return SparseSolver(spmat) - b is a BlockVector [F, g] where F = [f_u, f_ue] and g is the - constraint RHS. - """ + def update_A11(self, A11): + """Recompute A11 factorization when dt changes.""" + self._A11 = A11 + self._A11inv = self._factorize(A11) + + def solve(self, b, out=None): - A11 = self._A11 - A22 = self._A22 - B1 = self._B1 - B2 = self._B2 + A11 = self._A11 + A22 = self._A22 + B1 = self._B1 + B2 = self._B2 A11inv = self._A11inv A22inv = self._A22inv @@ -1961,8 +1962,8 @@ def solve(self, b, out=None): verbose = self._options["verbose"] recycle = self._options["recycle"] - F = b[0] - g = b[1] + F = b[0] + g = b[1] f_u = F[0] f_ue = F[1] @@ -1986,15 +1987,13 @@ def solve(self, b, out=None): for iteration in range(1, maxiter + 1): - # solve A11 * u = f_u - B1^T * p + # solve A11 * u = f_u - B1^T * p (exact direct solve) rhs_u = f_u - B1.T.dot(p) - rhs_u -= A11.dot(u) - u += A11inv.dot(rhs_u) + u = A11inv.solve(rhs_u.toarray()) - # solve A22 * ue = f_ue - B2^T * p + # solve A22 * ue = f_ue - B2^T * p (exact direct solve) rhs_ue = f_ue - B2.T.dot(p) - rhs_ue -= A22.dot(ue) - ue += A22inv.dot(rhs_ue) + ue = A22inv.solve(rhs_ue.toarray()) # constraint residual: R = B1*u + B2*ue - g R = B1.dot(u) + B2.dot(ue) - g @@ -2006,8 +2005,8 @@ def solve(self, b, out=None): if residual_norm < tol: break - # pressure update (steepest descent on Schur complement) TODO preconditioning - S_R = B1.dot(A11inv.dot(B1.T.dot(R))) + B2.dot(A22inv.dot(B2.T.dot(R))) + # pressure update: steepest descent step size + S_R = B1.dot(A11inv.solve(B1.T.dot(R).toarray())) + B2.dot(A22inv.solve(B2.T.dot(R).toarray())) alpha = R.inner(R).real / R.inner(S_R).real p += alpha * R @@ -2021,8 +2020,6 @@ def solve(self, b, out=None): } if recycle: - # store solution as next initial guess - from feectools.linalg.block import BlockVector, BlockVectorSpace block_u = BlockVector(BlockVectorSpace(A11.domain, A22.domain), blocks=[u, ue]) self._options["x0"] = BlockVector(self.domain, blocks=[block_u, p]) @@ -2033,8 +2030,7 @@ def solve(self, b, out=None): return out block_u = BlockVector(BlockVectorSpace(A11.domain, A22.domain), blocks=[u, ue]) - sol = BlockVector(self.domain, blocks=[block_u, p]) - return sol + 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 From 82151da2b0b7a31c5fc22c999eef7dfadda87d7a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?D=C3=A1vid=20Sz=C3=A9kedi?= Date: Tue, 14 Jul 2026 12:39:21 +0000 Subject: [PATCH 9/9] Added second saddle point solver from Restelli paper. --- feectools/linalg/solvers.py | 144 ++++++++++++++++++++++++++++++------ 1 file changed, 121 insertions(+), 23 deletions(-) diff --git a/feectools/linalg/solvers.py b/feectools/linalg/solvers.py index 648cacb85..dbbac8e58 100644 --- a/feectools/linalg/solvers.py +++ b/feectools/linalg/solvers.py @@ -66,7 +66,8 @@ def inverse(A, solver, **kwargs): 'minres' : MinimumResidual, 'lsmr' : LSMR, 'gmres' : GMRES, - 'uzawa' : UzawaSolver + 'uzawa' : UzawaSolver, + 'schur' : SchurComplementSolver } # Check solver input @@ -1918,7 +1919,9 @@ def dot(self, b, out=None): class UzawaSolver(InverseLinearOperator): def __init__(self, A, *, A11, A22, B1, B2, - x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle=False): + 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, @@ -1926,27 +1929,29 @@ def __init__(self, A, *, A11, A22, B1, B2, } 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._B1 = B1 + self._B2 = B2 - self._A11inv = self._factorize(A11) - self._A22inv = self._factorize(A22) + 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 _factorize(A): - """Compute sparse LU factorization of A.""" - from feectools.linalg.direct_solvers import SparseSolver - spmat = A.tosparse() - return SparseSolver(spmat) + 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 factorization when dt changes.""" + """Recompute A11 inverse when dt changes.""" self._A11 = A11 - self._A11inv = self._factorize(A11) + self._A11inv = self._inner_solve(A11, tol=self._inner_tol, maxiter=self._inner_maxiter) def solve(self, b, out=None): @@ -1978,6 +1983,7 @@ def solve(self, b, out=None): ue = A22.domain.zeros() p = B1.codomain.zeros() + if verbose: print("Uzawa solver:") print("+---------+---------------------+") @@ -1987,13 +1993,13 @@ def solve(self, b, out=None): for iteration in range(1, maxiter + 1): - # solve A11 * u = f_u - B1^T * p (exact direct solve) + # solve A11 * u = f_u - B1^T * p rhs_u = f_u - B1.T.dot(p) - u = A11inv.solve(rhs_u.toarray()) + u = A11inv.dot(rhs_u) - # solve A22 * ue = f_ue - B2^T * p (exact direct solve) + # solve A22 * ue = f_ue - B2^T * p rhs_ue = f_ue - B2.T.dot(p) - ue = A22inv.solve(rhs_ue.toarray()) + ue = A22inv.dot(rhs_ue) # constraint residual: R = B1*u + B2*ue - g R = B1.dot(u) + B2.dot(ue) - g @@ -2006,31 +2012,123 @@ def solve(self, b, out=None): break # pressure update: steepest descent step size - S_R = B1.dot(A11inv.solve(B1.T.dot(R).toarray())) + B2.dot(A22inv.solve(B2.T.dot(R).toarray())) + 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 + p += alpha * R if verbose: print("+---------+---------------------+") self._info = { - 'niter': iteration, - 'success': residual_norm < tol, + 'niter' : iteration, + 'success' : residual_norm < tol, 'res_norm': residual_norm, } if recycle: - block_u = BlockVector(BlockVectorSpace(A11.domain, A22.domain), blocks=[u, ue]) + 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 + 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