Add background flow functionality in SingleLayerQG model

[mncrowe]

Hi,

I'm looking to add a constant background flow into the single layer model in order to study steadily propagating features. This is equivalent to adding im * params.C * grid.kr to the linear operator L for speed C in the x-direction. My options seem to be:

1. Use the multi-layer model with a constant U and 1 layer, I assume this would work though I'd prefer to have the translation term (and the β term) dealt with in the linear term for stability etc.
2. Create my own module based off SingleLayerQG. This seems straightforward enough as it's just adding C to the various functions and structures and adding the extra term into Equation. However, I'm not familiar with creating modules so am not sure how to load this without modifying GeophyiscalFlows.jl directly.
3. Ask you very nicely if you'll add in this feature, either via a new module based off mine, or by modifying SingleLayerQG to include this functionality.

What would you suggest?

Thanks,

Matt

[navidcy]

How about:

4. You add this feature in SingleLayerQG module through a PR? We can help you with that!

We can discuss here how you'd go about and do it if you like.

[mncrowe]

How about:

4. You add this feature in SingleLayerQG module through a PR? We can help you with that!

We can discuss here how you'd go about and do it if you like.

This sounds ideal, yes please! I have no experience contributing to projects that aren't mine so any help would be greatly appreciated.

[navidcy]

Sure. What you would like help with?

Github-wise: You first you need to fork the repo, create a branch, modify the source code (as if it was your code!) and then submit a pull request...

Or you want to discuss what/where exactly to modify the module?

[mncrowe]

Sure. What you would like help with?

Github-wise: You first you need to fork the repo, create a branch, modify the source code (as if it was your code!) and then submit a pull request...

Or you want to discuss what/where exactly to modify the module?

Thanks, this pretty much answered my question; which was 'do I need to create my own fork to do pull requests'.

I've now created the fork, and unless I've missed something, I think the required changes have been made. Will test later today and submit a pull request if it works.

[mncrowe]

changes added and tested and pull request submitted

[navidcy]

Questions:

1. So you only meant to add a constant U flow? I thought you were going for a U(y)... If only a constant U then how does that make any difference -- isn't it just a Galilean transformation?

2. What do you mean "and tested"? I don't see any tests in the PR. But let's clarify Q1 and then we can discuss about tests in Add background zonal flow U or U(y) in the SingleLayerQG module #360.

[mncrowe]

Yes, my goal was to add a constant flow U, I figured U(y) was already implemented in the multi-layer module which could presumably be run with 1 layer so this was unnecessary. Specifically, I wanted to add the term - im * params.U * grid.kr to the linear operator. Conversely U(y) would have to be added as a nonlinear term which would presumably be slower, and potentially less accurate/stable depending on the timestepper.

You're correct that this is just a Galilean transformation. However, my goal is to simulate travelling vortices (such as the Lamb vortex you have implemented in lambdipole) with sponge layers to damp any waves or vorticity filaments. The options are to either move the sponge layers at the same speed as the vortex which requires updating the damping term at each timestep or to work in a frame where the vortex is stationary and the sponge layer is fixed in time. The second option is much quicker here as adding the translation into the linear operator doesn't slow down the simulation and the sponge layer is fixed, reducing the number of function evaluations.

By 'tested' I mean that I have created a new julia enviroment and added my forked version of GeophyiscalFlows.jl to it. All example scripts from before still run as expected and setting the parameter U gives the expected behaviour. Not sure what tests you can/should do in the PR or how to do them.

[glwagner]

Would U(y) need to be periodic for this to work? (Or perhaps dU/dy?)

One test is to make sure the code runs without erroring.

Another is to test correctness, ie use some metric + simulation to show that U causes things to propagate at the right speed.

If there's no test it could easily break without anyone knowing, so it's in your interest to add a test to make sure your code continues to work in the future

[navidcy]

@mncrowe thanks for the elaboration.

Fair enough for the constant U.

However, regarding testing we would need something more than taking your word for it. I mean, not that we don't believe you, but it's good if we have something in the test suite of the repo that uses the functionality you added and tests that the results as what being expected. This way, we would have a way to catch a potential bug introduced in the future.

E.g., for the SingleLayerQG module there is a tests that evolves a Rossby wave for a bit and checks whether the wave moves with the speed predicted by the dispersion relationship. It seems easy to have another case with a non-zero U and check again that the wave moves with what is expected then? Something like that.

[glwagner]

@navidcy I would just add that the test is actually most important to make sure the code keeps working in the future. I think if we manually verify the code works now, I think we can conclude that it does indeed work. The issue isn't figuring out whether the code works now. The thing we are actually trying to achieve with CI is to ensure that the code continues to work in the future as new PRs come in.

[apaloczy]

Just to add another thought about this, would it make sense to have a test comparing results between SingleLayerQG and MultiLayerQG (with nlayers = 1) simulations with a constant U? I don't know what kind of bug that could help catch, but it's a sanity check.

Yes, my goal was to add a constant flow U, I figured U(y) was already implemented in the multi-layer module which could presumably be run with 1 layer so this was unnecessary. Specifically, I wanted to add the term - im * params.U * grid.kr to the linear operator. Conversely U(y) would have to be added as a nonlinear term which would presumably be slower, and potentially less accurate/stable depending on the timestepper.

And about this, I wonder if it would be helpful to use multiple dispatch to put the background flow in the linear operator if params.U is a float, and put it in the nonlinear operator (as done the MultiLayerQG module) if it is a 1D Array?

[glwagner]

And about this, I wonder if it would be helpful to use multiple dispatch to put the background flow in the linear operator if params.U is a float, and put it in the nonlinear operator (as done the MultiLayerQG module) if it is a 1D Array?

Yeah for sure or just an if statement, presumably that'd be easy

[mncrowe]

And about this, I wonder if it would be helpful to use multiple dispatch to put the background flow in the linear operator if params.U is a float, and put it in the nonlinear operator (as done the MultiLayerQG module) if it is a 1D Array?

Yeah for sure or just an if statement, presumably that'd be easy

Agreed, will have a look at adding non-constant U to the nonlinear operator.

Regarding tests; how do I add them here and how can I locate the pre-existing ones to ensure I use a consistent approach?

[glwagner]

The tests are all in this file: https://github.com/FourierFlows/GeophysicalFlows.jl/blob/main/test/test_singlelayerqg.jl

[mncrowe]

And about this, I wonder if it would be helpful to use multiple dispatch to put the background flow in the linear operator if params.U is a float, and put it in the nonlinear operator (as done the MultiLayerQG module) if it is a 1D Array?

Yeah for sure or just an if statement, presumably that'd be easy

For consistency with the rest of the package, how would you prefer this if done? e.g. Comparing with typeof or checking size of array etc.

Also, would you prefer multiple Params structs for U as Float vs Array{Float}?

[glwagner]

And about this, I wonder if it would be helpful to use multiple dispatch to put the background flow in the linear operator if params.U is a float, and put it in the nonlinear operator (as done the MultiLayerQG module) if it is a 1D Array?

Yeah for sure or just an if statement, presumably that'd be easy

For consistency with the rest of the package, how would you prefer this if done? e.g. Comparing with typeof or checking size of array etc.

Also, would you prefer multiple Params structs for U as Float vs Array{Float}?

if U isa Number
   # add to linear operator
else
    # add to nonlinear RHS
end