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

Non-split physics for all explicit or IMEX time discretisations #578

Draft
wants to merge 10 commits into
base: main
Choose a base branch
from

Conversation

atb1995
Copy link
Collaborator

@atb1995 atb1995 commented Nov 26, 2024

Currently you can only have physics when using IMEX RK, Explicit RK (Predictor) or IMEX SDC by using SpitPhysicsTimestepper. I will add the capability to do this, and a unit test to show it has the same results as SplitPhysicsTimestepper.

@atb1995 atb1995 self-assigned this Nov 26, 2024
@atb1995 atb1995 marked this pull request as draft November 26, 2024 14:50
@tommbendall tommbendall added the enhancement Pull requests or issues relating to adding a new capability label Dec 2, 2024
Copy link
Contributor

@tommbendall tommbendall left a comment

Choose a reason for hiding this comment

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

Thanks for this Alex. I've got two main thoughts:

  • is the test demonstrating that the physics is correct?
  • the explicit 'predictor' physics isn't what I was expecting so needs checking and testing

gusto/time_discretisation/imex_runge_kutta.py Outdated Show resolved Hide resolved
raise NotImplementedError(
'Physics not implemented with RK schemes that use the '
+ 'predictor form')
evaluate(self.field_i[stage], self.dt)
Copy link
Contributor

Choose a reason for hiding this comment

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

Are we sure this is correct? It doesn't match the comment -- I thought we needed to sum over the source terms and this would be difficult to do

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I think this is correct, unless why understanding of evaluate is wrong.. We evaluate the physics terms in field_i from the previous stage, since field_i is a list we have all previous states for each stage. The physics terms (non-time derivative) are then in the rhs that goes into the solver.

Copy link
Contributor

Choose a reason for hiding this comment

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

So the equation we are trying to solve is:

y^m = y^0 - dt*Sum_{i=0}^{m-1} a_{m,i} * F[y^i], for m = 1 to M - 1

What the evaluate method does is perform F[y] and put the value of that into a source field S. The thing that we don't have is M different source fields -- we only have one. So rather than having M source fields containing F[y^i] (for i=1 to m-1), and doing m-1 evaluations, we are overwriting S with ``F[y^{m-1}]` and not bringing in all of the correct contributions.

I don't think coding this is straightforward which is why I think it hasn't been done yet! But hopefully a test could demonstrate what is correct or not...

return error


@pytest.mark.parametrize("timestepper", ["split", "nonsplit_imex_rk", "nonsplit_imex_sdc",
Copy link
Contributor

Choose a reason for hiding this comment

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

Could we also add the explicit with the 'predictor' method?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Done


stepper.run(0, tmax=tmax)

error = norm(stepper.fields('f') - f_end) / norm(f_end)
Copy link
Contributor

Choose a reason for hiding this comment

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

I don't know if this test is too slack? I think we are mainly testing the transport her and less testing whether the physics is correct.

Could we:
(a) check against the result from a split-physics time stepper (which would have the downside of being inefficient unless we make the split-physics value be a KGO)
or (b) find a coupled physics test with a known solution?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement Pull requests or issues relating to adding a new capability
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants