Skip to content
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

Cory Adams Moulton (V2) #10

Open
wants to merge 23 commits into
base: main
Choose a base branch
from
Open

Conversation

cmikida2
Copy link
Contributor

This pull request takes the single-rate component of the implicit Adams development associated with this pull request and splits it into a separate PR, while also attempting to fix a number of documentation and refactoring issues raised in the initial PR. Among the largest changes are

  • Adding a base class for Adams methods to remove duplicate code.
  • Adding a number of utilities to that base class to remove duplicate code, largely associated with history rotation.
  • Moving the implicit solve used by both RK and Adams methods to implicit.py.
  • Refactoring the generate_butcher function in the RK base class to allow for Adams methods to use this function for bootstrapping.
  • Modifying the simple convergence testing utility to allow for implicit testing, and folding Adams-Moulton testing into test_ab.py (test_ab.py -> test_adams.py).
  • Adding a lot of documentation, but I am still iterating on this and trying to get this part exactly right.

@cmikida2 cmikida2 mentioned this pull request Oct 27, 2020
@cmikida2 cmikida2 marked this pull request as draft November 2, 2020 18:05
@cmikida2 cmikida2 marked this pull request as ready for review November 2, 2020 18:06
leap/implicit.py Outdated

def generate_solve(cb, unknowns, equations, var_to_unknown, guess):

# got a square system, let's solve
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
# got a square system, let's solve
assert len(unknowns) == len(equations)

?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This change has been added in ae9d04f.

recycle_last_stage_coeff_set_names = ()


class DIRK2MethodBuilder(ImplicitButcherTableauMethodBuilder):
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are these methods being tested somewhere?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is addressed in 6d93fee, which adds accuracy tests for these methods to the RK testing suite.

@@ -219,7 +222,7 @@ class AdamsBashforthMethodBuilder(MethodBuilder):
"""

def __init__(self, component_id, function_family=None, state_filter_name=None,
hist_length=None, static_dt=False, order=None):
hist_length=None, static_dt=False, order=None, _extra_bootstrap=False):
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not liking this _extra_bootstrap thing, mostly because it's ill-defined.

  • What happens with the extra history after bootstrap?
  • What are its semantics in the explicit case?
  • The justification for it seems really weak and mostly focused on debugging?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is an option that was previously in place to pass a certain test that compared multi-rate IMEX to single-rate implicit. Given that that's not what we're talking about here, this is an artifact of outdated development and needs to be removed. I'm on it.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This option and its artifacts are removed by ae9d04f.

leap/multistep/__init__.py Show resolved Hide resolved

# {{{ emit solve if possible

if unknowns and len(unknowns) == len(equations):
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is the else to this if? Why is it OK to just go by here if we can't solve?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The else - as far as I'm aware - would be the fully explicit case, but at bare minimum in that case the if condition isn't indicative of that.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmm. I realize now, though, that I only have this in the Adams-Moulton class, so explicit isn't a concern. Let me think about contingencies/what should be done if we don't have a square system, or if it's even possible for that to happen with the current setup.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Addressed by 3803398.

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Non-square systems aren't allowed IMO.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When you say "not allowed" do you mean they should never happen? I added the elifs (and now an else) to this block to catch instances in which a system is non-square, but I'm not sure if you're saying that that's not sufficient, or if the check for a square system isn't needed at all because we should never not have one.

Comment on lines 572 to 581
if self.extra_bootstrap:
save_crit = i+1
else:
save_crit = i

with cb.if_(self.step, "==", save_crit):
cb(self.history[i], rhs_next_var)

if not self.static_dt:
cb(self.time_history[i], self.t + self.dt)
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This code seems like it's an adverse consequence of the awkward redefinition of "history".

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That and the extra bootstrapping. Given that it's a function of those two things, I'll need to think about how those things need to change before addressing this.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is addressed along with redefinition of history by 4a87ba9.

@cmikida2
Copy link
Contributor Author

I've addressed all of your comments above and made a number of changes to the Adams-Moulton code, and while most of them should explain themselves as direct responses to the comments above, I want to also mention that I had to very slightly twiddle the convergence tolerance for the accuracy tests in order for the implicit extended-history tests to pass. The pass criteria was changed from 0.9 * expected_order to 0.89 * expected_order, as the EOCs for these cases were one or two hundredths places away from the initial target. The one-line code change for this is in 4a87ba9.

I have checked that the addition of extended history has the same effect on EOC relative to the square system as in the case of Adams-Bashforth - in both cases, there's a slight drop in EOC, typically on the order of 0.08-0.1. I'm not overly concerned about this, but happy to discuss further, as I don't think you were all too thrilled with the fudge factor that was being applied to begin with.

Copy link
Owner

@inducer inducer left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks! Once the wrinkles below are addressed, I think this is in good shape.

leap/implicit.py Outdated Show resolved Hide resolved
Comment on lines 351 to 310
with cb_bootstrap.if_(self.step, "==", self.hist_length):
bootstrap_length = self.hist_length
with cb_bootstrap.if_(self.step, "==", bootstrap_length):
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, this was removed in 61277d8. An artifact of the old extend-bootstrap thingy that I overlooked.

for i in range(self.hist_length):
cb(time_data_var[i], time_data[i] - self.t)

relv_times = time_data_var
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
relv_times = time_data_var
relevant_times = time_data_var

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Addressed in 61277d8.

component_id=self.component_id,
time_id="", time=self.t)

def set_up_time_data(self, cb, new_t):
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Give me a docstring describing under which circumstances this is used, like pre/postcondition.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Added in 61277d8, but not sure if it was exactly what you're looking for. I went through the purpose of the function, described the argument, and what it takes in and returns. Happy to amend this docstring if it seems unsuitable.

dt_factor = self.dt
t_end = 1

return time_data, relv_times, dt_factor, t_end
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Explain difference between time_data and relevant_times.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I tried to add this to the docstring (added in 61277d8), but there are two main differences:

  • relevant_times subtracts out self.t, whereas time_data does not. I am truthfully not sure how critical this is.
  • more importantly, time_data is a list, whereas relevant_times is an array.

"than equations")
elif len(unknowns) < len(equations):
raise ValueError("Adams-Moulton implicit timestep has more equations "
"than unknowns")
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there an else to this last condition? If not:

Suggested change
"than unknowns")
"than unknowns")
else:
assert False

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Addressed in 61277d8.


# {{{ emit solve if possible

if unknowns and len(unknowns) == len(equations):
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Non-square systems aren't allowed IMO.

cb(my_rhs, last_rhss[name])
make_known(my_rhs)
def make_known(v):
unknowns.discard(v)
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is this here? Is there ever a scenario where a variable would be added to unknowns only to have to be removed later?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is something that I pulled in from the existing Butcher generator. Truthfully I'm not sure that such a scenario exists, though based on where this function is called, it looks like it might be possible in the case of stage recycling. Let me run some tests to determine whether this is even needed, because it's a fair question.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Addressed in 61277d8. I've removed this function and replaced calls to it with adding rhs's to knowns, but I agree with you where we should never need to remove things from unknowns, and getting rid of that has no ill effect that I can see.

"unknowns than equations")
elif len(unknowns) < len(equations):
raise ValueError("Runge-Kutta implicit timestep has more "
"equations than unknowns")
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

else assert False?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Addressed in 61277d8.

Comment on lines +389 to +391
def generate_butcher_finish(self, cb, stage_coeff_set_names,
stage_coeff_sets, rhs_funcs, estimate_coeff_set_names,
estimate_coeff_sets, stage_rhs_vars, estimate_vars):
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Explain the circumstances in which "finish" is invoked.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The reason that "finish" exists separately from "primary" essentially relates to the use of these functions in bootstrapping the multistep methods. This "finish" is invoked immediately after "primary" by regularly generated RK methods, because one of its primary responsibilities is yielding state and incrementing the time, but when bootstrapping a multistep method, this "finish" is not called, because other pieces of code are responsible for doing those things (not to mention there are instances in which other things need to occur in the "finish" of a bootstrap step for a multistep method pertaining to history updates).

cmikida2 and others added 2 commits December 4, 2020 13:51
Co-authored-by: Andreas Klöckner <[email protected]>
relv_times->relevant_times, else assert falses on solves
Base automatically changed from master to main March 8, 2021 05:09
@cmikida2
Copy link
Contributor Author

@inducer I've just looked this over again and brought it up to speed with the main leap branch. It remains, as far as I can see, ready for a further look, or at the very least is in a good position for someone else to pick it up and make any further tweaks necessary.

I'll also note that other PRs for adding adaptive implicit-explicit Adams schemes (#11 and #13) depend upon this work, at least in part.

@cmikida2 cmikida2 mentioned this pull request Aug 31, 2021
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