Cory implicit adams #4

cmikida2 · 2020-09-10T19:33:54Z

Adding implicit to Adams (Adams-Moulton), and incorporating the ability to treat certain RHS components implicitly in a multi-rate setting.

(Continuation of https://gitlab.tiker.net/inducer/leap/-/merge_requests/38)

Backward Euler tableau

avoid extra RHS evals - also uses implicit DIRK for bootstrap

IMEX and fully explicit

test to pass

tolerance

inducer · 2020-09-10T19:47:25Z

Thanks! A few checks failed. Could you take a look?

cmikida2 · 2020-09-10T19:49:02Z

Yeah, I'm already on it. It's small stuff, should be easily fixed soon. Will comment again when done.

cmikida2 · 2020-09-10T20:17:31Z

All checks are now passing, but a caveat (I'd say a major one) is that none of the implicit tests are actually being run, as they require scipy to run, and that was causing a large number of failures initially because the CI doesn't have scipy (I think that's what the issue is?). The other implicit tests that already existed prior to me (IMEX stuff) have importorskip statements that circumvent this.

Is there a way to add scipy to the CI so the implicit tests are run? Further, if there's a way for me to do this/if you want me to do it, let me know and I'm happy to handle it.

inducer

A quick pre-review since I'd like the implicit CI to actually run.

leap/rk/imex.py

test/test_multirate.py

cmikida2 · 2020-09-10T20:37:21Z

Sorry about that - CI is now actually running (and passing) the implicit tests with scipy.

inducer · 2020-10-24T18:14:50Z

Btw, I'm in the middle of reviewing this. If/when you go to add adaptivity, could you do that on a new branch off of this one?

inducer

Thanks for your work on this. Some comments below. I propose that we take this one step at a time. I've only reviewed the AM bits for now. Could you split those into a separate PR? Big PRs like this one are difficult and time-consuming to review.

inducer · 2020-09-18T01:10:34Z

leap/multistep/__init__.py

+    User-supplied context:
+        <state> + component_id: The value that is integrated
+        <func> + component_id: The right hand side


I suspect this won't render the way you intend. Check via

pip install sphinx cd doc make html open _build/html/index.html

inducer · 2020-10-22T00:28:09Z

leap/multistep/__init__.py

+            # but if we are comparing with IMEX MRAM, we need one more
+            # bootstrap step.


I'm not sure I fully understand this comment. What's being compared here?

inducer · 2020-10-22T00:38:26Z

leap/multistep/__init__.py

+    def rk_bootstrap(self, cb):
+        """Initialize the timestepper with an IMPLICIT RK method."""
+
+        equations = []
+        unknowns = set()
+        knowns = set()
+        rhs_var_to_unknown = {}
+
+        from leap.rk import IMPLICIT_ORDER_TO_RK_METHOD_BUILDER
+        rk_method = IMPLICIT_ORDER_TO_RK_METHOD_BUILDER[self.function_family.order]
+        rk_tableau = tuple(zip(rk_method.c, rk_method.a_implicit))
+        rk_coeffs = rk_method.output_coeffs
+
+        if self.extra_bootstrap:
+            first_save_step = 2
+        else:
+            first_save_step = 1
+
+        with cb.if_(self.step, "==", first_save_step):
+            # Save the first RHS to the AM history
+            rhs_var = var("rhs_var")
+
+            cb(rhs_var, self.eval_rhs(self.t, self.state))
+            cb(self.history[0], rhs_var)
+
+            if not self.static_dt:
+                cb(self.time_history[0], self.t)
+
+        # Stage loop
+        rhss = [var("rk_rhs_" + str(i)) for i in range(len(rk_tableau))]
+        for stage_num, (c, coeffs) in enumerate(rk_tableau):
+            stage = self.state + sum(self.dt * coeff * rhss[j]
+                                     for (j, coeff)
+                                     in enumerate(coeffs))
+
+            if self.state_filter is not None:
+                stage = self.state_filter(stage)
+
+            # In a DIRK setting, the unknown is always the same RHS
+            # as the stage number.
+            unknowns.add(rhss[stage_num])
+            unkvar = cb.fresh_var('unk_s%d' % (stage_num))
+            rhs_var_to_unknown[rhss[stage_num]] = unkvar
+            from dagrt.expression import collapse_constants
+            solve_expression = collapse_constants(
+                    rhss[stage_num] - self.eval_rhs(self.t + c*self.dt, stage),
+                    list(unknowns) + [self.state],
+                    cb.assign, cb.fresh_var)
+            equations.append(solve_expression)
+
+            # {{{ emit solve if possible
+
+            if unknowns and len(unknowns) == len(equations):
+                # got a square system, let's solve
+                assignees = [unk.name for unk in unknowns]
+
+                from pymbolic import substitute
+                subst_dict = dict(
+                        (rhs_var.name, rhs_var_to_unknown[rhs_var])
+                        for rhs_var in unknowns)
+
+                cb.assign_implicit(
+                        assignees=assignees,
+                        solve_components=[
+                            rhs_var_to_unknown[unk].name
+                            for unk in unknowns],
+                        expressions=[
+                            substitute(eq, subst_dict)
+                            for eq in equations],
+
+                        # TODO: Could supply a starting guess
+                        other_params={
+                            "guess": self.state},
+                        solver_id="solve")
+
+                del equations[:]
+                knowns.update(unknowns)
+                unknowns.clear()
+
+        # }}}
+
+        # Merge the values of the RHSs.
+        rk_comb = sum(coeff * rhss[j] for j, coeff in enumerate(rk_coeffs))
+
+        state_est = self.state + self.dt * rk_comb
+        if self.state_filter is not None:
+            state_est = self.state_filter(state_est)
+
+        # Assign the value of the new state.
+        cb(self.state, state_est)
+
+        # Save the "next" RHS to the AM history
+        rhs_next_var = var("rhs_next_var")
+
+        cb(rhs_next_var, self.eval_rhs(self.t + self.dt, self.state))
+
+        for i in range(1, len(self.history)):
+            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)


I expect that this is, at least to an extent, a version of the RK integration code, with some number of changes. Could you refactor this to use (an improved version of) generate_butcher to share code (or explain why that's not reasonable)?

inducer · 2020-10-22T00:38:48Z

leap/multistep/__init__.py

+                    },
+                initial_phase="initial")
+
+        # Bootstrap


Suggested change

# Bootstrap

# {{{ bootstrap

inducer · 2020-10-22T00:39:11Z

leap/multistep/__init__.py

+                    cb_bootstrap.switch_phase("primary")
+            else:
+                with cb_bootstrap.if_(self.step, "==", self.hist_length - 1):
+                    cb_bootstrap.switch_phase("primary")


Suggested change

cb_bootstrap.switch_phase("primary")

cb_bootstrap.switch_phase("primary")

# }}}

inducer · 2020-10-26T01:40:02Z

leap/multistep/__init__.py

+# {{{ am method
+


Suggested change

# {{{ am method

# {{{ am method

inducer · 2020-10-26T01:42:05Z

leap/multistep/__init__.py

+        if state_filter_name is not None:
+            self.state_filter = var("<func>" + state_filter_name)
+        else:
+            self.state_filter = None


It seems this constructor is largely the same as the AB one. Why is its code duplicated? (e.g. instead of inheriting it from a common base)

inducer · 2020-10-26T01:44:48Z

leap/multistep/__init__.py

+            # }}}
+
+            # Update the state now that we've solved.
+            cb_primary(self.state, state_est)


How come the solve result is not passed through the state filter?

inducer · 2020-10-26T01:46:43Z

leap/multistep/__init__.py

+            # Rotate history and time history.
+            for i in range(self.hist_length - 1):
+                cb_primary(self.history[i], history[i + 1])
+
+                if not self.static_dt:
+                    cb_primary(self.time_history[i], time_history_data[i + 1])


I imagine this time history management code occurs similarly somewhere else. Could you factor this out into a utility?

inducer · 2020-10-26T01:51:12Z

test/test_implicit_ab.py

+    print(eocrec.pretty_print())
+
+    orderest = eocrec.estimate_order_of_convergence()[0, 1]
+    assert orderest > expected_order * 0.9


It looks like this code largely duplicates what's already in test_imex.py. I would prefer if this could be expressed in one shared test.

cmikida2 · 2020-10-26T02:46:16Z

Hi - yes, I can split the AM only parts into a separate PR for easier review. I'll also address the comments above in that PR, try to slim down the code/combine things where I can, and ping you again with a fresh PR when that's done. As for adaptivity, I had already started that work by branching from this branch, as you suggested above. For now, I'll obviously prioritize cleaning this up more before continuing with adaptivity.

cmikida2 · 2020-10-27T20:32:32Z

Hi - I've moved the single-rate implicit Adams only to a new pull request here.. I'm not sure how you want to manage this - I believe I've responded to all of the comments you've put above, but I'm not sure what protocol is/whether I should answer the above comments with references to the code in that PR, or if we just start fresh there. Happy to do whatever, and thank you again/in advance for your continued feedback.

cmikida2 added 16 commits June 18, 2020 17:20

Initial stab at implicit Adams method generation

cb564e5

Adds implicit DIRK schemes for Adams-Moulton bootstrapping, fixes

e465177

Backward Euler tableau

Fixes single-rate AM to use history in a way more consistent with AB and

44baab1

avoid extra RHS evals - also uses implicit DIRK for bootstrap

Fix copyright line, add first order to tests

1d1e717

Simplifies solver hooks by applying distribution to solve expressions

1905ae4

Sets up IMEX bootstrapping for IMEX MR

f9f6619

First pass at implicit-explicit Adams-based multi-rate

62de9f8

Make minimal changes required to combine bootstrap step functions for

c1b563c

IMEX and fully explicit

Minor copyright/docstring updates

47d9b73

Placate flake8 by removing needless state update

4b8c9bf

Add extra bootstrap option to AM to get implicit single-rate identical

dd7990a

test to pass

Add extra bootstrap option to single-rate identical test, tweak

6945909

tolerance

PyLint fix (whoopsie)

b066fbc

More PyLint fixes

266d87d

Merge branch 'cory-adams-rename' into cory-implicit-adams

e854c74

Finish fixing Adams/AB renaming issues

83c7530

cmikida2 added 2 commits September 10, 2020 14:50

Fix errant class reference in documentation

cc24a36

Importorskip(scipy) for implicit tests - CI will now skip these

112b32b

inducer reviewed Sep 10, 2020

View reviewed changes

leap/rk/imex.py Show resolved Hide resolved

test/test_multirate.py Outdated Show resolved Hide resolved

test/test_multirate.py Outdated Show resolved Hide resolved

cmikida2 added 2 commits September 10, 2020 15:22

Add scipy to ci.yml so implicit tests run

75d3b7e

Remove importor skips so AM tests run in CI

571640e

cmikida2 added 3 commits October 23, 2020 14:30

Merge branch 'master' into cory-implicit-adams

aa19b84

Fix flake8 quotes

4d41c8e

Missed a bad quote

6d4d9ad

inducer reviewed Oct 26, 2020

View reviewed changes

cmikida2 mentioned this pull request Oct 27, 2020

Cory Adams Moulton (V2) #10

Open

Base automatically changed from master to main March 8, 2021 05:09

		# but if we are comparing with IMEX MRAM, we need one more
		# bootstrap step.

Cory implicit adams #4

Are you sure you want to change the base?

Cory implicit adams #4

Uh oh!

Conversation

cmikida2 commented Sep 10, 2020 • edited by inducer Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

inducer commented Sep 10, 2020

Uh oh!

cmikida2 commented Sep 10, 2020

Uh oh!

cmikida2 commented Sep 10, 2020

Uh oh!

inducer left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

cmikida2 commented Sep 10, 2020

Uh oh!

inducer commented Oct 24, 2020

Uh oh!

inducer left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

cmikida2 commented Oct 26, 2020

Uh oh!

cmikida2 commented Oct 27, 2020

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

cmikida2 commented Sep 10, 2020 •

edited by inducer

Loading