diff --git a/examples/compressible/dry_bryan_fritsch.py b/examples/compressible/dry_bryan_fritsch.py index a6dba5a22..653fc70bb 100644 --- a/examples/compressible/dry_bryan_fritsch.py +++ b/examples/compressible/dry_bryan_fritsch.py @@ -9,10 +9,8 @@ from firedrake import (IntervalMesh, ExtrudedMesh, SpatialCoordinate, conditional, cos, pi, sqrt, TestFunction, dx, TrialFunction, Constant, Function, - LinearVariationalProblem, LinearVariationalSolver, - FunctionSpace, VectorFunctionSpace) + LinearVariationalProblem, LinearVariationalSolver) import sys - # ---------------------------------------------------------------------------- # # Test case parameters # ---------------------------------------------------------------------------- # @@ -61,19 +59,14 @@ io = IO(domain, output, diagnostic_fields=diagnostic_fields) # Transport schemes -- set up options for using recovery wrapper -VDG1 = domain.spaces("DG1_equispaced") -VCG1 = FunctionSpace(mesh, "CG", 1) -Vu_DG1 = VectorFunctionSpace(mesh, VDG1.ufl_element()) -Vu_CG1 = VectorFunctionSpace(mesh, "CG", 1) - -u_opts = RecoveryOptions(embedding_space=Vu_DG1, - recovered_space=Vu_CG1, - boundary_method=BoundaryMethod.taylor) -rho_opts = RecoveryOptions(embedding_space=VDG1, - recovered_space=VCG1, - boundary_method=BoundaryMethod.taylor) -theta_opts = RecoveryOptions(embedding_space=VDG1, - recovered_space=VCG1) +boundary_methods = {'DG': BoundaryMethod.taylor, + 'HDiv': BoundaryMethod.taylor} + +recovery_spaces = RecoverySpaces(domain, boundary_methods, use_vector_spaces=True) + +u_opts = recovery_spaces.HDiv_options +rho_opts = recovery_spaces.DG_options +theta_opts = recovery_spaces.theta_options transported_fields = [SSPRK3(domain, "rho", options=rho_opts), SSPRK3(domain, "theta", options=theta_opts), diff --git a/examples/compressible/unsaturated_bubble.py b/examples/compressible/unsaturated_bubble.py index 074224cf8..23e38624f 100644 --- a/examples/compressible/unsaturated_bubble.py +++ b/examples/compressible/unsaturated_bubble.py @@ -6,11 +6,12 @@ Limiters are applied to the transport of the water species. """ from gusto import * +from gusto.equations import thermodynamics from firedrake import (PeriodicIntervalMesh, ExtrudedMesh, SpatialCoordinate, conditional, cos, pi, sqrt, exp, TestFunction, dx, TrialFunction, Constant, Function, LinearVariationalProblem, LinearVariationalSolver, - FunctionSpace, VectorFunctionSpace, errornorm) + errornorm) from firedrake.slope_limiter.vertex_based_limiter import VertexBasedLimiter import sys @@ -61,21 +62,17 @@ diagnostic_fields = [RelativeHumidity(eqns), Perturbation('theta'), Perturbation('water_vapour'), Perturbation('rho'), Perturbation('RelativeHumidity')] io = IO(domain, output, diagnostic_fields=diagnostic_fields) - # Transport schemes -- specify options for using recovery wrapper +boundary_methods = {'DG': BoundaryMethod.taylor, + 'HDiv': BoundaryMethod.taylor} + +recovery_spaces = RecoverySpaces(domain, boundary_method=boundary_methods, use_vector_spaces=True) + +u_opts = recovery_spaces.HDiv_options +rho_opts = recovery_spaces.DG_options +theta_opts = recovery_spaces.theta_options + VDG1 = domain.spaces("DG1_equispaced") -VCG1 = FunctionSpace(mesh, "CG", 1) -Vu_DG1 = VectorFunctionSpace(mesh, VDG1.ufl_element()) -Vu_CG1 = VectorFunctionSpace(mesh, "CG", 1) -Vt = domain.spaces("theta") - -u_opts = RecoveryOptions(embedding_space=Vu_DG1, - recovered_space=Vu_CG1, - boundary_method=BoundaryMethod.taylor) -rho_opts = RecoveryOptions(embedding_space=VDG1, - recovered_space=VCG1, - boundary_method=BoundaryMethod.taylor) -theta_opts = RecoveryOptions(embedding_space=VDG1, recovered_space=VCG1) limiter = VertexBasedLimiter(VDG1) transported_fields = [SSPRK3(domain, "u", options=u_opts), @@ -92,6 +89,7 @@ # Physics schemes # NB: can't yet use wrapper or limiter options with physics +Vt = domain.spaces('theta') rainfall_method = DGUpwind(eqns, 'rain', outflow=True) zero_limiter = MixedFSLimiter(eqns, {'water_vapour': ZeroLimiter(Vt), 'cloud_water': ZeroLimiter(Vt)}) diff --git a/gusto/recovery/__init__.py b/gusto/recovery/__init__.py index f12cb298d..16e23ba37 100644 --- a/gusto/recovery/__init__.py +++ b/gusto/recovery/__init__.py @@ -1,3 +1,4 @@ from gusto.recovery.averaging import * # noqa from gusto.recovery.recovery import * # noqa -from gusto.recovery.reversible_recovery import * # noqa \ No newline at end of file +from gusto.recovery.reversible_recovery import * # noqa +from gusto.recovery.recovery_options import * # noqa \ No newline at end of file diff --git a/gusto/recovery/recovery_options.py b/gusto/recovery/recovery_options.py new file mode 100644 index 000000000..936e30213 --- /dev/null +++ b/gusto/recovery/recovery_options.py @@ -0,0 +1,117 @@ +"""Contains an object to generate the set of spaces used for the recovery wrapper """ +from firedrake import (interval, FiniteElement, TensorProductElement, FunctionSpace, + VectorFunctionSpace) +from gusto.core.function_spaces import DeRhamComplex +from gusto.core.configuration import RecoveryOptions + + +class RecoverySpaces(object): + """ + Finds or builds necessary spaces to carry out recovery transport for lowest + and mixed order domains (0,0), (0,1) and (1,0) + """ + + def __init__(self, domain, boundary_method=None, use_vector_spaces=False): + """ + Args: + domain (:class:`Domain`): the model's domain object, containing the + mesh and the compatible function spaces. + + boundary_method (:variable:'dict', optional): A dictionary containing the space + the boundary method is to be applied to along with specified method. Acceptable keys are "DG", + "HDiv" and "theta". acceptable values are (BoundaryMethod.taylor/hcurl/extruded), + passed as ('space', 'boundary method'). Defaults to None + + use_vector_spaces (bool, optional):. Determines if we need to use DG / CG + space for the embedded and recovery space for the HDiv field instead of the usual + HDiv, HCurl spaces. Defaults to False + """ + family = domain.family + mesh = domain.mesh + # Need spaces from current deRham and a higher order deRham + self.de_Rham = DeRhamComplex(mesh, family, + horizontal_degree=1, + vertical_degree=1, + complex_name='recovery_de_Rham') + + valid_keys = ['DG', 'HDiv', 'theta'] + if boundary_method is not None: + for key in boundary_method: + if key not in valid_keys: + raise KeyError(f'Recovery spaces: boundary method key {key} not valid. Valid keys are DG, HDiv, theta') + + # ---------------------------------------------------------------------- + # Building theta options if on an extruded mesh + # ---------------------------------------------------------------------- + + # Check if extruded and if so builds theta spaces + if hasattr(mesh, "_base_mesh"): + # check if boundary method is present + if hasattr(boundary_method, 'theta'): + theta_boundary_method = boundary_method['theta'] + else: + theta_boundary_method = None + cell = mesh._base_mesh.ufl_cell().cellname() + DG_hori_ele = FiniteElement('DG', cell, 1, variant='equispaced') + DG_vert_ele = FiniteElement('DG', interval, (domain.vertical_degree + 1), variant='equispaced') + CG_hori_ele = FiniteElement('CG', cell, 1) + CG_vert_ele = FiniteElement('CG', interval, (domain.vertical_degree + 1)) + + VDG_ele = TensorProductElement(DG_hori_ele, DG_vert_ele) + VCG_ele = TensorProductElement(CG_hori_ele, CG_vert_ele) + VDG_theta = FunctionSpace(mesh, VDG_ele) + VCG_theta = FunctionSpace(mesh, VCG_ele) + + self.theta_options = RecoveryOptions(embedding_space=VDG_theta, + recovered_space=VCG_theta, + boundary_method=theta_boundary_method) + else: + cell = self.mesh.ufl_cell().cellname() + + # ---------------------------------------------------------------------- + # Building the DG options + # ---------------------------------------------------------------------- + if hasattr(boundary_method, 'DG'): + DG_boundary_method = boundary_method['DG'] + else: + DG_boundary_method = None + + DG_embedding_space = domain.spaces.DG1_equispaced + # Recovered space needs builing manually to avoid uneccesary DoFs + CG_hori_ele_DG = FiniteElement('CG', cell, 1) + CG_vert_ele_DG = FiniteElement('CG', interval, 1) + VCG_ele_DG = TensorProductElement(CG_hori_ele_DG, CG_vert_ele_DG) + DG_recovered_space = FunctionSpace(mesh, VCG_ele_DG) + + # DG_recovered_space = domain.spaces.H1 + self.DG_options = RecoveryOptions(embedding_space=DG_embedding_space, + recovered_space=DG_recovered_space, + boundary_method=DG_boundary_method) + # ---------------------------------------------------------------------- + # Building HDiv options + # ---------------------------------------------------------------------- + + if hasattr(boundary_method, 'HDiv'): + HDiv_boundary_method = boundary_method['HDiv'] + else: + HDiv_boundary_method = None + + if use_vector_spaces: + Vu_DG1 = VectorFunctionSpace(mesh, DG_embedding_space.ufl_element()) + Vu_CG1 = VectorFunctionSpace(mesh, "CG", 1) + + HDiv_embedding_Space = Vu_DG1 + HDiv_recovered_Space = Vu_CG1 + + else: + + HDiv_embedding_Space = self.de_Rham.HDiv + HDiv_recovered_Space = self.de_Rham.HCurl + + self.HDiv_options = RecoveryOptions(embedding_space=HDiv_embedding_Space, + recovered_space=HDiv_recovered_Space, + injection_method='recover', + project_high_method='project', + project_low_method='project', + broken_method='project', + boundary_method=HDiv_boundary_method) diff --git a/unit-tests/recovery_tests/test_mixed_spaces.py b/unit-tests/recovery_tests/test_mixed_spaces.py new file mode 100644 index 000000000..44d9aa546 --- /dev/null +++ b/unit-tests/recovery_tests/test_mixed_spaces.py @@ -0,0 +1,32 @@ +""" +Test that the Recovery_options object builds spaces of the correct degree +""" +from firedrake import UnitIntervalMesh, ExtrudedMesh +from gusto.core import Domain +from gusto.recovery.recovery_options import RecoverySpaces +import pytest + +order_correct_degree_dict = {(0, 0): {'DG': (1, 1), + 'HDiv': (2, 2), + 'theta': (1, 1)}, + (0, 1): {'DG': (1, 1), + 'HDiv': (2, 2), + 'theta': (1, 2)}, + (1, 0): {'DG': (1, 1), + 'HDiv': (2, 2), + 'theta': (1, 1)} + } + + +@pytest.mark.parametrize("order", [(0, 0), (1, 0), (0, 1)]) +@pytest.mark.parametrize('space', ['DG', 'HDiv', 'theta']) +def test_mixed_spaces(order, space): + mesh = UnitIntervalMesh(3) + emesh = ExtrudedMesh(mesh, 3, 0.33) + dt = 1 + domain = Domain(emesh, dt, family='CG', horizontal_degree=order[0], vertical_degree=order[1]) + recovery_spaces = RecoverySpaces(domain) + + tesing_space = getattr(recovery_spaces, f'{space}_options') + degree = tesing_space.recovered_space.finat_element.degree + assert degree == order_correct_degree_dict[order][space]