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

Implement simultaneous transport with SIQN #581

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

Conversation

ta440
Copy link
Collaborator

@ta440 ta440 commented Nov 28, 2024

Allow conservative transport with the SIQN timestepper. To do this, the variables involved in the conserved transport (mixing ratio(s) and the reference density) are passed as a list into a single time discretisation instance, e.g. for the moist, compressible Euler equations:
SSPRK3(domain, ["rho", "water_vapour", "cloud_water"], options=mixed_opts, rk_formulation=RungeKuttaFormulation.predictor)
If there is a mixed wrapper, any boundary conditions are applied to the new mixed function space.
In the transport step, the entire mixed function space is evolved, with any variables not in the list remaining untransported. In this example, transport is applied to the u_rho_theta_water_vapour_cloud_water field, with transport schemes applied to rho, water_vapour, and cloud_water, and no transport for u and theta. A new intermediate field of xsimult is used to store the output of the simultaneous transport. We then extract the fields of interest and assign them to xp. This avoids overwriting any previously transported variables (for example, transporting u or theta first).
This pull request also fixes the implementation of boundary conditions when using MixedFSOptions. This means that MixedFSOptions can now be extended to the prognostic equations sets of Shallow Water, Boussinesq, Compressible Euler, when it was previously breaking; this addresses issue #549 .

@ta440 ta440 added the enhancement Pull requests or issues relating to adding a new capability label Dec 2, 2024
@ta440 ta440 marked this pull request as ready for review December 2, 2024 20:19
for i in range(len(active_tracers) - 1):
tracer = active_tracers[i]
if tracer.transport_eqn == TransportEquationType.tracer_conservative:
ref_density = next(x for x in active_tracers if x.name == tracer.density_name)
ref_density = next((x for x in active_tracers if x.name == tracer.density_name), tracer)
Copy link
Contributor

Choose a reason for hiding this comment

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

Do we want tracer to be a default value, or would it be better to fail if we can't find the density name?

Copy link
Collaborator Author

@ta440 ta440 Dec 3, 2024

Choose a reason for hiding this comment

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

There is a reason for this default value. We might need to reorder the active tracers if the density is ordered after a conservatively transported tracer. However, this previously failed when rho is defined as a field not an activate tracer, like with the CompressibleEulerEquations, as next() won't find any match within the active tracer list. Hence, I set the default to be tracer, so that i=j and it won't try to reorder the active tracers. I'll add a comment to this effect.

@@ -175,7 +207,9 @@ def setup(self, equation, apply_bcs=True, *active_labels):
# timestepper should be used instead.
if len(field_terms.label_map(lambda t: t.has_label(mass_weighted), map_if_false=drop)) > 0:
if len(field_terms.label_map(lambda t: not t.has_label(mass_weighted), map_if_false=drop)) > 0:
raise ValueError(f"Mass-weighted and non-mass-weighted terms are present in a timestepping equation for {field}. As these terms cannot be solved for simultaneously, a split timestepping method should be used instead.")
raise ValueError('Mass-weighted and non-mass-weighted terms are present in a timestepping '
Copy link
Contributor

Choose a reason for hiding this comment

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

Can we alter these comments to change how they are split over multiple lines?

Copy link
Contributor

Choose a reason for hiding this comment

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

Ah I realise you're already reducing the length of the lines here -- but can we go even further so that each line is 80 characters or less? If it looks messier then please ignore me!

logger.info(f'Semi-implicit Quasi Newton: Transport {outer}: {name}')
# transports a field from xstar and puts result in xp
self.transport_field(name, scheme, xstar, xp)
if isinstance(name, list):
Copy link
Contributor

Choose a reason for hiding this comment

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

Can we find a way to move this if statement and the comments into the transport_field routine? I would quite like to keep this part of the timeloop looking as simple as:
self.transport_field(...)

For instance maybe we should move the iteration into transport_field so this call becomes:

self.transport_field(xstar, xp)

and self.transport_field is defined so that:

def transport_field(self, xstar, xp):
    for name, scheme in self.active_transport:
        ... # do everything we previously did

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 agree that this is a nice way of doing it. I've implemented this change, which assumes that we won't want to use a divergence predictor with simultaneous transport.


def run_simult_SIQN(tmpdir, order):

ncolumns = 50
Copy link
Contributor

Choose a reason for hiding this comment

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

Can we reduce the resolution, since this is just for testing. I think 10 columns/layers should be enough for the higher-order case, and 20 for the lower-order case

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