Merge pull request #464 from firedrakeproject/ufl_update_fixes

Fix issues from update to ufl (deleting brackets!)

gusto/coordinates.py

@@ -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 '

gusto/kernels.py

@@ -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
                   """)

gusto/limiters.py

@@ -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')
@@ -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)):

gusto/preconditioners.py

@@ -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
@@ -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)

gusto/recovery/recovery.py

@@ -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)
@@ -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]
 

gusto/wrappers.py

@@ -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):
@@ -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)]

integration-tests/model/test_checkpointing.py

@@ -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'
 
     # ------------------------------------------------------------------------ #