Pipe bounds into Galerkin-in-time methods #157
Conversation
rckirby
left a comment
rckirby
left a comment
There was a problem hiding this comment.
This looks mostly good, but we either need to figure out how to fix continuous Galerkin with non-Bernstein or remove the failing tests (for now) and just add warnings when users try to do lagrange-type bases with bounds. If the Bernstein thing works now, we shouldn't let maths (not understanding failure of Lagrange) hold up a good feature.
| sub = Function(Vbig).assign(PETSc.INFINITY) | ||
|
|
||
| if bounds_type == "stage": | ||
| if (bounds_type == "stage") or (bounds_type == "galerkin"): |
There was a problem hiding this comment.
Not sure what bounds_type should read for Galerkin and DiscontinuousGalerkin. Maybe stage is fine tag, as there could be multiple ways of imposing Galerkin bounds (e.g. at all the dofs in time or just the last one).
There was a problem hiding this comment.
I really wasn't sure about this one. We could add a new tag or just work with what we already have. The tags stage and last_stage still make sense and seem to work as one would expect.
There was a problem hiding this comment.
last_stage is only going to make sense if you're in a basis nodal w.r.t. to the end points. But stage surely makes sense -- just apply the bounds to all the degrees of freedom.
| assert np.size(quadrature.get_points()) >= order | ||
|
|
||
| super().__init__(F, t, dt, u0, order, bcs=bcs, **kwargs) | ||
| if bounds is not None: |
There was a problem hiding this comment.
See above discussion in DG stepper for how this logic should go. Once we resolve that, make this the same.
|
|
||
| bc = BoundsConstrainedDirichletBC(V, uexact, "on_boundary", (lb, ub), solver_parameters=vi_params) | ||
|
|
||
| if quad_rule is not None: |
There was a problem hiding this comment.
This clobbers the argument quad_rule and uses it both as a function name and the result of calling that function.
| @@ -0,0 +1,328 @@ | |||
| import numpy as np | |||
There was a problem hiding this comment.
Several of the tests on CI are failing. Since we don't have good theory (yet!) for when these methods should pass/fail, I think we should be finding a set of tests that pass (e.g. just using Bernstein) for the CI. If bounds imposed on the GalerkinTimeStepper seem sufficiently fiddly that they frequently fail, maybe we want to allow non-Bernstein for "research purposes" but issue a stern warning that this is developmental/poorly understood and if the user wants something robust they should switch to Bernstein and/or DG?
There was a problem hiding this comment.
If you remove the Lagrange basis, all the tests pass. The heat equation is okay with DG + Lagrange. Interestingly, the wave equation is currently failing with DG + Lagrange for order 2 and higher--it appears to be a failed linear solve somewhere, not sure what that is about yet.
There was a problem hiding this comment.
Again, it could be the result of non-collocated DOFs and quadrature making the solver fail. It would help if we could exactly test that case for DG as well.
| return bounds_violations | ||
|
|
||
|
|
||
| @pytest.mark.parametrize('order', (1, 2, 3)) |
There was a problem hiding this comment.
I think this may be the issue -- I don't see us putting in the gll option anywhere into the Lagrange basis. This plus the same Gauss-Lobatto quadrature rule together is what I would hope might actually pass the test.
There was a problem hiding this comment.
Ah, to Pablo's point, to get Lagrange variants such as gll, we need to check for something like basis_type != "integral" when accepting bounds.
I'm still having trouble with convergence when using gll and the associated quadrature. The SNES gets stuck in a loop and runs into the MAX_IT.
There was a problem hiding this comment.
You're not testing GLL -- you should have basis_type = "gll" instead of just "Lagrange" to try that out. My conjecture is that this with GLL quadrature should work, Bernstein with any (accurate enough) quadrature should work, and that anything else might or might not. I would be mostly happy with this PR if that were the case.
There was a problem hiding this comment.
Sorry, I forgot to push the updated tests...new failing tests should be incoming.
Co-authored-by: Pablo Brubeck <brubeck@protonmail.com>
…oject/Irksome into jdstephens/DG_bounds
Initial work on piping bounds into Galerkin-in-time methods