Horizontal & vertical splitting of advective transport terms, and a DG_Upwind version #616
Conversation
atb1995
added
enhancement
gusto/core/labels.py
Outdated
| @@ -92,6 +92,8 @@ def __call__(self, target, value=None): | |||
| # ---------------------------------------------------------------------------- # | |||
| implicit = Label("implicit") | |||
| explicit = Label("explicit") | |||
| horizontal = Label("horizontal") | |||
There was a problem hiding this comment.
I think we might want to slightly tweak the names of these labels. In #584 I would introduce labels for the vertical and horizontal components of the transported wind, so we don't want to confuse things.
What would you think to horizontal_transport and vertical_transport? I think these labels will only apply to transport, so is that reasonable?
| @@ -44,7 +44,7 @@ def __init__(self, equation, variable, term_label): | |||
| map_if_true=keep, map_if_false=drop) | |||
|
|
|||
| num_terms = len(self.original_form.terms) | |||
| assert num_terms == 1, f'Unable to find {term_label.label} term ' \ | |||
| assert num_terms >= 1, f'Unable to find {term_label.label} term ' \ | |||
There was a problem hiding this comment.
I'm nervous about this change! Is this because we are replacing a single 3D term with horizontal and vertical discretisation terms? I'd be worried that we could unintentionally cause multiple terms to be replaced when we don't mean to -- is there a way we can make this check safer?
There was a problem hiding this comment.
I've changed this to loop through a list of term_labels, but I'm not sure if this is the best solution!
There was a problem hiding this comment.
I've actually changed it to pass through a list of term labels, then check if number of terms per label is correct!
| u""" | ||
| Splits advective term into horizontal and vertical terms. | ||
| This describes splitting u.∇(q) terms into u.(∇_h)q and w dq/dz, | ||
| for transporting velocity u and transported q. |
There was a problem hiding this comment.
I think this all looks correct. It's very cumbersome having to go through the whole of this process -- do you think there are any shortcuts we can take, e.g. in adding the linearisations?
There was a problem hiding this comment.
I've tried to tidy it a bit by splitting out the advection forms, hopefully this helps!
| @@ -111,8 +111,41 @@ def replace_form(self, equation): | |||
| map_if_true=lambda t: new_term) | |||
|
|
|||
| else: | |||
| raise RuntimeError('Found multiple transport terms for ' | |||
There was a problem hiding this comment.
Again I'm a bit worried about just removing this! Can we have a special check for the horizontal/vertical case?
|
|
||
| # Check if this is a conservative transport | ||
| if horizontal_term.has_label(mass_weighted) or vertical_term.has_label(mass_weighted): | ||
| raise RuntimeError('Mass weighted transport terms not yet supported for multiple terms') |
There was a problem hiding this comment.
Should this be a NotImplementedError? Do you think it would be worth raising an issue to capture this debt once implemented?
There was a problem hiding this comment.
Done and issue raised #627
| @@ -277,9 +310,133 @@ def __init__(self, equation, variable, ibp=IntegrateByParts.ONCE, | |||
| self.form = form | |||
|
|
|||
|
|
|||
| class Split_DGUpwind(TransportMethod): | |||
There was a problem hiding this comment.
| class Split_DGUpwind(TransportMethod): | |
| class SplitDGUpwind(TransportMethod): |
| @@ -15,19 +15,20 @@ class SpatialMethod(object): | |||
| The base object for describing a spatial discretisation of some term. | |||
| """ | |||
|
|
|||
| def __init__(self, equation, variable, term_label): | |||
| def __init__(self, equation, variable, term_labels): | |||
There was a problem hiding this comment.
We discussed this offline but putting it here for posterity: could we avoid some of the complicated logic below by using *term_labels here, or similar?
Can we still enforce that the original labelled form still only has one valid term to be replaced?
| return norm(timestepper.fields("f") - f_end) / norm(f_end) | ||
|
|
||
|
|
||
| @pytest.mark.parametrize("geometry", ["slice", "sphere"]) |
There was a problem hiding this comment.
Actually this sphere is only 2D so not really testing the splitting, so shall we remove this?
No description provided.