-
Notifications
You must be signed in to change notification settings - Fork 12
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
Conversation
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) |
There was a problem hiding this comment.
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?
There was a problem hiding this comment.
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 ' |
There was a problem hiding this comment.
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?
There was a problem hiding this comment.
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): |
There was a problem hiding this comment.
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
There was a problem hiding this comment.
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 |
There was a problem hiding this comment.
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
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thanks Tim!!
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 torho
,water_vapour
, andcloud_water
, and no transport foru
andtheta
. A new intermediate field ofxsimult
is used to store the output of the simultaneous transport. We then extract the fields of interest and assign them toxp
. This avoids overwriting any previously transported variables (for example, transportingu
ortheta
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 .