Skip to content

Support blocked linear systems. - #73

Merged
finsberg merged 27 commits into
mainfrom
dokken/blocked-linear
Aug 28, 2026
Merged

Support blocked linear systems.#73
finsberg merged 27 commits into
mainfrom
dokken/blocked-linear

Conversation

@jorgensd

@jorgensd jorgensd commented Aug 26, 2026

Copy link
Copy Markdown
Member

Summary

Adds support for blocked forms [[a00, ..., a0N], ...., [aN0, ..., aNN]] as input to LinearProblem.

Bug-fixes for other problems

  • For the time-dependent BC use a direct solver to ensure we solve to machine precision and are independent of initial guess.

Various lessons learnt along the way

  • Using compute_adjoint action on problems that mixes ufl.MixedFunctionSpace arguments and non-mixed arguments is a bad idea (looking at you prepare_hessian). Therefore we use the trick of creating a Lagrangian and differentiate it.
  • Various nested type-hinting to make the recursive code mypy friendly.
  • Action on a form does not like lists, only place that has a bit of split handling at the momment (_compute_residual).
    -As pyadjoint eagerly creates new functions, we have to recompile the forms (maybe we should use the data independent form creators, would align with redoxnics ideas).
  • Add convenience function to get arguments sorted by part.
  • Minor bug-fixes on using PETSc/non-PETSc vectors.
  • Correct usage of idx in evaluate_hessian_component. idx only relates to inputs. Output is in control space (always). Added check of this.
  • Instead of recreating vectors to attach to the tlm/adjoint solvers, instead use the vector that is created once and zero it out. Better for evyerone.

jorgensd and others added 27 commits August 25, 2026 14:06
After long dicussion with Gemini (August 2026), I disregarded all advice it gave, and figured out that one of the two issues were:
- Since we use SNES and setCTX in the setFunction and setJacobian operations, we can't use the classical python-esque replacement of compiled forms. This invalidated the recompute for non-linear problems.

Another issue covered by Claude (in another PR spawned by @finsberg) is that `u` was added as its own dependency if the problem was truely nonlinear. This is now fixed.
Added tests to check this.
…f my own thoughts/changes:

-As pyadjoint eagerly creates new functions, we have to recompile the forms (maybe we should use the data independent form creators, would align with redoxnics ideas).
- Add convenience function to get arguments sorted by part.
- Minor bug-fixes on using PETSc/non-PETSc vectors.
- Correct usage of idx in evaluate_hessian_component. idx only relates to inputs. Output is in control space (always). Added check of this.
- Instead of recreating vectors to attach to the tlm/adjoint solvers, instead use the vector that is created once and zero it out. Better for evyerone.
- Using compute_adjoint action on problems that mixes ufl.MixedFunctionSpace arguments and non-mixed arguments is a bad idea (looking at you prepare_hessian). Therefore we use the trick of creating a Lagrangian and differentiate it.
- Various nested type-hinting to make the recursive code mypy friendly.
- Action on a form does not like lists, only place that has a bit of split handling at the momment (_compute_residual).
Simplify the code for adjoint in linear solver.
Time-distributed control + blocked rewrite.
…o let adjoint_solver have its own internal state and then copy to correct array post solve.
@finsberg
finsberg self-requested a review August 28, 2026 07:31
@finsberg
finsberg merged commit 23d8042 into main Aug 28, 2026
5 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants