-
Notifications
You must be signed in to change notification settings - Fork 11
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
PR #524: from firedrakeproject/mixed-order-recovery
Implement an object to encapsulate the choice of recovered spaces
- Loading branch information
Showing
5 changed files
with
172 additions
and
31 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,3 +1,4 @@ | ||
| from gusto.recovery.averaging import * # noqa | ||
| from gusto.recovery.recovery import * # noqa | ||
| from gusto.recovery.reversible_recovery import * # noqa | ||
| from gusto.recovery.reversible_recovery import * # noqa | ||
| from gusto.recovery.recovery_options import * # noqa |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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] |