-
Notifications
You must be signed in to change notification settings - Fork 11
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
Butcher Tableau #413
Butcher Tableau #413
Conversation
…/gusto into butcher_tableau
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This looks really neatly done to me. I think the main thing is to make sure that we're including the physics source evaluations (which were added in after you started this branch)
@@ -315,123 +331,55 @@ def apply(self, x_out, x_in): | |||
x_out.assign(self.x1) | |||
|
|||
|
|||
class ForwardEuler(ExplicitTimeDiscretisation): | |||
class ExplicitMultistage(ExplicitTimeDiscretisation): |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I don't mind ExplicitMultistage
, but have also wondered if ExplicitRungeKutta
is clearer?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think I may veto this.. I'm keener on multistage (especially in contrast to multistep)
self.k1.assign(self.x_out) | ||
self.x1.assign(x_in + 0.5 * self.dt * self.k1) | ||
def solve_stage(self, x0, stage): | ||
self.x1.assign(x0) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think we need to include the call to evaluate physics terms here (this was included after you branched off from trunk!). But you should check if makes sense to evaluate the physics sources at exactly this point...
self.x1.assign(x0) | |
self.x1.assign(x0) | |
for evaluate in self.evaluate_source: | |
evaluate(self.x1, self.dt) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I've added this!
self.x1.assign(x_in + 0.5 * self.dt * self.k2) | ||
if (stage == self.nStages - 1): | ||
self.x1.assign(x0) | ||
for i in range(self.nStages): |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I guess we need to add the physics source evaluation here too?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Not required here, as this is the final timestep evaluation of field (after solve)
super().__init__(domain, field_name=field_name, subcycles=subcycles, | ||
solver_parameters=solver_parameters, | ||
limiter=limiter, options=options, butcher_matrix=butcher_matrix) | ||
self.butcher_matrix = np.array([[1., 0., 0.], [1./4., 1./4., 0.], [1./6., 1./6., 2./3.]]) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I really like how clear this is now
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'm happy with this, but suggest including "Runge-Kutta" in the docstring for the new class
Co-authored-by: tommbendall <[email protected]>
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The only thing is to take care with line length, otherwise this looks great - thanks @atb1995 !
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Looks great - thanks @atb1995 !
We have created a new time stepper class
ExplicitMultistage
which takes in a Butcher Tableau in a numpy array of the form: [a_00 0 . 0 ][a_10 a_11 . 0 ]
[ . . . . ]
[ b_0 b_1 . b_{s-1}]
to run single time step multistage explicit time discretisation's. This allows time steppers such as
ForwardEuler
,RK4
,SSPRK3
andHeun
to be rewritten by the use of a single Butcher Tableau. It will also allow explicit multistage schemes to be easily written in future.Unfortunately two bugs were found in the process. The first is that the
Fallout
physics scheme has a few bugs in, meaningtest_precipitation.py
fails. I have xfailed this pytest for now. This problem is detailed in issue #334.The second was that
RK4
on main fails in some pytests when using theNonlinearVariationalSolver
. This failure translated to allExplicitMultistage
as they are written in the same way asRK4
. I have changed toLinearVariationalSolver
inExplicitMultistage
, and have created a minimal failing example outside of Gusto, which I will raise with Firedrake team.