Skip to content
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

Fix issues from update to ufl (deleting brackets!) #464

Merged
merged 6 commits into from
Nov 28, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion gusto/coordinates.py
Original file line number Diff line number Diff line change
Expand Up @@ -99,7 +99,7 @@ def register_space(self, domain, space_name):
space = domain.spaces(space_name)

# Use the appropriate scalar function space if the space is vector
if np.prod(space.ufl_element().value_shape()) > 1:
if np.prod(space.ufl_element().value_shape) > 1:
# TODO: get scalar space, and only compute coordinates if necessary
logger.warning(f'Space {space_name} has more than one dimension, '
+ 'and coordinates used for netCDF output have not '
Expand Down
12 changes: 6 additions & 6 deletions gusto/kernels.py
Original file line number Diff line number Diff line change
Expand Up @@ -41,16 +41,16 @@ def __init__(self, Vt_brok):
<float64> min_value = 0.0
for i
for j
field_hat[i*3+2*j] = field_DG1[i*2+j]
field_hat[i*3+j] = field_DG1[i*2+j]
end
max_value = fmax(field_DG1[i*2], field_DG1[i*2+1])
min_value = fmin(field_DG1[i*2], field_DG1[i*2+1])
if field_old[i*3+1] > max_value
field_hat[i*3+1] = 0.5 * (field_DG1[i*2] + field_DG1[i*2+1])
elif field_old[i*3+1] < min_value
field_hat[i*3+1] = 0.5 * (field_DG1[i*2] + field_DG1[i*2+1])
if field_old[i*3+2] > max_value
field_hat[i*3+2] = 0.5 * (field_DG1[i*2] + field_DG1[i*2+1])
elif field_old[i*3+2] < min_value
field_hat[i*3+2] = 0.5 * (field_DG1[i*2] + field_DG1[i*2+1])
else
field_hat[i*3+1] = field_old[i*3+1]
field_hat[i*3+2] = field_old[i*3+2]
end
end
""")
Expand Down
4 changes: 2 additions & 2 deletions gusto/limiters.py
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@ def __init__(self, space, subspace=None):

# check that space is DG1
degree = space.ufl_element().degree()
if (space.ufl_element().sobolev_space().name != 'L2'
if (space.ufl_element().sobolev_space.name != 'L2'
or ((type(degree) is tuple and np.any([deg != 1 for deg in degree]))
and degree != 1)):
raise ValueError('DG1 limiter can only be applied to DG1 space')
Expand Down Expand Up @@ -111,7 +111,7 @@ def __init__(self, space):
raise ValueError('The Theta Limiter can only be used on an extruded mesh')

# check that horizontal degree is 1 and vertical degree is 2
sub_elements = space.ufl_element().sub_elements()
sub_elements = space.ufl_element().sub_elements
if (sub_elements[0].family() not in ['Discontinuous Lagrange', 'DQ']
or sub_elements[1].family() != 'Lagrange'
or space.ufl_element().degree() != (1, 2)):
Expand Down
6 changes: 3 additions & 3 deletions gusto/preconditioners.py
Original file line number Diff line number Diff line change
Expand Up @@ -72,7 +72,7 @@ def initialize(self, pc):
for i, Vi in enumerate(V):

# Vector-valued spaces will have a non-empty value_shape
if Vi.ufl_element().value_shape():
if Vi.ufl_element().value_shape:
self.vidx = i
else:
self.pidx = i
Expand All @@ -88,11 +88,11 @@ def initialize(self, pc):
deg, _ = Vv.ufl_element().degree()

# Assumes a tensor product cell (quads, triangular-prisms, cubes)
if not isinstance(Vp.ufl_element().cell(), TensorProductCell):
if not isinstance(Vp.ufl_element().cell, TensorProductCell):
raise NotImplementedError("Currently only implemented for tensor product discretizations")

# Only want the horizontal cell
cell, _ = Vp.ufl_element().cell()._cells
cell, _ = Vp.ufl_element().cell._cells

DG = FiniteElement("DG", cell, deg)
CG = FiniteElement("CG", interval, 1)
Expand Down
10 changes: 5 additions & 5 deletions gusto/recovery/recovery.py
Original file line number Diff line number Diff line change
Expand Up @@ -207,7 +207,7 @@ def __init__(self, x_in, x_out, method='interpolate', boundary_method=None):
if self.vector_function_space:
# VectorElement has to be on the outside
# so first need to get underlying finite element
brok_elt = VectorElement(BrokenElement(rec_elt.sub_elements()[0]))
brok_elt = VectorElement(BrokenElement(rec_elt.sub_elements[0]))
else:
# Otherwise we can immediately get broken element
brok_elt = BrokenElement(rec_elt)
Expand Down Expand Up @@ -332,20 +332,20 @@ def find_eff_coords(V0):
vec_DG1 = VectorFunctionSpace(mesh, DG1_element)
x = SpatialCoordinate(mesh)

if isinstance(V0.ufl_element(), VectorElement) or V0.ufl_element().value_size() > 1:
if isinstance(V0.ufl_element(), VectorElement) or V0.ufl_element().value_size > 1:
eff_coords_list = []
V0_coords_list = []

# treat this separately for each component
for i in range(V0.ufl_element().value_size()):
for i in range(V0.ufl_element().value_size):
# fill an d-dimensional list with i-th coordinate
x_list = [x[i] for j in range(V0.ufl_element().value_size())]
x_list = [x[i] for j in range(V0.ufl_element().value_size)]

# the i-th element in V0_coords_list is a vector with all components the i-th coord
ith_V0_coords = Function(V0).project(as_vector(x_list))
V0_coords_list.append(ith_V0_coords)

for i in range(V0.ufl_element().value_size()):
for i in range(V0.ufl_element().value_size):
# slice through V0_coords_list to obtain the coords of the DOFs for that component
x_list = [V0_coords[i] for V0_coords in V0_coords_list]

Expand Down
6 changes: 3 additions & 3 deletions gusto/wrappers.py
Original file line number Diff line number Diff line change
Expand Up @@ -253,9 +253,9 @@ def is_cg(V):
if isinstance(ele, BrokenElement):
return False
elif type(ele) == VectorElement:
return all([e.sobolev_space().name == "H1" for e in ele._sub_elements])
return all([e.sobolev_space.name == "H1" for e in ele._sub_elements])
else:
return V.ufl_element().sobolev_space().name == "H1"
return V.ufl_element().sobolev_space.name == "H1"


class SUPGWrapper(Wrapper):
Expand Down Expand Up @@ -293,7 +293,7 @@ def setup(self):
if is_cg(self.function_space):
vals = default_vals
else:
space = self.function_space.ufl_element().sobolev_space()
space = self.function_space.ufl_element().sobolev_space
if space.name in ["HDiv", "DirectionalH"]:
vals = [default_vals[i] if space[i].name == "H1"
else 0. for i in range(dim)]
Expand Down
Binary file modified integration-tests/data/dry_compressible_chkpt.h5
Binary file not shown.
Binary file modified integration-tests/data/linear_sw_wave_chkpt.h5
Binary file not shown.
Binary file modified integration-tests/data/moist_compressible_chkpt.h5
Binary file not shown.
2 changes: 1 addition & 1 deletion integration-tests/model/test_checkpointing.py
Original file line number Diff line number Diff line change
Expand Up @@ -175,7 +175,7 @@ def test_checkpointing(tmpdir, stepper_type, checkpoint_method):
for field_name in ['rho', 'theta', 'u']:
diff_array = stepper_2.fields(field_name).dat.data - stepper_3.fields(field_name).dat.data
error = np.linalg.norm(diff_array) / np.linalg.norm(stepper_2.fields(field_name).dat.data)
assert error < 5e-16, \
assert error < 5e-15, \
f'Checkpointed and picked up field {field_name} is not equal'

# ------------------------------------------------------------------------ #
Expand Down
Loading