DAE support in StageDerivativeTimeStepper - #251
Conversation
| try: | ||
| stage_update = numpy.linalg.solve(butcher_tableau.A.T, butcher_tableau.b) | ||
| self.aux_updateb = vecconst(numpy.linalg.solve(A2.T, stage_update), backend=backend) | ||
| self.aux_update_scale = 1 - numpy.sum(stage_update) |
There was a problem hiding this comment.
I didn't have to do the maths, codex figured out that the correct linear combination for the update should match that of the stage-value formulation.
There was a problem hiding this comment.
I recall there being problems with updates in the non-stiffly accurate case. I suspect this may still fail silently and compute the wrong answer if you have things not correctly labeled as auxiliary. So some tests might be in order.
| stage_update = numpy.linalg.solve(butcher_tableau.A.T, butcher_tableau.b) | ||
| self.aux_updateb = vecconst(numpy.linalg.solve(A2.T, stage_update), backend=backend) |
There was a problem hiding this comment.
However, given that A = A1 * A2, I think this could be simplified into a single linear system solve involving A1
| stage_update = numpy.linalg.solve(butcher_tableau.A.T, butcher_tableau.b) | |
| self.aux_updateb = vecconst(numpy.linalg.solve(A2.T, stage_update), backend=backend) | |
| stage_update = numpy.linalg.solve(butcher_tableau.A1.T, butcher_tableau.b) | |
| self.aux_updateb = vecconst(stage_update, backend=backend) |
Something along those lines
There was a problem hiding this comment.
Please confirm and get that correct?
| :kwarg aux_indices: Valid for continuous Petrov Galerkin schemes and | ||
| derivative-stage RK schemes. It specifies that some variables in `u0` |
There was a problem hiding this comment.
| :kwarg aux_indices: Valid for continuous Petrov Galerkin schemes and | |
| derivative-stage RK schemes. It specifies that some variables in `u0` | |
| :kwarg aux_indices: Valid for stage-derivative RK and continuous Petrov | |
| Galerkin schemes. It specifies that some variables in `u0` |
|
This should be a no-op for a stiffly accurate scheme, shouldn't it? What's the motivation for pushing things through for the GL (or other non-stiffly accurate) case? |
ScottMacLachlan
left a comment
ScottMacLachlan
left a comment
There was a problem hiding this comment.
This does not make sense to me.
| u0bit += sum(self.stages.subfunctions[nf * s + i] * (self.aux_updateb[s] * dt) | ||
| for s in range(ns)) | ||
| else: | ||
| u0bit += sum(self.stages.subfunctions[nf * s + i] * (b[s] * dt) for s in range(ns)) |
There was a problem hiding this comment.
This introduces an if statement into the general _update formula which is only needed for an edge case. I would argue this deserves a special update function that is selected at __init__ time and not foisted on every other use case
There was a problem hiding this comment.
It is what we have in Galerkin, happy to split it into its own function
There was a problem hiding this comment.
It's probably a minor ding on run-time given the cost of solving systems, but if we can lift if statements out of do loops, I'm in favor.
It's a partitioned method. Stage-derivative for differential variables and stage-value for algebraic ones |
GiT with integrated Lagrange trial functions and Lagrange test functions becomes stage-derivative RK if you underintegrate correctly. GiT may do structure-preservation, but one wants (say, for NSE) pressure in the auxiliary test/DG space to make this happen. So, the motivation here would be to allow IRK to emulate the Right Thing for Galerkin-in-time, and presumably to allow preconditioning? |
The motivation is to allow IRK to do the Right Thing for IRK. Secondly to allow preconditioning. |
IRK already does the right thing for IRK, though. IRK for DAEs is perfectly well-defined, and it's what we already have implemented. If there's a use case for this variant, I have no problem with it being added, but I think it's better as an independent |
Upon reflection: no. I don't know what this is, but it isn't a partitioned method. (Or, at least, if it is, it's only by some coincidental backdoor that isn't clear.) To do a partitioned method, you'd need changes to |
I don't know what the right thing is here. A reference to the literature might be in order saying "we're implementing the method of X" |
rckirby
left a comment
rckirby
left a comment
There was a problem hiding this comment.
See comments already made. Also, please put an actual description of what problem you're fixing/feature you're adding in the PR, as it will make getting started on reviewing much easier!
3 participants
No description provided.