Merge branch 'main' into nonsplit_physics

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)
 

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 = []

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
 

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)

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"]