Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

BUG: LinearSolver.solve fails when passed two Functions. #3203

Closed
APaganini opened this issue Nov 3, 2023 · 8 comments · May be fixed by #3922
Closed

BUG: LinearSolver.solve fails when passed two Functions. #3203

APaganini opened this issue Nov 3, 2023 · 8 comments · May be fixed by #3922
Labels

Comments

@APaganini
Copy link
Contributor

APaganini commented Nov 3, 2023

Describe the bug
LinearSolver.solve fails when passed two Functions. The issue may be that here b could be a function (which is allowed in solve, see here), but blift returned by _rhs (see here) is always a cofunction.

Steps to Reproduce

from firedrake import *

mesh = UnitSquareMesh(3,3)
V = FunctionSpace(mesh, "CG", 1)
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u),grad(v))*dx
bcs = DirichletBC(V, 0, "on_boundary")
A = assemble(a, mat_type='aij', bcs=bcs)
L = LinearSolver(A)

uout = Function(V)
uin = Function(V)
L.solve(uout, uin)

Error message

Traceback (most recent call last):
 File "/Users/alberto/Documents/FIREDRAKE/fireshape/tests/secondminfail.py", line 14, in <module>
  L.solve(uout, uin)
 File "petsc4py/PETSc/Log.pyx", line 188, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
 File "petsc4py/PETSc/Log.pyx", line 189, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
 File "/Users/alberto/Documents/FIREDRAKE/firedrake/src/firedrake/firedrake/linear_solver.py", line 163, in solve
  b = self._lifted(b)
    ^^^^^^^^^^^^^^^
 File "/Users/alberto/Documents/FIREDRAKE/firedrake/src/firedrake/firedrake/linear_solver.py", line 133, in _lifted
  blift += b
TypeError: unsupported operand type(s) for +=: 'Cofunction' and 'Function'

Environment:
MacOS,
Python 3.11.6

Firedrake Configuration:

    package_manager: True
    minimal_petsc: False
    mpicc: None
    mpicxx: None
    mpif90: None
    mpiexec: None
    disable_ssh: True
    honour_petsc_dir: False
    with_parmetis: False
    slepc: False
    packages: []
    honour_pythonpath: False
    opencascade: False
    tinyasm: False
    torch: False
    petsc_int_type: int32
    cache_dir: /Users/alberto/Documents/FIREDRAKE/firedrake/.cache
    complex: False
    remove_build_files: False
    with_blas: download
    netgen: False
Additions:
    None
Environment:
    PYTHONPATH: None
    PETSC_ARCH: None
    PETSC_DIR: None
Status of components:
---------------------------------------------------------------------------
|Package             |Branch                        |Revision  |Modified  |
---------------------------------------------------------------------------
|FInAT               |master                        |f493084   |False     |
|PyOP2               |master                        |d017d594  |False     |
|fiat                |master                        |e440a06   |False     |
|firedrake           |master                        |1692a70c7 |False     |
|h5py                |firedrake                     |6cc4c912  |False     |
|libspatialindex     |master                        |4768bf3   |True      |
|libsupermesh        |master                        |dbe226b   |False     |
|loopy               |main                          |8158afdb  |False     |
|petsc               |firedrake                     |98a9995288|False     |
|pyadjoint           |master                        |cefba3f   |False     |
|pytest-mpi          |main                          |a478bc8   |False     |
|tsfc                |master                        |47c1324   |False     |
|ufl                 |master                        |cf66c4b3  |False     |
---------------------------------------------------------------------------
@APaganini APaganini added the bug label Nov 3, 2023
@APaganini APaganini changed the title BUG: LinearSolver.solve fails when passed to Functions. BUG: LinearSolver.solve fails when passed two Functions. Nov 5, 2023
@connorjward
Copy link
Contributor

@nbouziani you added the previous dual code to LinearSolver. Is it possible that you just missed a case?

@nbouziani
Copy link
Contributor

I am not sure it makes sense for both to be Functions. I think the rhs needs to be a Cofunction, which can either be the result of assembling a 1-form or a Cofunction specified by the user. For instance, your example aims at solving $Au_{out} = u_{in}$, but the action of the 2-form a on uout is a 1-form and therefore lives in $V^{*}$, so uin should be a Cofunction.

@wence-
Copy link
Contributor

wence- commented Nov 8, 2023

Yeah, $A : V \to U$, so $u_{\text{in}}$ should live in in $U$, not $V$. That isn't to say it should by necessity be a cofunction. That solver code should probably be redone to check that $u_\text{out}$ lives in the trial space and $u_\text{in}$ lives in the dual of the test space.

@APaganini
Copy link
Contributor Author

APaganini commented Nov 8, 2023

I agree with the theory, but I'm unsure about the practice: how strongly is the cofunction business enforced? linear_solver allows passing functions (see here), which leads to the error being raised in the wrong place. And looking at the manual doesn't help either: in the section "Solving linear systems" it says that "In Firedrake, A is represented as a Matrix, while b and x are both Function".

Update: indeed replacing uin = Function(V) with uin = Cofunction(V.dual()) make the error disappear, but I had to read the source code of assemble to find out how to define a cofunction...

to sum up, it looks like the bug is that linear_solve still accepts Functions where it shouldn't

@wence-
Copy link
Contributor

wence- commented Nov 8, 2023

how strongly is the cofunction business in enforced?

Places where it is not enforced correctly are bugs.

@nbouziani
Copy link
Contributor

to sum up, it looks like the bug is that linear_solve still accepts Functions where it shouldn't

The linear solver can accept Functions and Cofunctions, which would be the case for your problem. The bug is that checking that the solution $u_{out}$ and the rhs $u_{in}$ are functions or cofunctions is not enough, their spaces need to be compatible with the matrix A as well. For $A : V -&gt; U$, we must have $u_{out} \in V$ and $u_{in} \in U$. Now the actual type of both depends on whether U and V are primal or dual spaces.

As for the manual, this part seems out-of-date and should be updated. I will open up a PR to fix both issues.

@nbouziani
Copy link
Contributor

See #3214

@APaganini
Copy link
Contributor Author

I guess this can be closed.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging a pull request may close this issue.

4 participants