From 8e14e7a00e17ce366d8ec65f17e174e7a3f6198b Mon Sep 17 00:00:00 2001 From: Tim Andrews Date: Fri, 10 Nov 2023 23:48:12 +0000 Subject: [PATCH 01/10] start MixedFSLimiter --- gusto/limiters.py | 52 +++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 50 insertions(+), 2 deletions(-) diff --git a/gusto/limiters.py b/gusto/limiters.py index e30a1200c..fc7ec04a5 100644 --- a/gusto/limiters.py +++ b/gusto/limiters.py @@ -6,13 +6,13 @@ """ from firedrake import (BrokenElement, Function, FunctionSpace, interval, - FiniteElement, TensorProductElement) + FiniteElement, TensorProductElement, split) from firedrake.slope_limiter.vertex_based_limiter import VertexBasedLimiter from gusto.kernels import LimitMidpoints import numpy as np -__all__ = ["DG1Limiter", "ThetaLimiter", "NoLimiter"] +__all__ = ["DG1Limiter", "ThetaLimiter", "NoLimiter", "MixedFSLimiter"] class DG1Limiter(object): @@ -178,3 +178,51 @@ def apply(self, field): applied, if this was not a blank limiter. """ pass + + +class MixedFSLimiter(object): + """A dictionary that holds limiters for individual components. + If no limiter is specified, then we apply a blank limiter""" + + + def __init__(self, tracers, sublimiters): + self.tracers = tracers + self.sublimiters = sublimiters + + for tracer in tracers: + #print(tracer.name) + #sublimiter = sublimiters.get(tracer.name) + #print(sublimiter) + if tracer.name not in sublimiters: + sublimiters[tracer.name] = NoLimiter() + sublimiter = sublimiters.get(tracer.name) + print(sublimiter) + + self.sublimiters = sublimiters + print(len(self.tracers)) + print(len(self.sublimiters)) + + # Check that each tracer has a defined sublimiter + # Else, give it a blank limiter instead. + + + def apply(self, fields): + """ + Apply individual limiter + """ + + #Iterate over each tracer + #Apply the sublimiter + + for tracer in self.tracers: + sublimiter = self.sublimiters[tracer.name] + print(sublimiter) + + # split specific tracer from x_in + # then apply the limiter to this. + + sublimiter.apply(tracer.field) + + + + \ No newline at end of file From 6062ebfde05e087adc013368e66fa53ce5cb2ca3 Mon Sep 17 00:00:00 2001 From: Tim Andrews Date: Tue, 14 Nov 2023 02:14:52 +0000 Subject: [PATCH 02/10] MixedFSLimiter is operational --- gusto/limiters.py | 75 +++++++++++++++++++++++------------------------ 1 file changed, 37 insertions(+), 38 deletions(-) diff --git a/gusto/limiters.py b/gusto/limiters.py index fc7ec04a5..b8d97382f 100644 --- a/gusto/limiters.py +++ b/gusto/limiters.py @@ -181,48 +181,47 @@ def apply(self, field): class MixedFSLimiter(object): - """A dictionary that holds limiters for individual components. - If no limiter is specified, then we apply a blank limiter""" - + """An object to hold a dictionary that defines limiters for transported prognostic + variables. Different limiters may be applied to different fields and not every transported variable needs a defined limiter. + """ - def __init__(self, tracers, sublimiters): - self.tracers = tracers - self.sublimiters = sublimiters - - for tracer in tracers: - #print(tracer.name) - #sublimiter = sublimiters.get(tracer.name) - #print(sublimiter) - if tracer.name not in sublimiters: - sublimiters[tracer.name] = NoLimiter() - sublimiter = sublimiters.get(tracer.name) - print(sublimiter) - + def __init__(self, equation, sublimiters): + """ + Args: + equation (:class: 'PrognosticEquationSet'): the prognostic equation(s) + sublimiters (dict): A dictionary holding limiters defined for individual prognostic variables + Raises: + ValueError: If a limiter is defined for a field that is not in the prognostic variable set + """ self.sublimiters = sublimiters - print(len(self.tracers)) - print(len(self.sublimiters)) - - # Check that each tracer has a defined sublimiter - # Else, give it a blank limiter instead. - + for tracer, sublimiter in sublimiters.items(): + # Check that the tracer is being solved in the equation: + if tracer not in equation.field_names: + raise ValueError(f"The limiter for {tracer} is for a field that does not exist in the equation set") + else: + self.sublimiters[tracer].idx = equation.field_names.index(tracer) + + #for tracer in tracers: + # if tracer.name in sublimiters: + # sublimiter = sublimiters.get(tracer.name) + # print(sublimiter) + # self.sublimiters[tracer.name].idx = equation.field_names.index(tracer.name + def apply(self, fields): """ - Apply individual limiter + Apply the individual limiters to specific prognostic variables """ - #Iterate over each tracer - #Apply the sublimiter - - for tracer in self.tracers: - sublimiter = self.sublimiters[tracer.name] - print(sublimiter) - - # split specific tracer from x_in - # then apply the limiter to this. - - sublimiter.apply(tracer.field) - - - - \ No newline at end of file + #for tracer in self.tracers: + # if tracer.name in self.sublimiters: + # sublimiter = self.sublimiters[tracer.name] + # field = fields.subfunctions[sublimiter.idx] + # sublimiter.apply(field) + # print(f'Applying sublimiter on {tracer.name} field') + + for _, sublimiter in self.sublimiters.items(): + field = fields.subfunctions[sublimiter.idx] + sublimiter.apply(field) + #print(f'Applying sublimiter on {tracer} field') + \ No newline at end of file From 18867a5a9163b94692272d3609ff7ff6db740c68 Mon Sep 17 00:00:00 2001 From: Tim Andrews Date: Tue, 14 Nov 2023 23:38:22 +0000 Subject: [PATCH 03/10] added an integration test for a mixed function space limiter --- integration-tests/transport/test_limiters.py | 141 +++++++++++++++---- 1 file changed, 113 insertions(+), 28 deletions(-) diff --git a/integration-tests/transport/test_limiters.py b/integration-tests/transport/test_limiters.py index 8255ad2cd..0e507dcba 100644 --- a/integration-tests/transport/test_limiters.py +++ b/integration-tests/transport/test_limiters.py @@ -50,17 +50,36 @@ def setup_limiters(dirname, space): VDG1 = domain.spaces('DG1_equispaced') elif space == 'Vtheta_degree_1': V = domain.spaces('theta') + elif space == 'mixed_FS': + V = domain.spaces('HDiv') + VA = domain.spaces('DG') + #VB = domain.spaces('theta') + VB = domain.spaces('DG') else: raise NotImplementedError Vpsi = domain.spaces('H1') # Equation - eqn = AdvectionEquation(domain, V, 'tracer') - - # I/O - output = OutputParameters(dirname=dirname+'/limiters', + if space == 'mixed_FS': + tracerA = ActiveTracer(name='tracerA', space='DG', + variable_type=TracerVariableType.mixing_ratio, + transport_eqn=TransportEquationType.advective) + #tracerB = ActiveTracer(name='tracerB', space='theta', + # variable_type=TracerVariableType.mixing_ratio, + # transport_eqn=TransportEquationType.advective) + tracerB = ActiveTracer(name='tracerB', space='DG', + variable_type=TracerVariableType.mixing_ratio, + transport_eqn=TransportEquationType.advective) + tracers = [tracerA, tracerB] + eqn = CoupledTransportEquation(domain, active_tracers=tracers, Vu = V) + output = OutputParameters(dirname=dirname+'/limiters', + dumpfreq=1, dumplist=['u', 'tracerA', 'tracerB', 'true_tracerA', 'true_tracerB']) + else: + eqn = AdvectionEquation(domain, V, 'tracer') + output = OutputParameters(dirname=dirname+'/limiters', dumpfreq=1, dumplist=['u', 'tracer', 'true_tracer']) + io = IO(domain, output) # ------------------------------------------------------------------------ # @@ -84,10 +103,24 @@ def setup_limiters(dirname, space): elif space == 'Vtheta_degree_1': opts = EmbeddedDGOptions() transport_schemes = SSPRK3(domain, options=opts, limiter=ThetaLimiter(V)) + + elif space == 'mixed_FS': + sublimiters = {'tracerA': DG1Limiter(domain.spaces('DG')), + 'tracerB': DG1Limiter(domain.spaces('DG'))} + #'tracerB': ThetaLimiter(domain.spaces('theta'))} + # Need Embedded DG if wanting to test theta limiter. + MixedLimiter = MixedFSLimiter(eqn, sublimiters) + transport_schemes = SSPRK3(domain, limiter=MixedLimiter) + + #theta_opts = EmbeddedDGOptions() + else: raise NotImplementedError - transport_method = DGUpwind(eqn, "tracer") + if space == 'mixed_FS': + transport_method = [DGUpwind(eqn, 'tracerA'), DGUpwind(eqn, 'tracerB')] + else: + transport_method = DGUpwind(eqn, "tracer") # Build time stepper stepper = PrescribedTransport(eqn, transport_schemes, io, transport_method) @@ -96,8 +129,14 @@ def setup_limiters(dirname, space): # Initial condition # ------------------------------------------------------------------------ # - tracer0 = stepper.fields('tracer', V) - true_field = stepper.fields('true_tracer', space=V) + if space == 'mixed_FS': + tracerA_0 = stepper.fields('tracerA', space=VA) + tracerB_0 = stepper.fields('tracerB', space=VB) + true_fieldA = stepper.fields('true_tracerA', space=VA) + true_fieldB = stepper.fields('true_tracerB', space=VB) + else: + tracer0 = stepper.fields('tracer', V) + true_field = stepper.fields('true_tracer', space=V) x, z = SpatialCoordinate(mesh) @@ -146,9 +185,17 @@ def setup_limiters(dirname, space): 0.0) if i == 0: - tracer0.interpolate(Constant(tracer_min) + expr_1 + expr_2) + if space=='mixed_FS': + tracerA_0.interpolate(Constant(tracer_min) + expr_1 + expr_2) + tracerB_0.interpolate(Constant(tracer_min) + expr_1 + expr_2) + else: + tracer0.interpolate(Constant(tracer_min) + expr_1 + expr_2) elif i == 1: - true_field.interpolate(Constant(tracer_min) + expr_1 + expr_2) + if space=='mixed_FS': + true_fieldA.interpolate(Constant(tracer_min) + expr_1 + expr_2) + true_fieldB.interpolate(Constant(tracer_min) + expr_1 + expr_2) + else: + true_field.interpolate(Constant(tracer_min) + expr_1 + expr_2) else: raise ValueError @@ -181,29 +228,67 @@ def setup_limiters(dirname, space): gradperp = lambda v: as_vector([-v.dx(1), v.dx(0)]) u.project(gradperp(psi)) - return stepper, tmax, true_field - + if space == 'mixed_FS': + return stepper, tmax, true_fieldA, true_fieldB + else: + return stepper, tmax, true_field -@pytest.mark.parametrize('space', ['Vtheta_degree_0', 'Vtheta_degree_1', - 'DG0', 'DG1', 'DG1_equispaced']) +#move mixed_FS to the end once finished debugging +@pytest.mark.parametrize('space', ['Vtheta_degree_0', 'mixed_FS'])#, 'Vtheta_degree_1', + #'DG0', 'DG1', 'DG1_equispaced']) def test_limiters(tmpdir, space): # Setup and run dirname = str(tmpdir) - stepper, tmax, true_field = setup_limiters(dirname, space) + + if space == 'mixed_FS': + stepper, tmax, true_fieldA, true_fieldB = setup_limiters(dirname, space) + else: + stepper, tmax, true_field = setup_limiters(dirname, space) + stepper.run(t=0, tmax=tmax) - final_field = stepper.fields('tracer') - - # Check tracer is roughly in the correct place - assert norm(true_field - final_field) / norm(true_field) < 0.05, \ - 'Something appears to have gone wrong with transport of tracer using a limiter' - + tol = 1e-9 - - # Check for no new overshoots - assert np.max(final_field.dat.data) <= np.max(true_field.dat.data) + tol, \ - 'Application of limiter has not prevented overshoots' - - # Check for no new undershoots - assert np.min(final_field.dat.data) >= np.min(true_field.dat.data) - tol, \ - 'Application of limiter has not prevented undershoots' + + if space == 'mixed_FS': + final_fieldA = stepper.fields('tracerA') + final_fieldB = stepper.fields('tracerB') + + # Check tracer is roughly in the correct place + assert norm(true_fieldA - final_fieldA) / norm(true_fieldA) < 0.05, \ + 'Something is wrong with the DG space tracer using a mixed limiter' + + # Check tracer is roughly in the correct place + assert norm(true_fieldB - final_fieldB) / norm(true_fieldB) < 0.05, \ + 'Something is wrong with the theta space tracer using a mixed limiter' + + # Check for no new overshoots in A + assert np.max(final_fieldA.dat.data) <= np.max(true_fieldA.dat.data) + tol, \ + 'Application of the DG space limiter in the mixed limiter has not prevented overshoots' + + # Check for no new undershoots in A + assert np.min(final_fieldA.dat.data) >= np.min(true_fieldA.dat.data) - tol, \ + 'Application of the DG space limiter in the mixed limiter has not prevented undershoots' + + # Check for no new overshoots in B + assert np.max(final_fieldB.dat.data) <= np.max(true_fieldB.dat.data) + tol, \ + 'Application of the theta space limiter in the mixed limiter has not prevented overshoots' + + # Check for no new undershoots in B + assert np.min(final_fieldB.dat.data) >= np.min(true_fieldB.dat.data) - tol, \ + 'Application of the theta space limiter in the mixed limiter has not prevented undershoots' + + else: + final_field = stepper.fields('tracer') + + # Check tracer is roughly in the correct place + assert norm(true_field - final_field) / norm(true_field) < 0.05, \ + 'Something appears to have gone wrong with transport of tracer using a limiter' + + # Check for no new overshoots + assert np.max(final_field.dat.data) <= np.max(true_field.dat.data) + tol, \ + 'Application of limiter has not prevented overshoots' + + # Check for no new undershoots + assert np.min(final_field.dat.data) >= np.min(true_field.dat.data) - tol, \ + 'Application of limiter has not prevented undershoots' From 021d9d7f017830159c05015f55a4fc9d64afc59e Mon Sep 17 00:00:00 2001 From: Tim Andrews Date: Wed, 15 Nov 2023 01:00:46 +0000 Subject: [PATCH 04/10] changed mixed_fs_limiter test to DG and DG1_equispaced --- gusto/limiters.py | 37 +++++--------- integration-tests/transport/test_limiters.py | 51 ++++++++------------ 2 files changed, 32 insertions(+), 56 deletions(-) diff --git a/gusto/limiters.py b/gusto/limiters.py index b8d97382f..2011bc680 100644 --- a/gusto/limiters.py +++ b/gusto/limiters.py @@ -6,7 +6,7 @@ """ from firedrake import (BrokenElement, Function, FunctionSpace, interval, - FiniteElement, TensorProductElement, split) + FiniteElement, TensorProductElement) from firedrake.slope_limiter.vertex_based_limiter import VertexBasedLimiter from gusto.kernels import LimitMidpoints @@ -184,44 +184,29 @@ class MixedFSLimiter(object): """An object to hold a dictionary that defines limiters for transported prognostic variables. Different limiters may be applied to different fields and not every transported variable needs a defined limiter. """ - + def __init__(self, equation, sublimiters): """ Args: equation (:class: 'PrognosticEquationSet'): the prognostic equation(s) - sublimiters (dict): A dictionary holding limiters defined for individual prognostic variables + sublimiters (dict): A dictionary holding limiters defined for individual prognostic variables Raises: ValueError: If a limiter is defined for a field that is not in the prognostic variable set """ self.sublimiters = sublimiters - - for tracer, sublimiter in sublimiters.items(): - # Check that the tracer is being solved in the equation: - if tracer not in equation.field_names: - raise ValueError(f"The limiter for {tracer} is for a field that does not exist in the equation set") + + for field, sublimiter in sublimiters.items(): + # Check that the field is in the prognostic variable set: + if field not in equation.field_names: + raise ValueError(f"The limiter defined for {field} is for a field that does not exist in the equation set") else: - self.sublimiters[tracer].idx = equation.field_names.index(tracer) - - #for tracer in tracers: - # if tracer.name in sublimiters: - # sublimiter = sublimiters.get(tracer.name) - # print(sublimiter) - # self.sublimiters[tracer.name].idx = equation.field_names.index(tracer.name - + self.sublimiters[field].idx = equation.field_names.index(field) + def apply(self, fields): """ Apply the individual limiters to specific prognostic variables """ - - #for tracer in self.tracers: - # if tracer.name in self.sublimiters: - # sublimiter = self.sublimiters[tracer.name] - # field = fields.subfunctions[sublimiter.idx] - # sublimiter.apply(field) - # print(f'Applying sublimiter on {tracer.name} field') - + for _, sublimiter in self.sublimiters.items(): field = fields.subfunctions[sublimiter.idx] sublimiter.apply(field) - #print(f'Applying sublimiter on {tracer} field') - \ No newline at end of file diff --git a/integration-tests/transport/test_limiters.py b/integration-tests/transport/test_limiters.py index 0e507dcba..d992dfded 100644 --- a/integration-tests/transport/test_limiters.py +++ b/integration-tests/transport/test_limiters.py @@ -53,8 +53,7 @@ def setup_limiters(dirname, space): elif space == 'mixed_FS': V = domain.spaces('HDiv') VA = domain.spaces('DG') - #VB = domain.spaces('theta') - VB = domain.spaces('DG') + VB = domain.spaces('DG1_equispaced') else: raise NotImplementedError @@ -65,10 +64,7 @@ def setup_limiters(dirname, space): tracerA = ActiveTracer(name='tracerA', space='DG', variable_type=TracerVariableType.mixing_ratio, transport_eqn=TransportEquationType.advective) - #tracerB = ActiveTracer(name='tracerB', space='theta', - # variable_type=TracerVariableType.mixing_ratio, - # transport_eqn=TransportEquationType.advective) - tracerB = ActiveTracer(name='tracerB', space='DG', + tracerB = ActiveTracer(name='tracerB', space='DG1_equispaced', variable_type=TracerVariableType.mixing_ratio, transport_eqn=TransportEquationType.advective) tracers = [tracerA, tracerB] @@ -106,14 +102,9 @@ def setup_limiters(dirname, space): elif space == 'mixed_FS': sublimiters = {'tracerA': DG1Limiter(domain.spaces('DG')), - 'tracerB': DG1Limiter(domain.spaces('DG'))} - #'tracerB': ThetaLimiter(domain.spaces('theta'))} - # Need Embedded DG if wanting to test theta limiter. + 'tracerB': VertexBasedLimiter(domain.spaces('DG1_equispaced'))} MixedLimiter = MixedFSLimiter(eqn, sublimiters) transport_schemes = SSPRK3(domain, limiter=MixedLimiter) - - #theta_opts = EmbeddedDGOptions() - else: raise NotImplementedError @@ -234,8 +225,8 @@ def setup_limiters(dirname, space): return stepper, tmax, true_field #move mixed_FS to the end once finished debugging -@pytest.mark.parametrize('space', ['Vtheta_degree_0', 'mixed_FS'])#, 'Vtheta_degree_1', - #'DG0', 'DG1', 'DG1_equispaced']) +@pytest.mark.parametrize('space', ['Vtheta_degree_0', 'Vtheta_degree_1', 'DG0', + 'DG1', 'DG1_equispaced', 'mixed_FS']) def test_limiters(tmpdir, space): # Setup and run @@ -245,50 +236,50 @@ def test_limiters(tmpdir, space): stepper, tmax, true_fieldA, true_fieldB = setup_limiters(dirname, space) else: stepper, tmax, true_field = setup_limiters(dirname, space) - + stepper.run(t=0, tmax=tmax) - + tol = 1e-9 - + if space == 'mixed_FS': final_fieldA = stepper.fields('tracerA') final_fieldB = stepper.fields('tracerB') - + # Check tracer is roughly in the correct place assert norm(true_fieldA - final_fieldA) / norm(true_fieldA) < 0.05, \ 'Something is wrong with the DG space tracer using a mixed limiter' - + # Check tracer is roughly in the correct place assert norm(true_fieldB - final_fieldB) / norm(true_fieldB) < 0.05, \ - 'Something is wrong with the theta space tracer using a mixed limiter' - + 'Something is wrong with the DG1 equispaced tracer using a mixed limiter' + # Check for no new overshoots in A assert np.max(final_fieldA.dat.data) <= np.max(true_fieldA.dat.data) + tol, \ 'Application of the DG space limiter in the mixed limiter has not prevented overshoots' - + # Check for no new undershoots in A assert np.min(final_fieldA.dat.data) >= np.min(true_fieldA.dat.data) - tol, \ 'Application of the DG space limiter in the mixed limiter has not prevented undershoots' - + # Check for no new overshoots in B assert np.max(final_fieldB.dat.data) <= np.max(true_fieldB.dat.data) + tol, \ - 'Application of the theta space limiter in the mixed limiter has not prevented overshoots' - + 'Application of the DG1 equispaced limiter in the mixed limiter has not prevented overshoots' + # Check for no new undershoots in B assert np.min(final_fieldB.dat.data) >= np.min(true_fieldB.dat.data) - tol, \ - 'Application of the theta space limiter in the mixed limiter has not prevented undershoots' - + 'Application of the DG1 equispaced limiter in the mixed limiter has not prevented undershoots' + else: final_field = stepper.fields('tracer') - + # Check tracer is roughly in the correct place assert norm(true_field - final_field) / norm(true_field) < 0.05, \ 'Something appears to have gone wrong with transport of tracer using a limiter' - + # Check for no new overshoots assert np.max(final_field.dat.data) <= np.max(true_field.dat.data) + tol, \ 'Application of limiter has not prevented overshoots' - + # Check for no new undershoots assert np.min(final_field.dat.data) >= np.min(true_field.dat.data) - tol, \ 'Application of limiter has not prevented undershoots' From 4186efc22c01e810966e83f4286fc23d31a913a5 Mon Sep 17 00:00:00 2001 From: Tim Andrews Date: Thu, 16 Nov 2023 01:53:02 +0000 Subject: [PATCH 05/10] lint --- gusto/limiters.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gusto/limiters.py b/gusto/limiters.py index 2011bc680..41adde50d 100644 --- a/gusto/limiters.py +++ b/gusto/limiters.py @@ -181,7 +181,7 @@ def apply(self, field): class MixedFSLimiter(object): - """An object to hold a dictionary that defines limiters for transported prognostic + """An object to hold a dictionary that defines limiters for transported prognostic variables. Different limiters may be applied to different fields and not every transported variable needs a defined limiter. """ From 792bb45d3de950e76d0d1c4f0b282c9f86d5f0e4 Mon Sep 17 00:00:00 2001 From: Tim Andrews Date: Thu, 16 Nov 2023 02:41:12 +0000 Subject: [PATCH 06/10] lint --- integration-tests/transport/test_limiters.py | 25 ++++++++++---------- 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/integration-tests/transport/test_limiters.py b/integration-tests/transport/test_limiters.py index d992dfded..b4635cc1c 100644 --- a/integration-tests/transport/test_limiters.py +++ b/integration-tests/transport/test_limiters.py @@ -68,14 +68,14 @@ def setup_limiters(dirname, space): variable_type=TracerVariableType.mixing_ratio, transport_eqn=TransportEquationType.advective) tracers = [tracerA, tracerB] - eqn = CoupledTransportEquation(domain, active_tracers=tracers, Vu = V) - output = OutputParameters(dirname=dirname+'/limiters', - dumpfreq=1, dumplist=['u', 'tracerA', 'tracerB', 'true_tracerA', 'true_tracerB']) + eqn = CoupledTransportEquation(domain, active_tracers=tracers, Vu=V) + output = OutputParameters(dirname=dirname+'/limiters', dumpfreq=1, + dumplist=['u', 'tracerA', 'tracerB', 'true_tracerA', 'true_tracerB']) else: eqn = AdvectionEquation(domain, V, 'tracer') output = OutputParameters(dirname=dirname+'/limiters', - dumpfreq=1, dumplist=['u', 'tracer', 'true_tracer']) - + dumpfreq=1, dumplist=['u', 'tracer', 'true_tracer']) + io = IO(domain, output) # ------------------------------------------------------------------------ # @@ -99,7 +99,7 @@ def setup_limiters(dirname, space): elif space == 'Vtheta_degree_1': opts = EmbeddedDGOptions() transport_schemes = SSPRK3(domain, options=opts, limiter=ThetaLimiter(V)) - + elif space == 'mixed_FS': sublimiters = {'tracerA': DG1Limiter(domain.spaces('DG')), 'tracerB': VertexBasedLimiter(domain.spaces('DG1_equispaced'))} @@ -176,13 +176,13 @@ def setup_limiters(dirname, space): 0.0) if i == 0: - if space=='mixed_FS': + if space == 'mixed_FS': tracerA_0.interpolate(Constant(tracer_min) + expr_1 + expr_2) tracerB_0.interpolate(Constant(tracer_min) + expr_1 + expr_2) else: tracer0.interpolate(Constant(tracer_min) + expr_1 + expr_2) elif i == 1: - if space=='mixed_FS': + if space == 'mixed_FS': true_fieldA.interpolate(Constant(tracer_min) + expr_1 + expr_2) true_fieldB.interpolate(Constant(tracer_min) + expr_1 + expr_2) else: @@ -224,14 +224,15 @@ def setup_limiters(dirname, space): else: return stepper, tmax, true_field -#move mixed_FS to the end once finished debugging @pytest.mark.parametrize('space', ['Vtheta_degree_0', 'Vtheta_degree_1', 'DG0', 'DG1', 'DG1_equispaced', 'mixed_FS']) + + def test_limiters(tmpdir, space): # Setup and run dirname = str(tmpdir) - + if space == 'mixed_FS': stepper, tmax, true_fieldA, true_fieldB = setup_limiters(dirname, space) else: @@ -247,11 +248,11 @@ def test_limiters(tmpdir, space): # Check tracer is roughly in the correct place assert norm(true_fieldA - final_fieldA) / norm(true_fieldA) < 0.05, \ - 'Something is wrong with the DG space tracer using a mixed limiter' + 'Something is wrong with the DG space tracer using a mixed limiter' # Check tracer is roughly in the correct place assert norm(true_fieldB - final_fieldB) / norm(true_fieldB) < 0.05, \ - 'Something is wrong with the DG1 equispaced tracer using a mixed limiter' + 'Something is wrong with the DG1 equispaced tracer using a mixed limiter' # Check for no new overshoots in A assert np.max(final_fieldA.dat.data) <= np.max(true_fieldA.dat.data) + tol, \ From 126ce4e935dd17d88612b302db68c89ab2a992a1 Mon Sep 17 00:00:00 2001 From: Tim Andrews Date: Thu, 16 Nov 2023 02:46:44 +0000 Subject: [PATCH 07/10] lint --- integration-tests/transport/test_limiters.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/integration-tests/transport/test_limiters.py b/integration-tests/transport/test_limiters.py index b4635cc1c..9853b11c3 100644 --- a/integration-tests/transport/test_limiters.py +++ b/integration-tests/transport/test_limiters.py @@ -224,10 +224,9 @@ def setup_limiters(dirname, space): else: return stepper, tmax, true_field + @pytest.mark.parametrize('space', ['Vtheta_degree_0', 'Vtheta_degree_1', 'DG0', 'DG1', 'DG1_equispaced', 'mixed_FS']) - - def test_limiters(tmpdir, space): # Setup and run From c7edb1be7611739bfb3f3c2cc98b1f73ed812d4f Mon Sep 17 00:00:00 2001 From: Tim Andrews Date: Tue, 21 Nov 2023 18:27:26 +0000 Subject: [PATCH 08/10] removed fabs from minimum diagnostic --- gusto/diagnostics.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gusto/diagnostics.py b/gusto/diagnostics.py index 024a42fa9..f53d047d1 100644 --- a/gusto/diagnostics.py +++ b/gusto/diagnostics.py @@ -73,7 +73,7 @@ def min(f): fmin = op2.Global(1, np.finfo(float).max, dtype=float, comm=f._comm) op2.par_loop(op2.Kernel(""" static void minify(double *a, double *b) { - a[0] = a[0] > fabs(b[0]) ? fabs(b[0]) : a[0]; + a[0] = a[0] > b[0] ? b[0] : a[0]; } """, "minify"), f.dof_dset.set, fmin(op2.MIN), f.dat(op2.READ)) return fmin.data[0] @@ -91,7 +91,7 @@ def max(f): fmax = op2.Global(1, np.finfo(float).min, dtype=float, comm=f._comm) op2.par_loop(op2.Kernel(""" static void maxify(double *a, double *b) { - a[0] = a[0] < fabs(b[0]) ? fabs(b[0]) : a[0]; + a[0] = a[0] < b[0] ? b[0] : a[0]; } """, "maxify"), f.dof_dset.set, fmax(op2.MAX), f.dat(op2.READ)) return fmax.data[0] From dd18af9c539fc408bbb45d83b8546d6cac58c97a Mon Sep 17 00:00:00 2001 From: Tim Andrews Date: Wed, 6 Dec 2023 21:21:48 +0000 Subject: [PATCH 09/10] tom changes --- gusto/limiters.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/gusto/limiters.py b/gusto/limiters.py index bba56d10f..a39c6688d 100644 --- a/gusto/limiters.py +++ b/gusto/limiters.py @@ -181,14 +181,16 @@ def apply(self, field): class MixedFSLimiter(object): - """An object to hold a dictionary that defines limiters for transported prognostic - variables. Different limiters may be applied to different fields and not every transported variable needs a defined limiter. + """ + An object to hold a dictionary that defines limiters for transported prognostic + variables. Different limiters may be applied to different fields and not every + transported variable needs a defined limiter. """ def __init__(self, equation, sublimiters): """ Args: - equation (:class: 'PrognosticEquationSet'): the prognostic equation(s) + equation (:class: `PrognosticEquationSet`): the prognostic equation(s) sublimiters (dict): A dictionary holding limiters defined for individual prognostic variables Raises: ValueError: If a limiter is defined for a field that is not in the prognostic variable set From 3a7942d89d2ab365045e9f32f081d97d6f83c8ec Mon Sep 17 00:00:00 2001 From: Tim Andrews Date: Wed, 6 Dec 2023 21:36:31 +0000 Subject: [PATCH 10/10] revert diagnostics file --- gusto/diagnostics.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gusto/diagnostics.py b/gusto/diagnostics.py index ca3837016..b0969dfd8 100644 --- a/gusto/diagnostics.py +++ b/gusto/diagnostics.py @@ -73,7 +73,7 @@ def min(f): fmin = op2.Global(1, np.finfo(float).max, dtype=float, comm=f._comm) op2.par_loop(op2.Kernel(""" static void minify(double *a, double *b) { - a[0] = a[0] > b[0] ? b[0] : a[0]; + a[0] = a[0] > fabs(b[0]) ? fabs(b[0]) : a[0]; } """, "minify"), f.dof_dset.set, fmin(op2.MIN), f.dat(op2.READ)) return fmin.data[0] @@ -91,7 +91,7 @@ def max(f): fmax = op2.Global(1, np.finfo(float).min, dtype=float, comm=f._comm) op2.par_loop(op2.Kernel(""" static void maxify(double *a, double *b) { - a[0] = a[0] < b[0] ? b[0] : a[0]; + a[0] = a[0] < fabs(b[0]) ? fabs(b[0]) : a[0]; } """, "maxify"), f.dof_dset.set, fmax(op2.MAX), f.dat(op2.READ)) return fmax.data[0]