Reduce Firedrake deprectation warnings

1) Don't use Constant(value, domain=mesh). New utility function
   domain_constant(value, mesh) that creates Function in R space
2) Don't use interpolate(expr, V) (without assemble) instead use
   Function(V).interpolate(expr)

demos/demo_2d_multiple_tracers.py

@@ -45,9 +45,8 @@
 options.use_lax_friedrichs_tracer = False
 options.use_limiter_for_tracers = False
 
-vP1_2d = VectorFunctionSpace(mesh2d, "CG", 1)
 x, y = SpatialCoordinate(mesh2d)
-uv_init = interpolate(as_vector([0.5 - y, x - 0.5]), vP1_2d)
+uv_init = as_vector([0.5 - y, x - 0.5])
 
 # Initial conditions for each tracer are defined as before,
 # but must be assigned separately. ::
@@ -65,9 +64,9 @@
     0.0
 )
 
-bell_init = interpolate(1.0 + bell, P1_2d)
-cone_init = interpolate(1.0 + cone, P1_2d)
-slot_cyl_init = interpolate(1.0 + slot_cyl, P1_2d)
+bell_init = Function(P1_2d).interpolate(1.0 + bell)
+cone_init = Function(P1_2d).interpolate(1.0 + cone)
+slot_cyl_init = Function(P1_2d).interpolate(1.0 + slot_cyl)
 solver_obj.assign_initial_conditions(
     uv=uv_init, bell_2d=bell_init, cone_2d=cone_init, slot_cyl_2d=slot_cyl_init
 )

demos/demo_2d_tracer.py

@@ -87,9 +87,7 @@
 
 # The velocity field is set up using a simple analytic expression. ::
 
-vP1_2d = VectorFunctionSpace(mesh2d, "CG", 1)
-x, y = SpatialCoordinate(mesh2d)
-uv_init = interpolate(as_vector([0.5 - y, x - 0.5]), vP1_2d)
+uv_init = as_vector([0.5 - y, x - 0.5])
 
 # Now, we set up the cosine-bell--cone--slotted-cylinder initial condition. The
 # first four lines declare various parameters relating to the positions of these
@@ -115,7 +113,7 @@
 # neglecting the inflow boundary condition.  We also save the initial state so
 # that we can check the :math:`L^2`-norm error at the end. ::
 
-q_init = interpolate(1.0 + bell + cone + slot_cyl, P1_2d)
+q_init = Function(P1_2d).interpolate(1.0 + bell + cone + slot_cyl)
 solver_obj.assign_initial_conditions(uv=uv_init, tracer_2d=q_init)
 
 # Now we are in a position to run the time loop. ::

examples/baroclinic_eddies/diagnostics.py

@@ -28,8 +28,8 @@ class RPECalculator(DiagnosticCallback):
 
     def _initialize(self):
         # compute area of 2D mesh
-        one_2d = Constant(1.0, domain=self.solver_obj.mesh2d.coordinates.ufl_domain())
-        self.area_2d = assemble(one_2d*dx)
+        one_2d = Constant(1.0)
+        self.area_2d = assemble(one_2d*dx(domain=self.solver_obj.mesh2d))
         self.rho = self.solver_obj.fields.density_3d
         fs = self.rho.function_space()
         self.test = TestFunction(fs)

examples/channel_inversion/inverse_problem.py

@@ -93,7 +93,7 @@
 sta_manager.set_model_field(solver_obj.fields.elev_2d)
 
 # Define the scaling for the cost function so that J ~ O(1)
-J_scalar = Constant(solver_obj.dt / options.simulation_end_time, domain=mesh2d)
+J_scalar = domain_constant(solver_obj.dt / options.simulation_end_time, mesh2d)
 
 # Create inversion manager and add controls
 inv_manager = inversion_tools.InversionManager(

examples/discrete_turbines/channel-optimisation.py

@@ -93,7 +93,7 @@
 r = farm_options.turbine_options.diameter/2.
 
 # a list contains the coordinates of all turbines: regular, non-staggered 4 x 2 layout
-farm_options.turbine_coordinates = [[Constant(x+cos(y), domain=mesh2d), Constant(y+numpy.sin(x), domain=mesh2d)]
+farm_options.turbine_coordinates = [[domain_constant(x+cos(y), mesh2d), domain_constant(y+numpy.sin(x), mesh2d)]
                                     for x in np.linspace(site_x_start + 4*r, site_x_start + site_x - 4*r, 4)
                                     for y in np.linspace(site_y_start + 0.5*site_y-2*r,
                                                          site_y_start + 0.5*site_y+2*r, 2)]
@@ -199,11 +199,11 @@ def perturbation(r):
         # ensure we use the same perturbation on all processes when parallel:
         return mesh2d.comm.bcast(dx, 0)
 
-    m0 = [Constant(float(x) + perturbation(r), domain=mesh2d) for xy in farm_options.turbine_coordinates for x in xy]
+    m0 = [domain_constant(float(x) + perturbation(r), mesh2d) for xy in farm_options.turbine_coordinates for x in xy]
 
     # the perturbation over which we test the Taylor approximation
     # (the taylor test below starts with a 1/100th of that, followed by a series of halvings
-    h0 = [Constant(perturbation(1), domain=mesh2d) for xy in farm_options.turbine_coordinates for x in xy]
+    h0 = [domain_constant(perturbation(1), mesh2d) for xy in farm_options.turbine_coordinates for x in xy]
 
     # this tests whether the above Taylor series residual indeed converges to zero at 2nd order in h as h->0
     minconv = taylor_test(rf, m0, h0)

examples/lockExchange/diagnostics.py

@@ -18,8 +18,7 @@ class FrontLocationCalculator(DiagnosticCallback):
     variable_names = ['front_bot', 'front_top']
 
     def _initialize(self):
-        one_2d = Constant(1.0, domain=self.solver_obj.mesh2d.coordinates.ufl_domain())
-        self.area_2d = assemble(one_2d*dx)
+        self.area_2d = assemble(Constant(1.0)*dx(domain=self.solver_obj.mesh2d))
         # store density difference
         self.rho = self.solver_obj.fields.density_3d
         self.rho_lim = [self.rho.dat.data.min(), self.rho.dat.data.max()]
@@ -105,8 +104,8 @@ class RPECalculator(DiagnosticCallback):
 
     def _initialize(self):
         # compute area of 2D mesh
-        one_2d = Constant(1.0, domain=self.solver_obj.mesh2d.coordinates.ufl_domain())
-        self.area_2d = assemble(one_2d*dx)
+        one_2d = Constant(1.0)
+        self.area_2d = assemble(one_2d*dx(domain=self.solver_obj.mesh2d))
         self.rho = self.solver_obj.fields.density_3d
         fs = self.rho.function_space()
         self.test = TestFunction(fs)

test/operations/test_diagnostics.py

@@ -27,7 +27,7 @@ def test_hessian_recovery2d(interp, tol=1.0e-08):
     bowl = 0.5*(x**2 + y**2)
     if interp:
         P2_2d = get_functionspace(mesh2d, 'CG', 2)
-        bowl = interpolate(bowl, P2_2d)
+        bowl = Function(P2_2d).interpolate(bowl)
     P1v_2d = get_functionspace(mesh2d, 'CG', 1, vector=True)
     P1t_2d = get_functionspace(mesh2d, 'CG', 1, tensor=True)
 
@@ -63,7 +63,7 @@ def test_vorticity_calculation2d(interp, tol=1.0e-08):
     uv = 0.5*as_vector([y, -x])
     if interp:
         P1v_2d = get_functionspace(mesh2d, 'CG', 1, vector=True)
-        uv = interpolate(uv, P1v_2d)
+        uv = Function(P1v_2d).interpolate(uv)
     P1_2d = get_functionspace(mesh2d, 'CG', 1)
 
     # Recover vorticity

test/swe2d/test_rossby_wave.py

@@ -39,7 +39,7 @@ def asymptotic_expansion_uv(U_2d, order=1, time=0.0, soliton_amplitude=0.395):
     u_terms = phi*0.25*(-9 + 6*y*y)*psi
     v_terms = 2*y*dphidx*psi
     if order == 0:
-        return interpolate(as_vector([u_terms, v_terms]), U_2d)
+        return Function(U_2d).interpolate(as_vector([u_terms, v_terms]))
 
     # Unnormalised Hermite series coefficients for u
     u = numpy.zeros(28)
@@ -83,7 +83,7 @@ def asymptotic_expansion_uv(U_2d, order=1, time=0.0, soliton_amplitude=0.395):
     u_terms += phi*phi*psi*sum(u[i]*polynomials[i] for i in range(28))
     v_terms += dphidx*phi*psi*sum(v[i]*polynomials[i] for i in range(28))
 
-    return interpolate(as_vector([u_terms, v_terms]), U_2d)
+    return Function(U_2d).interpolate(as_vector([u_terms, v_terms]))
 
 
 def asymptotic_expansion_elev(H_2d, order=1, time=0.0, soliton_amplitude=0.395):
@@ -105,7 +105,7 @@ def asymptotic_expansion_elev(H_2d, order=1, time=0.0, soliton_amplitude=0.395):
     # Zeroth order terms
     eta_terms = phi*0.25*(3 + 6*y*y)*psi
     if order == 0:
-        return interpolate(eta_terms, H_2d)
+        return Function(H_2d).interpolate(eta_terms)
 
     # Unnormalised Hermite series coefficients for eta
     eta = numpy.zeros(28)
@@ -133,7 +133,7 @@ def asymptotic_expansion_elev(H_2d, order=1, time=0.0, soliton_amplitude=0.395):
     eta_terms += C*phi*0.5625*(-5 + 2*y*y)*psi
     eta_terms += phi*phi*psi*sum(eta[i]*polynomials[i] for i in range(28))
 
-    return interpolate(eta_terms, H_2d)
+    return Function(H_2d).interpolate(eta_terms)
 
 
 def run(refinement_level, **model_options):
@@ -172,7 +172,8 @@ def run(refinement_level, **model_options):
     options.use_grad_depth_viscosity_term = False
     options.horizontal_viscosity = None
     solver_obj.create_function_spaces()
-    options.coriolis_frequency = interpolate(y, solver_obj.function_spaces.P1_2d)
+    P1_2d = solver_obj.function_spaces.P1_2d
+    options.coriolis_frequency = Function(P1_2d).interpolate(y)
     options.swe_timestepper_options.solver_parameters['ksp_rtol'] = 1.0e-04
     options.no_exports = True
     options.fields_to_export = ['uv_2d', 'elev_2d', 'vorticity_2d']
@@ -202,8 +203,8 @@ def run(refinement_level, **model_options):
     physical_constants['g_grav'].assign(g)  # Revert g_grav value
 
     # Get mean peak heights
-    elev = interpolate(sign(y)*solver_obj.fields.elev_2d, P1_2d)  # Flip sign in southern hemisphere
-    xcoords = interpolate(mesh2d.coordinates[0], P1_2d)
+    elev = Function(P1_2d).interpolate(sign(y)*solver_obj.fields.elev_2d)  # Flip sign in southern hemisphere
+    xcoords = Function(P1_2d).interpolate(mesh2d.coordinates[0])
     with elev.dat.vec_ro as v:
         i_n, h_n = v.max()
         i_s, h_s = v.min()

test/swe2d/test_steady_state_channel.py

@@ -58,7 +58,7 @@ def test_steady_state_channel(do_export=False):
     uv, eta = solver_obj.fields.solution_2d.subfunctions
 
     x_2d, y_2d = SpatialCoordinate(mesh2d)
-    eta_ana = interpolate(1 - x_2d/lx, p1_2d)
+    eta_ana = Function(p1_2d).interpolate(1 - x_2d/lx)
     area = lx*ly
     l2norm = errornorm(eta_ana, eta)/math.sqrt(area)
     print_output(l2norm)

test/tracerEq/test_bcs_2d.py

@@ -38,7 +38,7 @@ def fourier_series_solution(mesh, lx, diff_flux, **model_options):
     # Initial condition and source term for two homogeneous Neumann problems
     P1 = get_functionspace(mesh, 'CG', 1)
     ic = Function(P1).interpolate(diff_flux*0.5*(lx - x)*(lx - x)/lx)
-    source = Constant(-nu*diff_flux/lx, domain=mesh)
+    source = domain_constant(-nu*diff_flux/lx, mesh)
 
     # The solution uses truncated Fourier expansions, meaning we need the following...
 

thetis/inversion_tools.py

@@ -3,7 +3,7 @@
 import ufl
 from .configuration import FrozenHasTraits
 from .solver2d import FlowSolver2d
-from .utility import create_directory, print_function_value_range, get_functionspace, unfrozen
+from .utility import create_directory, print_function_value_range, get_functionspace, unfrozen, domain_constant
 from .log import print_output
 from .diagnostics import HessianRecoverer2D
 from .exporter import HDF5Exporter
@@ -476,7 +476,7 @@ def construct_evaluator(self):
         self.mod_values_0d = fd.Function(self.fs_points_0d, name='model values')
         self.indicator_0d = fd.Function(self.fs_points_0d, name='station use indicator')
         self.indicator_0d.assign(1.0)
-        self.cost_function_scaling_0d = fd.Constant(0.0, domain=mesh0d)
+        self.cost_function_scaling_0d = domain_constant(0.0, mesh0d)
         self.cost_function_scaling_0d.assign(self.cost_function_scaling)
         self.station_weight_0d = fd.Function(self.fs_points_0d, name='station-wise weighting')
         self.station_weight_0d.assign(1.0)
@@ -611,7 +611,7 @@ def __init__(self, function, scaling=1.0):
         self.regularization_expr = 0
         self.mesh = function.function_space().mesh()
         # calculate mesh area (to scale the cost function)
-        self.mesh_area = fd.assemble(fd.Constant(1.0, domain=self.mesh) * fd.dx)
+        self.mesh_area = fd.assemble(fd.Constant(1.0) * fd.dx(domain=self.mesh))
         self.name = function.name()
 
     def eval_cost_function(self):

thetis/utility.py

@@ -136,6 +136,16 @@ def __setattr__(self, key, value):
         super(FieldDict, self).__setattr__(key, value)
 
 
+def domain_constant(value, mesh, **kwargs):
+    """Create constant over a domain
+
+    Returns what used to be Constant(mesh, domain=mesh)"""
+    R = FunctionSpace(mesh, "R", 0)
+    c = Function(R, **kwargs)
+    c.assign(value)
+    return c
+
+
 def get_functionspace(mesh, h_family, h_degree, v_family=None, v_degree=None,
                       vector=False, tensor=False, hdiv=False, variant=None, v_variant=None,
                       **kwargs):
@@ -404,8 +414,7 @@ def comp_volume_2d(eta, bath):
 @PETSc.Log.EventDecorator("thetis.comp_volume_3d")
 def comp_volume_3d(mesh):
     """Computes volume of the 3D domain as an integral"""
-    one = Constant(1.0, domain=mesh.coordinates.ufl_domain())
-    val = assemble(one*dx)
+    val = assemble(Constant(1.0)*dx(domain=mesh))
     return val