BUG: linear solver fails if fields are initialised in a different order #4001

tommbendall · 2025-01-29T17:47:00Z

Describe the bug
Changing the order of declaration of a field leads to a solver failure. This is a Gusto-independent version of one of the several different errors currently breaking Gusto tests that most likely originated with Firedrake PR #3947

For the Gusto errors, see for instance: https://github.com/firedrakeproject/gusto/actions/runs/12864886794/job/36134821475?pr=610

Steps to Reproduce
Steps to reproduce the behavior:

Here is a failing example. Sorry that it's not more minimal, I couldn't find a smaller example that would fail!

This is clearly very odd. I can also make the test pass by making any number of other trivial changes, but changing the point at which the field N in the example below is declared seems to be enough. Otherwise the code below is a Gusto-independent version of the thermal shallow water solver which is currently exhibiting these failures.

"""
MFE for broken thermal shallow water solver.
"""

from firedrake import (
    IcosahedralSphereMesh, FunctionSpace, MixedFunctionSpace, Function,
    TestFunction, TestFunctions, TrialFunction, TrialFunctions, Constant,
    split, dx, dS, inner, dot, grad, div, jump, avg, VectorSpaceBasis,
    LinearVariationalProblem, LinearVariationalSolver, lhs, rhs, FacetNormal,
    SpatialCoordinate, VectorFunctionSpace, cross, CellNormal
)

# ============================================================================ #
# FLAG TO CONTROL WHETHER THIS TEST PASSES OR FAILS
set_up_normal_earlier = True   # False passes, True fails
# ============================================================================ #

# ---------------------------------------------------------------------------- #
# Set up core objects
# ---------------------------------------------------------------------------- #

degree = 1
radius = Constant(6371220.)
dt = Constant(3000.0)
H = Constant(3e4/9.80616)
g = Constant(9.80616)
Omega = Constant(7.292e-5)
alpha = Constant(0.5)
beta_u = alpha*dt
beta_d = alpha*dt

mesh = IcosahedralSphereMesh(float(radius), refinement_level=3, degree=2)
x = SpatialCoordinate(mesh)

# Function spaces
VDG = FunctionSpace(mesh, "DG", degree)
VHDiv = FunctionSpace(mesh, "BDM", degree+1)
mesh.init_cell_orientations(x)

VCG1 = FunctionSpace(mesh, 'CG', 1)
Vmixed = MixedFunctionSpace((VHDiv, VDG, VDG))
Vreduced = MixedFunctionSpace((VHDiv, VDG))

# Things for setting up perp
sphere_degree = mesh.coordinates.function_space().ufl_element().degree()
VecDG = VectorFunctionSpace(mesh, 'DG', sphere_degree)

# Mixed RHS and LHS functions
x_rhs = Function(Vmixed)
x_lhs = Function(Vmixed)

# Components of various mixed functions
u_rhs, D_rhs, b_rhs = split(x_rhs)
b_lhs = x_lhs.subfunctions[2]

uD_lhs = Function(Vreduced)
u_test, D_test = TestFunctions(Vreduced)
u_trial, D_trial = TrialFunctions(Vreduced)
u_lhsr, D_lhsr = uD_lhs.subfunctions

# Buoyancy trial/test functions
b_test = TestFunction(VDG)
b_trial = TrialFunction(VDG)

if set_up_normal_earlier:
    N = Function(VecDG).interpolate(CellNormal(mesh))

# Reference profiles
uDb_bar = Function(Vmixed)

if not set_up_normal_earlier:
    N = Function(VecDG).interpolate(CellNormal(mesh))

D_bar = split(uDb_bar)[1]
b_bar = split(uDb_bar)[2]
D_bar_subfunc = uDb_bar.subfunctions[1]
b_bar_subfunc = uDb_bar.subfunctions[2]

# Set up perp function
perp = lambda u: cross(N, u)

# Coriolis
f = Function(VCG1)

# ---------------------------------------------------------------------------- #
# Set up problem and solvers
# ---------------------------------------------------------------------------- #

# We have a 3-function mixed space for wind, depth and buoyancy
# Strategy is to approximately eliminate buoyancy, leaving a wind-depth problem
# The wind-depth problem is solved using the hydridization preconditioner

# Approximate elimination of b
b = -0.5 * dt * dot(u_trial, grad(b_bar)) + b_rhs

solver_parameters = {
    'ksp_type': 'preonly',
    'mat_type': 'matfree',
    'pc_type': 'python',
    'pc_python_type': 'firedrake.HybridizationPC',
    'hybridization': {
        'ksp_type': 'cg',
        'pc_type': 'gamg',
        'ksp_rtol': 1e-8,
        'mg_levels': {
            'ksp_type': 'chebyshev',
            'ksp_max_it': 2,
            'pc_type': 'bjacobi',
            'sub_pc_type': 'ilu'
        }
    }
}

n = FacetNormal(mesh)

# Linear thermal shallow water problem
uD_eqn = (
    inner(u_test, (u_trial - u_rhs)) * dx
    - beta_u * (D_trial - D_bar) * div(u_test * b_bar) * dx
    + beta_u * jump(u_test * b_bar, n) * avg(D_trial - D_bar) * dS
    - beta_u * 0.5 * D_bar * b_bar * div(u_test) * dx
    - beta_u * 0.5 * D_bar * b * div(u_test) * dx
    - beta_u * 0.5 * b_bar * div(u_test*(D_trial - D_bar)) * dx
    + beta_u * 0.5 * jump((D_trial - D_bar)*u_test, n) * avg(b_bar) * dS
    + inner(D_test, (D_trial - D_rhs)) * dx
    + beta_d * D_test * div(D_bar*u_trial) * dx
    + beta_u * f * inner(u_test, perp(u_trial)) * dx
)

# Boundary conditions
bcs = []

# Solver for u, D
uD_problem = LinearVariationalProblem(
    lhs(uD_eqn), rhs(uD_eqn), uD_lhs, bcs=bcs, constant_jacobian=True
)

# Provide callback for the nullspace of the trace system
def trace_nullsp(T):
    return VectorSpaceBasis(constant=True)

appctx = {"trace_nullspace": trace_nullsp}
uD_solver = LinearVariationalSolver(
    uD_problem, solver_parameters=solver_parameters, appctx=appctx
)

# Problem for reconstructing buoyancy
b_eqn = b_test*(b_trial - b_rhs + 0.5*dt*inner(u_lhsr, grad(b_bar))) * dx

b_problem = LinearVariationalProblem(
    lhs(b_eqn), rhs(b_eqn), b_lhs, constant_jacobian=True
)
b_solver = LinearVariationalSolver(b_problem)

# ---------------------------------------------------------------------------- #
# Set some initial conditions
# ---------------------------------------------------------------------------- #

f.interpolate(2*Omega*x[2]/radius)
D_bar_subfunc.assign(H)
b_bar_subfunc.assign(g)

# Simplest test is to give a right-hand side that is zero
x_rhs.assign(1.0)

# ---------------------------------------------------------------------------- #
# Run
# ---------------------------------------------------------------------------- #

print('wind-depth solve')
uD_solver.solve()
print('buoyancy solve')
b_solver.solve()
print('completed')

Expected behavior
I expect this to complete independent of the point at which the cell normal field N is first set.

Error message
A short back trace for the error is here:

Traceback (most recent call last):
  File "/home/thomas/odds_and_ends/thermal_swe_solver_mfe.py", line 172, in <module>
    uD_solver.solve()
  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 "/home/thomas/firedrake/src/firedrake/firedrake/adjoint_utils/variational_solver.py", line 104, in wrapper
    out = solve(self, **kwargs)
  File "/home/thomas/firedrake/src/firedrake/firedrake/variational_solver.py", line 325, in solve
    self.snes.solve(None, work)
  File "petsc4py/PETSc/SNES.pyx", line 1724, in petsc4py.PETSc.SNES.solve
petsc4py.PETSc.Error: error code -1
[0] SNESSolve() at /home/thomas/firedrake/src/petsc/src/snes/interface/snes.c:4839
[0] SNESSolve_KSPONLY() at /home/thomas/firedrake/src/petsc/src/snes/impls/ksponly/ksponly.c:49
[0] KSPSolve() at /home/thomas/firedrake/src/petsc/src/ksp/ksp/interface/itfunc.c:1075
[0] KSPSolve_Private() at /home/thomas/firedrake/src/petsc/src/ksp/ksp/interface/itfunc.c:900
[0] KSPSolve_PREONLY() at /home/thomas/firedrake/src/petsc/src/ksp/ksp/impls/preonly/preonly.c:27
[0] KSP_PCApply() at /home/thomas/firedrake/src/petsc/include/petsc/private/kspimpl.h:411
[0] PCApply() at /home/thomas/firedrake/src/petsc/src/ksp/pc/interface/precon.c:522

I will also attach the full backtrace, which prints out nearly 2000 lines of loopy code.

Environment:

OS: Linux (Ubuntu
Python version: Python 3.10.12
Output of firedrake-status:

Firedrake Configuration:
    package_manager: True
    minimal_petsc: False
    mpicc: None
    mpicxx: None
    mpif90: None
    mpiexec: None
    disable_ssh: False
    honour_petsc_dir: False
    with_parmetis: False
    slepc: False
    packages: ['git+ssh://github.com/firedrakeproject/gusto.git@main#egg=gusto']
    honour_pythonpath: False
    opencascade: False
    tinyasm: False
    torch: False
    petsc_int_type: int32
    cache_dir: /home/thomas/firedrake/.cache
    complex: False
    remove_build_files: False
    with_blas: None
    netgen: False
    jax: False
Additions:
    None
Environment:
    PYTHONPATH: None
    PETSC_ARCH: None
    PETSC_DIR: None
Status of components:
---------------------------------------------------------------------------
|Package             |Branch                        |Revision  |Modified  |
---------------------------------------------------------------------------
|FInAT_old           |master                        |d8836d5   |False     |
|PyOP2_old           |master                        |31471a60  |False     |
|fiat                |master                        |8839f876  |False     |
|firedrake           |master                        |58f3d6f26 |False     |
|gusto               |TBendall/FixCourant           |b1144cdc  |True      |
|h5py                |firedrake                     |06d2da80  |False     |
|libsupermesh        |master                        |f87cbfd   |True      |
|loopy               |main                          |ad07454c  |False     |
|petsc               |firedrake                     |4354988f8b6|False     |
|pyadjoint           |master                        |da58d45   |False     |
|pytest-mpi          |main                          |f5668e4   |False     |
|tsfc_old            |master                        |e06308d   |False     |
|ufl                 |master                        |5da0ab98  |False     |
---------------------------------------------------------------------------

Any relevant environment variables or modifications: Not that I know of

Additional Info
Unfortunately this is only one of several different errors currently afflicting Gusto since we picked up the changes from #3947

The text was updated successfully, but these errors were encountered:

tommbendall · 2025-01-29T17:49:10Z

thermal_swe_error.log

ksagiyam · 2025-01-30T02:50:10Z

I think the following reproduces the issue. So far I observed that it happens when:

we use lhs function, and
g.count() < c.count().

Example:

from firedrake import *


mesh = UnitTriangleMesh()
BDM = FunctionSpace(mesh, "BDM", 2)
DG = FunctionSpace(mesh, "DG", 1)
V = BDM * DG
VectorDG = VectorFunctionSpace(mesh, 'DG', 0)
u = TrialFunction(V)
v = TestFunction(V)
g = Function(VectorDG, count=10001)
c = Function(DG, count=10002)
a = (
    inner(u, v) * dx
    + inner(g[0] * u[0], v[0]) * dx
    + inner(c * grad(u[0]), grad(v[0])) * dx
)
sol = Function(V)
solver_parameters = {
    'mat_type': 'matfree',
    'ksp_type': 'preonly',
    'pc_type': 'python',
    'pc_python_type': 'firedrake.HybridizationPC',
}
# a is two-form.
# Does not fail.
solve(a == inner(as_vector([1., 1., 1.]), v) * dx, sol, solver_parameters=solver_parameters)
# Fails if g.count() < c.count().
solve(lhs(a) == rhs(a), sol, solver_parameters=solver_parameters)

Output:

  ...

  File "/home/sagiyama999/firedrake/src/loopy/loopy/check.py", line 776, in map_subscript
    raise LoopyIndexError("'%s' in instruction '%s' "
loopy.diagnostic.LoopyIndexError: 'w_1[2]' in instruction 'insn_65' accesses out-of-bounds array element (could not establish '{ [2] }' is a subset of '{ [i0] : 0 <= i0 <= 1 }').

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/home/sagiyama999/firedrake/src/firedrake/temp4.py", line 28, in <module>
    solve(lhs(a) == rhs(a), sol, solver_parameters=solver_parameters)
  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 "/home/sagiyama999/firedrake/src/firedrake/firedrake/adjoint_utils/solving.py", line 57, in wrapper
    output = solve(*args, **kwargs)
             ^^^^^^^^^^^^^^^^^^^^^^
  File "/home/sagiyama999/firedrake/src/firedrake/firedrake/solving.py", line 141, in solve
    _solve_varproblem(*args, **kwargs)
  File "/home/sagiyama999/firedrake/src/firedrake/firedrake/solving.py", line 178, in _solve_varproblem
    solver.solve()
  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 "/home/sagiyama999/firedrake/src/firedrake/firedrake/adjoint_utils/variational_solver.py", line 104, in wrapper
    out = solve(self, **kwargs)
          ^^^^^^^^^^^^^^^^^^^^^
  File "/home/sagiyama999/firedrake/src/firedrake/firedrake/variational_solver.py", line 325, in solve
    self.snes.solve(None, work)
  File "petsc4py/PETSc/SNES.pyx", line 1724, in petsc4py.PETSc.SNES.solve
petsc4py.PETSc.Error: error code -1
[0] SNESSolve() at /home/sagiyama999/firedrake/src/petsc/src/snes/interface/snes.c:4839
[0] SNESSolve_KSPONLY() at /home/sagiyama999/firedrake/src/petsc/src/snes/impls/ksponly/ksponly.c:49
[0] KSPSolve() at /home/sagiyama999/firedrake/src/petsc/src/ksp/ksp/interface/itfunc.c:1075
[0] KSPSolve_Private() at /home/sagiyama999/firedrake/src/petsc/src/ksp/ksp/interface/itfunc.c:900
[0] KSPSolve_PREONLY() at /home/sagiyama999/firedrake/src/petsc/src/ksp/ksp/impls/preonly/preonly.c:27
[0] KSP_PCApply() at /home/sagiyama999/firedrake/src/petsc/include/petsc/private/kspimpl.h:411
[0] PCApply() at /home/sagiyama999/firedrake/src/petsc/src/ksp/pc/interface/precon.c:522

tommbendall · 2025-01-30T08:33:42Z

Thanks @ksagiyam, that's a much simpler example!

pbrubeck · 2025-01-30T10:31:17Z

Maybe this is due to the fact that lhs and rhs make a call to expand_derivatives

pbrubeck · 2025-01-30T11:39:44Z

lhs(a) expands grad(u[0]) into grad(u)[0, :]

In [2]: print(a)
{ conj(((v_0) : (v_1))) } * dx(<Mesh #0>[everywhere], {})
  +  { v_1[0] * w⃗₁₀₀₀₁[0] * (conj((v_0[0]))) } * dx(<Mesh #0>[everywhere], {})
  +  { ({ A | A_{i_8} = (grad(v_1[0]))[i_8] * w₁₀₀₀₂ }) : (grad(v_0[0])) } * dx(<Mesh #0>[everywhere], {})

In [3]: print(lhs(a))
{ conj((sum_{i_{52}} v_0[i_{52}] * (conj((v_1[i_{52}]))) )) } * dx(<Mesh #0>[everywhere], {})
  +  { v_1[0] * w⃗₁₀₀₀₁[0] * (conj((v_0[0]))) } * dx(<Mesh #0>[everywhere], {})
  +  { sum_{i_{53}} ({ A | A_{i_8} = (grad(v_1))[0, i_8] * w₁₀₀₀₂ })[i_{53}] * (conj(((grad(v_0))[0, i_{53}])))  } * dx(<Mesh #0>[everywhere], {})

pbrubeck · 2025-01-30T11:50:14Z

This fails without lhs

from firedrake import *


mesh = UnitTriangleMesh()
BDM = FunctionSpace(mesh, "BDM", 2)
DG = FunctionSpace(mesh, "DG", 1)
V = BDM * DG
VectorDG = VectorFunctionSpace(mesh, 'DG', 0)
u = TrialFunction(V)
v = TestFunction(V)
g = Function(VectorDG, count=10001)
c = Function(DG, count=10002)
a = (
    inner(u, v) * dx
    + inner(g[0] * u[0], v[0]) * dx
    + inner(c * grad(u)[0, :], grad(v)[0, :]) * dx
)
sol = Function(V)
solver_parameters = {
    'mat_type': 'matfree',
    'ksp_type': 'preonly',
    'pc_type': 'python',
    'pc_python_type': 'firedrake.HybridizationPC',
}
# a is two-form.
# Fails
solve(a == inner(as_vector([1., 1., 1.]), v) * dx, sol, solver_parameters=solver_parameters)

pbrubeck · 2025-01-30T11:57:43Z

The issue also goes away if the TestFunction is split via TestFunctions

from firedrake import *


mesh = UnitTriangleMesh()
BDM = FunctionSpace(mesh, "BDM", 2)
DG = FunctionSpace(mesh, "DG", 1)
V = BDM * DG
VectorDG = VectorFunctionSpace(mesh, 'DG', 0)
u = TrialFunctions(V)
v = TestFunctions(V)
g = Function(VectorDG, count=10001)
c = Function(DG, count=10002)
a = (
     (inner(u[0], v[0]) + inner(u[1], v[1])) * dx
     + inner(g[0] * u[0][0], v[0][0]) * dx
    + inner(c * grad(u[0][0]), grad(v[0][0])) * dx
)
 sol = Function(V)
 solver_parameters = {
     'mat_type': 'matfree',
     'ksp_type': 'preonly',
     'pc_type': 'python',
     'pc_python_type': 'firedrake.HybridizationPC',
 }
 # a is two-form.
 # Does not fail.
 solve(a == inner(as_vector([1., 1.]), v[0]) * dx, sol, solver_parameters=solver_parameters)

solve(lhs(a) == rhs(a), sol, solver_parameters=solver_parameters)

pbrubeck · 2025-01-30T17:15:37Z

By inspecting @ksagiyam's MFE, with grad(v[0]) vs grad(v)[0, :], I noticed the following:

There is no mechanism to know that the unsplit mass term inner(u, v)*dx does not couple the BDM and the DG fields, so splitting a into submatrices gives nonzero offdiagonal blocks. In the case that passes this "zero" offdiagonal takes the form:

{ conj((([v_0[0], v_0[1], 0]) : ([0, 0, v_1]))) } * dx(<Mesh #31>[everywhere], {})

vs the failing case with

{ conj((([v_0[0], v_0[1], 0]) : ([0, 0, v_1]))) } * dx(<Mesh #0>[everywhere], {})
  +  { conj((({ A | A_{i_{10}} = (grad([v_0[0], v_0[1], 0]))[0, i_{10}] }) : ({ A | A_{i_9} = (grad([0, 0, v_1]))[0, i_9] * w₁₀₀₀₂ }))) } * dx(<Mesh #0>[everywhere], {})

This suggests that we should zero-simplify more eagerly. I'm not entirely sure whether the gusto test can be fixed in this way.

ksagiyam · 2025-01-30T17:26:07Z

That could work around the issue, but I think the code should work whether or not we zero-simplify.

pbrubeck · 2025-01-30T17:31:22Z

I agree. I also just want to point out that it could be that the zero-simplification that has been added so far is just exposing bugs that were not triggered by the rather complicated expressions that we previously allowed.

ksagiyam · 2025-01-30T17:34:53Z

Yes, I know.

ksagiyam · 2025-01-31T14:47:56Z

@tommbendall Can you see if #4003 fixes your issue?

tommbendall · 2025-01-31T16:47:44Z

This should fix about half of the failing Gusto tests -- and it does fix almost all of those. However our test_moist_gravity_wave_equivb_gw test (https://github.com/firedrakeproject/gusto/blob/main/examples/shallow_water/test_shallow_water_examples.py#L149) still seems to fail in the same way

ksagiyam · 2025-01-31T16:59:15Z

Can you show me the error message?

tommbendall · 2025-01-31T17:05:41Z

The brief part is this:

>       moist_thermal_gw(
            ncells_per_edge=4,
            dt=1800.,
            tmax=18000.,
            dumpfreq=10,
            dirname=make_dirname(test_name)
        )

examples/shallow_water/test_shallow_water_examples.py:152: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
examples/shallow_water/moist_thermal_equivb_gw.py:184: in moist_thermal_gw
    stepper.run(t=0, tmax=tmax)
gusto/timestepping/semi_implicit_quasi_newton.py:482: in run
    super().run(t, tmax, pick_up=pick_up)
gusto/timestepping/timestepper.py:216: in run
    self.timestep()
gusto/timestepping/semi_implicit_quasi_newton.py:433: in timestep
    self.linear_solver.solve(xrhs, dy)  # solves linear system and places result in dy
<decorator-gen-40>:2: in solve
    ???
../firedrake/pyop2/profiling.py:60: in wrapper
    return f(*args, **kwargs)
gusto/solvers/linear_solvers.py:758: in solve
    self.uD_solver.solve()
petsc4py/PETSc/Log.pyx:188: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
    ???
petsc4py/PETSc/Log.pyx:189: in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
    ???
../firedrake/firedrake/adjoint_utils/variational_solver.py:104: in wrapper
    out = solve(self, **kwargs)
../firedrake/firedrake/variational_solver.py:325: in solve
    self.snes.solve(None, work)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

>   ???
E   petsc4py.PETSc.Error: error code -1
E   [0] SNESSolve() at /home/thomas/firedrake/src/petsc/src/snes/interface/snes.c:4839
E   [0] SNESSolve_KSPONLY() at /home/thomas/firedrake/src/petsc/src/snes/impls/ksponly/ksponly.c:49
E   [0] KSPSolve() at /home/thomas/firedrake/src/petsc/src/ksp/ksp/interface/itfunc.c:1075
E   [0] KSPSolve_Private() at /home/thomas/firedrake/src/petsc/src/ksp/ksp/interface/itfunc.c:900
E   [0] KSPSolve_PREONLY() at /home/thomas/firedrake/src/petsc/src/ksp/ksp/impls/preonly/preonly.c:27
E   [0] KSP_PCApply() at /home/thomas/firedrake/src/petsc/include/petsc/private/kspimpl.h:411
E   [0] PCApply() at /home/thomas/firedrake/src/petsc/src/ksp/pc/interface/precon.c:522

petsc4py/PETSc/SNES.pyx:1724: Error
----------------------------- Captured stdout call -----------------------------

But I've attached the whole log in case that it also useful
moist_thermal_equivb_error.log

ksagiyam · 2025-01-31T22:09:31Z

If I install gusto and run ../gusto/examples/shallow_water/test_shallow_water_examples.py::test_moist_thermal_eqiuvb_gw, I get:

self = <firedrake.assemble.OneFormAssembler object at 0x7e04cf70f4a0>
tensor = Coefficient(WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x7e04d4ada840>, FiniteElement('Discontinuous Lagrange', triangle, 0), name=None), Mesh(VectorElement(FiniteElement('Lagrange', triangle, 2), dim=3), 4)), 432)

    def _check_tensor(self, tensor):
        if tensor.function_space() != self._form.arguments()[0].function_space().dual():
>           raise ValueError("Form's argument does not match provided result tensor")
E           ValueError: Form's argument does not match provided result tensor

Is this expected? You might want to fix these first.

jshipton · 2025-01-31T22:24:57Z

See this PR for fixes firedrakeproject/gusto#610 Sent from Outlook for iOS<https://aka.ms/o0ukef>

________________________________ From: ksagiyam ***@***.***> Sent: Friday, January 31, 2025 10:09:56 PM To: firedrakeproject/firedrake ***@***.***> Cc: Subscribed ***@***.***> Subject: Re: [firedrakeproject/firedrake] BUG: linear solver fails if fields are initialised in a different order (Issue #4001) CAUTION: This email originated from outside of the organisation. Do not click links or open attachments unless you recognise the sender and know the content is safe. If I install gusto and run ../gusto/examples/shallow_water/test_shallow_water_examples.py::test_moist_thermal_eqiuvb_gw, I get: self = <firedrake.assemble.OneFormAssembler object at 0x7e04cf70f4a0> tensor = Coefficient(WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x7e04d4ada840>, FiniteElement('Discontinuous Lagrange', triangle, 0), name=None), Mesh(VectorElement(FiniteElement('Lagrange', triangle, 2), dim=3), 4)), 432) def _check_tensor(self, tensor): if tensor.function_space() != self._form.arguments()[0].function_space().dual():

raise ValueError("Form's argument does not match provided result tensor")

E ValueError: Form's argument does not match provided result tensor Is this expected? You might want to fix these first. — Reply to this email directly, view it on GitHub<#4001 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ABUOTSVVSOKKX3FKDDE4FHL2NPYDJAVCNFSM6AAAAABWDNZU3OVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDMMRYGQ4DGNBSGI>. You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

ksagiyam · 2025-01-31T23:09:19Z

I do not see issues running pytest ../gusto/examples/shallow_water/test_shallow_water_examples.py if I use those branches.

---------------------------------------------------------------------------
|Package             |Branch                        |Revision  |Modified  |
---------------------------------------------------------------------------
|asQ                 |master                        |76cb14a   |False     |
|fiat                |master                        |8839f876  |False     |
|firedrake           |ksagiyam/fix_slate            |5aa2a35ee |False     |
|gusto               |TBendall/FixCourant           |b1144cdc  |False     |
|h5py                |firedrake                     |06d2da80  |False     |
|libsupermesh        |master                        |f87cbfd   |True      |
|libsupermesh_old    |master                        |1eccc99   |False     |
|loopy               |main                          |ad07454c  |False     |
|petsc               |firedrake                     |4354988f8b6|False     |
|pyadjoint           |master                        |da58d45   |False     |
|pytest-mpi          |main                          |f5668e4   |False     |
|ufl                 |master                        |5da0ab98  |False     |
---------------------------------------------------------------------------

tommbendall · 2025-02-01T09:45:24Z

Thanks @ksagiyam, I've just re-run that test with a clean Firedrake cache and it now passes, so I'd be happy that this is fixed by #4003

tommbendall added the bug label Jan 29, 2025

ksagiyam mentioned this issue Jan 31, 2025

slate: fix kernel argument ordering #4003

Merged

JHopeCollins closed this as completed in #4003 Feb 6, 2025

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

BUG: linear solver fails if fields are initialised in a different order #4001

BUG: linear solver fails if fields are initialised in a different order #4001

tommbendall commented Jan 29, 2025 •

edited

Loading

tommbendall commented Jan 29, 2025

ksagiyam commented Jan 30, 2025

tommbendall commented Jan 30, 2025

pbrubeck commented Jan 30, 2025

pbrubeck commented Jan 30, 2025 •

edited

Loading

pbrubeck commented Jan 30, 2025 •

edited

Loading

pbrubeck commented Jan 30, 2025 •

edited

Loading

pbrubeck commented Jan 30, 2025

ksagiyam commented Jan 30, 2025

pbrubeck commented Jan 30, 2025

ksagiyam commented Jan 30, 2025

ksagiyam commented Jan 31, 2025

tommbendall commented Jan 31, 2025

ksagiyam commented Jan 31, 2025

tommbendall commented Jan 31, 2025

ksagiyam commented Jan 31, 2025

jshipton commented Jan 31, 2025 via email

ksagiyam commented Jan 31, 2025

tommbendall commented Feb 1, 2025

BUG: linear solver fails if fields are initialised in a different order #4001

BUG: linear solver fails if fields are initialised in a different order #4001

Comments

tommbendall commented Jan 29, 2025 • edited Loading

tommbendall commented Jan 29, 2025

ksagiyam commented Jan 30, 2025

tommbendall commented Jan 30, 2025

pbrubeck commented Jan 30, 2025

pbrubeck commented Jan 30, 2025 • edited Loading

pbrubeck commented Jan 30, 2025 • edited Loading

pbrubeck commented Jan 30, 2025 • edited Loading

pbrubeck commented Jan 30, 2025

ksagiyam commented Jan 30, 2025

pbrubeck commented Jan 30, 2025

ksagiyam commented Jan 30, 2025

ksagiyam commented Jan 31, 2025

tommbendall commented Jan 31, 2025

ksagiyam commented Jan 31, 2025

tommbendall commented Jan 31, 2025

ksagiyam commented Jan 31, 2025

jshipton commented Jan 31, 2025 via email

ksagiyam commented Jan 31, 2025

tommbendall commented Feb 1, 2025

tommbendall commented Jan 29, 2025 •

edited

Loading

pbrubeck commented Jan 30, 2025 •

edited

Loading

pbrubeck commented Jan 30, 2025 •

edited

Loading

pbrubeck commented Jan 30, 2025 •

edited

Loading