diff --git a/gusto/diagnostics/diagnostics.py b/gusto/diagnostics/diagnostics.py index 00272edac..982c2a5e9 100644 --- a/gusto/diagnostics/diagnostics.py +++ b/gusto/diagnostics/diagnostics.py @@ -7,7 +7,7 @@ ds_b, ds_v, ds_t, dS_h, dS_v, ds, dS, div, avg, pi, TensorFunctionSpace, SpatialCoordinate, as_vector, Projector, assemble, FunctionSpace, FiniteElement, - TensorProductElement) + TensorProductElement, CellVolume, Cofunction) from firedrake.assign import Assigner from firedrake.__future__ import interpolate from ufl.domain import extract_unique_domain @@ -410,12 +410,8 @@ def setup(self, domain, state_fields): V = FunctionSpace(domain.mesh, "DG", 0) test = TestFunction(V) - cell_volume = Function(V) - self.cell_flux = Function(V) - - # Calculate cell volumes - One = Function(V).assign(1) - assemble(One*test*dx, tensor=cell_volume) + cell_volume = Function(V).interpolate(CellVolume(domain.mesh)) + self.cell_flux = Cofunction(V.dual()) # Get the velocity that is being used if type(self.velocity) is str: @@ -450,7 +446,10 @@ def setup(self, domain, state_fields): self.cell_flux_form = 2*avg(un*test)*dS_calc + un*test*ds_calc # Final Courant number expression - self.expr = self.cell_flux * domain.dt / cell_volume + cell_flux = self.cell_flux.riesz_representation( + 'l2', solver_options={'function_space': V} + ) + self.expr = cell_flux * domain.dt / cell_volume super().setup(domain, state_fields) diff --git a/gusto/spatial_methods/augmentation.py b/gusto/spatial_methods/augmentation.py index 42aedaa64..0a893394d 100644 --- a/gusto/spatial_methods/augmentation.py +++ b/gusto/spatial_methods/augmentation.py @@ -84,12 +84,14 @@ def __init__( F, Z = split(self.X) test_F, test_Z = self.tests + quad = domain.max_quad_degree + if hasattr(domain.mesh, "_base_mesh"): - self.ds = ds_b + ds_t + ds_v - self.dS = dS_v + dS_h + self.ds = ds_b(degree=quad) + ds_t(degree=quad) + ds_v(degree=quad) + self.dS = dS_v(degree=quad) + dS_h(degree=quad) else: - self.ds = ds - self.dS = dS + self.ds = ds(degree=quad) + self.dS = dS(degree=quad) # Add boundary conditions self.bcs = [] diff --git a/gusto/spatial_methods/diffusion_methods.py b/gusto/spatial_methods/diffusion_methods.py index c36fe1be3..768b7532e 100644 --- a/gusto/spatial_methods/diffusion_methods.py +++ b/gusto/spatial_methods/diffusion_methods.py @@ -44,7 +44,8 @@ def interior_penalty_diffusion_form(domain, test, q, parameters): :class:`ufl.Form`: the diffusion form. """ - dS_ = (dS_v + dS_h) if domain.mesh.extruded else dS + quad = domain.max_quad_degree + dS_ = (dS_v(degree=quad) + dS_h(degree=quad)) if domain.mesh.extruded else dS(degree=quad) kappa = parameters.kappa mu = parameters.mu diff --git a/gusto/spatial_methods/transport_methods.py b/gusto/spatial_methods/transport_methods.py index f075eeadd..8d194796f 100644 --- a/gusto/spatial_methods/transport_methods.py +++ b/gusto/spatial_methods/transport_methods.py @@ -310,7 +310,8 @@ def upwind_advection_form(domain, test, q, ibp=IntegrateByParts.ONCE, outflow=Fa if outflow and ibp == IntegrateByParts.NEVER: raise ValueError("outflow is True and ibp is None are incompatible options") Vu = domain.spaces("HDiv") - dS_ = (dS_v + dS_h) if Vu.extruded else dS + quad = domain.max_quad_degree + dS_ = (dS_v(degree=quad) + dS_h(degree=quad)) if Vu.extruded else dS(degree=quad) ubar = Function(Vu) if ibp == IntegrateByParts.ONCE: @@ -372,10 +373,12 @@ def upwind_advection_form_1d(domain, test, q, ibp=IntegrateByParts.ONCE, ubar = Function(Vu) n = FacetNormal(domain.mesh) un = 0.5*(ubar * n[0] + abs(ubar * n[0])) + quad = domain.max_quad_degree + dS_ = dS(degree=quad) if ibp == IntegrateByParts.ONCE: L = -(test * ubar).dx(0) * q * dx + \ - jump(test) * (un('+')*q('+') - un('-')*q('-'))*dS + jump(test) * (un('+')*q('+') - un('-')*q('-'))*dS_ else: raise NotImplementedError("1d advection form only implemented with option ibp=IntegrateByParts.ONCE") @@ -414,7 +417,8 @@ def upwind_continuity_form(domain, test, q, ibp=IntegrateByParts.ONCE, outflow=F if outflow and ibp == IntegrateByParts.NEVER: raise ValueError("outflow is True and ibp is None are incompatible options") Vu = domain.spaces("HDiv") - dS_ = (dS_v + dS_h) if Vu.extruded else dS + quad = domain.max_quad_degree + dS_ = (dS_v(degree=quad) + dS_h(degree=quad)) if Vu.extruded else dS(degree=quad) ubar = Function(Vu) if ibp == IntegrateByParts.ONCE: @@ -476,10 +480,12 @@ def upwind_continuity_form_1d(domain, test, q, ibp=IntegrateByParts.ONCE, ubar = Function(Vu) n = FacetNormal(domain.mesh) un = 0.5*(ubar * n[0] + abs(ubar * n[0])) + quad = domain.max_quad_degree + dS_ = dS(degree=quad) if ibp == IntegrateByParts.ONCE: L = -test.dx(0) * q * ubar * dx \ - + jump(test) * (un('+')*q('+') - un('-')*q('-')) * dS + + jump(test) * (un('+')*q('+') - un('-')*q('-')) * dS_ else: raise NotImplementedError("1d continuity form only implemented with option ibp=IntegrateByParts.ONCE") @@ -520,7 +526,8 @@ def upwind_tracer_conservative_form(domain, test, q, rho, if outflow and ibp == IntegrateByParts.NEVER: raise ValueError("outflow is True and ibp is None are incompatible options") Vu = domain.spaces("HDiv") - dS_ = (dS_v + dS_h) if Vu.extruded else dS + quad = domain.max_quad_degree + dS_ = (dS_v(degree=quad) + dS_h(degree=quad)) if Vu.extruded else dS(degree=quad) ubar = Function(Vu) if ibp == IntegrateByParts.ONCE: @@ -576,7 +583,8 @@ def vector_manifold_advection_form(domain, test, q, ibp=IntegrateByParts.ONCE, o # TODO: there should maybe be a restriction on IBP here Vu = domain.spaces("HDiv") - dS_ = (dS_v + dS_h) if Vu.extruded else dS + quad = domain.max_quad_degree + dS_ = (dS_v(degree=quad) + dS_h(degree=quad)) if Vu.extruded else dS(degree=quad) ubar = Function(Vu) n = FacetNormal(domain.mesh) un = 0.5*(dot(ubar, n) + abs(dot(ubar, n))) @@ -613,7 +621,8 @@ def vector_manifold_continuity_form(domain, test, q, ibp=IntegrateByParts.ONCE, L = upwind_continuity_form(domain, test, q, ibp, outflow) Vu = domain.spaces("HDiv") - dS_ = (dS_v + dS_h) if Vu.extruded else dS + quad = domain.max_quad_degree + dS_ = (dS_v(degree=quad) + dS_h(degree=quad)) if Vu.extruded else dS(degree=quad) ubar = Function(Vu) n = FacetNormal(domain.mesh) un = 0.5*(dot(ubar, n) + abs(dot(ubar, n))) @@ -653,7 +662,8 @@ def upwind_circulation_form(domain, test, q, ibp=IntegrateByParts.ONCE): """ Vu = domain.spaces("HDiv") - dS_ = (dS_v + dS_h) if Vu.extruded else dS + quad = domain.max_quad_degree + dS_ = (dS_v(degree=quad) + dS_h(degree=quad)) if Vu.extruded else dS(degree=quad) ubar = Function(Vu) n = FacetNormal(domain.mesh) Upwind = 0.5*(sign(dot(ubar, n))+1) diff --git a/integration-tests/data/boussinesq_compressible_chkpt.h5 b/integration-tests/data/boussinesq_compressible_chkpt.h5 index 2ad2cab9d..13819691c 100644 Binary files a/integration-tests/data/boussinesq_compressible_chkpt.h5 and b/integration-tests/data/boussinesq_compressible_chkpt.h5 differ diff --git a/integration-tests/data/boussinesq_incompressible_chkpt.h5 b/integration-tests/data/boussinesq_incompressible_chkpt.h5 index c2aadff89..a4612c676 100644 Binary files a/integration-tests/data/boussinesq_incompressible_chkpt.h5 and b/integration-tests/data/boussinesq_incompressible_chkpt.h5 differ diff --git a/integration-tests/data/dry_compressible_chkpt.h5 b/integration-tests/data/dry_compressible_chkpt.h5 index 523aef482..c10668ea1 100644 Binary files a/integration-tests/data/dry_compressible_chkpt.h5 and b/integration-tests/data/dry_compressible_chkpt.h5 differ diff --git a/integration-tests/data/linear_sw_wave_chkpt.h5 b/integration-tests/data/linear_sw_wave_chkpt.h5 index 313448078..70c7c0244 100644 Binary files a/integration-tests/data/linear_sw_wave_chkpt.h5 and b/integration-tests/data/linear_sw_wave_chkpt.h5 differ diff --git a/integration-tests/data/moist_compressible_chkpt.h5 b/integration-tests/data/moist_compressible_chkpt.h5 index bd43cdadc..60381d2b6 100644 Binary files a/integration-tests/data/moist_compressible_chkpt.h5 and b/integration-tests/data/moist_compressible_chkpt.h5 differ diff --git a/integration-tests/data/simult_SIQN_order0_chkpt.h5 b/integration-tests/data/simult_SIQN_order0_chkpt.h5 index ae942022f..32e48d376 100644 Binary files a/integration-tests/data/simult_SIQN_order0_chkpt.h5 and b/integration-tests/data/simult_SIQN_order0_chkpt.h5 differ diff --git a/integration-tests/data/simult_SIQN_order1_chkpt.h5 b/integration-tests/data/simult_SIQN_order1_chkpt.h5 index d8210cdd3..f13913a48 100644 Binary files a/integration-tests/data/simult_SIQN_order1_chkpt.h5 and b/integration-tests/data/simult_SIQN_order1_chkpt.h5 differ diff --git a/integration-tests/data/sw_fplane_chkpt.h5 b/integration-tests/data/sw_fplane_chkpt.h5 index 03cde643b..d340df20d 100644 Binary files a/integration-tests/data/sw_fplane_chkpt.h5 and b/integration-tests/data/sw_fplane_chkpt.h5 differ diff --git a/pyproject.toml b/pyproject.toml index cbe693e26..031d28e95 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -24,19 +24,5 @@ classifiers = [ Homepage = "http://www.firedrakeproject.org/gusto/" Repository = "https://github.com/firedrakeproject/gusto.git" -[tool.setuptools] -packages = [ - "gusto", - "gusto.complex_proxy", - "gusto.core", - "gusto.diagnostics", - "gusto.equations", - "gusto.initialisation", - "gusto.physics", - "gusto.recovery", - "gusto.rexi", - "gusto.solvers", - "gusto.spatial_methods", - "gusto.time_discretisation", - "gusto.timestepping", -] +[tool.setuptools.packages.find] +include = ["gusto"]