From 6eec7a13bb1733d774d11f697d7772370df115ba Mon Sep 17 00:00:00 2001 From: Stephan Kramer Date: Fri, 7 Jun 2024 11:54:43 +0100 Subject: [PATCH] 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 | 9 ++++----- demos/demo_2d_tracer.py | 6 ++---- examples/baroclinic_eddies/diagnostics.py | 4 ++-- examples/channel_inversion/inverse_problem.py | 2 +- .../discrete_turbines/channel-optimisation.py | 6 +++--- examples/lockExchange/diagnostics.py | 7 +++---- test/operations/test_diagnostics.py | 4 ++-- test/swe2d/test_rossby_wave.py | 15 ++++++++------- test/swe2d/test_steady_state_channel.py | 2 +- test/tracerEq/test_bcs_2d.py | 2 +- thetis/inversion_tools.py | 6 +++--- thetis/utility.py | 13 +++++++++++-- 12 files changed, 41 insertions(+), 35 deletions(-) diff --git a/demos/demo_2d_multiple_tracers.py b/demos/demo_2d_multiple_tracers.py index 4cfdd3ea9..c60f3138a 100644 --- a/demos/demo_2d_multiple_tracers.py +++ b/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 ) diff --git a/demos/demo_2d_tracer.py b/demos/demo_2d_tracer.py index fa4694f9c..2308b06dc 100644 --- a/demos/demo_2d_tracer.py +++ b/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. :: diff --git a/examples/baroclinic_eddies/diagnostics.py b/examples/baroclinic_eddies/diagnostics.py index f5194e4bd..7176000d6 100644 --- a/examples/baroclinic_eddies/diagnostics.py +++ b/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) diff --git a/examples/channel_inversion/inverse_problem.py b/examples/channel_inversion/inverse_problem.py index 8a636d614..e652955ef 100644 --- a/examples/channel_inversion/inverse_problem.py +++ b/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( diff --git a/examples/discrete_turbines/channel-optimisation.py b/examples/discrete_turbines/channel-optimisation.py index 8b16efc5e..3c84c482b 100644 --- a/examples/discrete_turbines/channel-optimisation.py +++ b/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) diff --git a/examples/lockExchange/diagnostics.py b/examples/lockExchange/diagnostics.py index cfafb31c6..e3529d988 100644 --- a/examples/lockExchange/diagnostics.py +++ b/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) diff --git a/test/operations/test_diagnostics.py b/test/operations/test_diagnostics.py index cbe4701cd..6d098a59b 100644 --- a/test/operations/test_diagnostics.py +++ b/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 diff --git a/test/swe2d/test_rossby_wave.py b/test/swe2d/test_rossby_wave.py index 33a47675f..27c6e27d2 100644 --- a/test/swe2d/test_rossby_wave.py +++ b/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() diff --git a/test/swe2d/test_steady_state_channel.py b/test/swe2d/test_steady_state_channel.py index 29daca10f..e69936087 100644 --- a/test/swe2d/test_steady_state_channel.py +++ b/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) diff --git a/test/tracerEq/test_bcs_2d.py b/test/tracerEq/test_bcs_2d.py index 5052866fe..904386144 100644 --- a/test/tracerEq/test_bcs_2d.py +++ b/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... diff --git a/thetis/inversion_tools.py b/thetis/inversion_tools.py index 36d762d0a..86acdf11e 100644 --- a/thetis/inversion_tools.py +++ b/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): diff --git a/thetis/utility.py b/thetis/utility.py index f787f8f5d..3a341f06b 100755 --- a/thetis/utility.py +++ b/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