Merge pull request #65 from pyroteus/61_fix_adjoint

Catch up with Firedrake updates

.github/workflows/build.yml

@@ -27,13 +27,6 @@ jobs:
         uses: actions/setup-python@v2
         with:
           python-version: 3.8
-      - name: Install Animate
-        run: |
-          . /home/firedrake/firedrake/bin/activate
-          cd ..
-          git clone https://github.com/pyroteus/animate.git
-          cd animate
-          python -m pip install -e .
       - name: Install Goalie
         run: |
           . /home/firedrake/firedrake/bin/activate

demos/burgers-hessian.py

@@ -27,10 +27,11 @@ def get_form(mesh_seq):
     def form(index, solutions):
         u, u_ = solutions["u"]
         P = mesh_seq.time_partition
-        dt = Constant(P.timesteps[index])
 
-        # Specify viscosity coefficient
-        nu = Constant(0.0001)
+        # Define constants
+        R = FunctionSpace(mesh_seq[index], "R", 0)
+        dt = Function(R).assign(P.timesteps[index])
+        nu = Function(R).assign(0.0001)
 
         # Setup variational problem
         v = TestFunction(u.function_space())

demos/burgers.py

@@ -50,17 +50,19 @@ def get_function_spaces(mesh):
 #
 # Timestepping information associated with a given subinterval
 # can be accessed via the :attr:`TimePartition` attribute of
-# the :class:`MeshSeq`. ::
+# the :class:`MeshSeq`. For technical reasons, we need to create a :class:`Function`
+# in the `'R'` space (of real numbers) to hold constants. ::
 
 
 def get_form(mesh_seq):
     def form(index, solutions):
         u, u_ = solutions["u"]
         P = mesh_seq.time_partition
-        dt = Constant(P.timesteps[index])
 
-        # Specify viscosity coefficient
-        nu = Constant(0.0001)
+        # Define constants
+        R = FunctionSpace(mesh_seq[index], "R", 0)
+        dt = Function(R).assign(P.timesteps[index])
+        nu = Function(R).assign(0.0001)
 
         # Setup variational problem
         v = TestFunction(u.function_space())

demos/burgers1.py

@@ -28,10 +28,11 @@ def get_form(mesh_seq):
     def form(index, solutions):
         u, u_ = solutions["u"]
         P = mesh_seq.time_partition
-        dt = Constant(P.timesteps[index])
 
-        # Specify viscosity coefficient
-        nu = Constant(0.0001)
+        # Define constants
+        R = FunctionSpace(mesh_seq[index], "R", 0)
+        dt = Function(R).assign(P.timesteps[index])
+        nu = Function(R).assign(0.0001)
 
         # Setup variational problem
         v = TestFunction(u.function_space())

demos/burgers2.py

@@ -28,10 +28,11 @@ def get_form(mesh_seq):
     def form(index, solutions):
         u, u_ = solutions["u"]
         P = mesh_seq.time_partition
-        dt = Constant(P.timesteps[index])
 
-        # Specify viscosity coefficient
-        nu = Constant(0.0001)
+        # Define constants
+        R = FunctionSpace(mesh_seq[index], "R", 0)
+        dt = Function(R).assign(P.timesteps[index])
+        nu = Function(R).assign(0.0001)
 
         # Setup variational problem
         v = TestFunction(u.function_space())

demos/burgers_ee.py

@@ -43,10 +43,11 @@ def get_form(mesh_seq):
     def form(index, solutions):
         u, u_ = solutions["u"]
         P = mesh_seq.time_partition
-        dt = Constant(P.timesteps[index])
 
-        # Specify viscosity coefficient
-        nu = Constant(0.0001)
+        # Define constants
+        R = FunctionSpace(mesh_seq[index], "R", 0)
+        dt = Function(R).assign(P.timesteps[index])
+        nu = Function(R).assign(0.0001)
 
         # Setup variational problem
         v = TestFunction(u.function_space())

demos/burgers_oo.py

@@ -33,10 +33,11 @@ def get_form(self):
         def form(index, solutions):
             u, u_ = solutions["u"]
             P = self.time_partition
-            dt = Constant(P.timesteps[index])
 
-            # Specify viscosity coefficient
-            nu = Constant(0.0001)
+            # Define constants
+            R = FunctionSpace(mesh_seq[index], "R", 0)
+            dt = Function(R).assign(P.timesteps[index])
+            nu = Function(R).assign(0.0001)
 
             # Setup variational problem
             v = TestFunction(u.function_space())
@@ -84,7 +85,8 @@ def get_initial_condition(self):
 
     @annotate_qoi
     def get_qoi(self, solutions, i):
-        dt = Constant(self.time_partition[i].timestep)
+        R = FunctionSpace(self[i], "R", 0)
+        dt = Function(R).assign(self.time_partition[i].timestep)
 
         def end_time_qoi():
             u = solutions["u"]

demos/burgers_time_integrated.py

@@ -23,10 +23,11 @@ def get_form(mesh_seq):
     def form(index, solutions):
         u, u_ = solutions["u"]
         P = mesh_seq.time_partition
-        dt = Constant(P.timesteps[index])
 
-        # Specify viscosity coefficient
-        nu = Constant(0.0001)
+        # Define constants
+        R = FunctionSpace(mesh_seq[index], "R", 0)
+        dt = Function(R).assign(P.timesteps[index])
+        nu = Function(R).assign(0.0001)
 
         # Setup variational problem
         v = TestFunction(u.function_space())
@@ -93,12 +94,13 @@ def solver(index, ic):
 #    \;\mathrm dy\;\mathrm dt.
 #
 # Note that in this case we multiply by the timestep.
-# It is wrapped in a :class:`Constant` to avoid
+# It is wrapped in a :class:`Function` from `'R'` space to avoid
 # recompilation if the value is changed. ::
 
 
 def get_qoi(mesh_seq, solutions, i):
-    dt = Constant(mesh_seq.time_partition[i].timestep)
+    R = FunctionSpace(mesh_seq[i], "R", 0)
+    dt = Function(R).assign(mesh_seq.time_partition[i].timestep)
 
     def time_integrated_qoi(t):
         u = solutions["u"]

demos/gray_scott.py

@@ -57,11 +57,12 @@ def form(index, sols):
         psi_a, psi_b = TestFunctions(mesh_seq.function_spaces["ab"][index])
 
         # Define constants
-        dt = Constant(mesh_seq.time_partition[index].timestep)
-        D_a = Constant(8.0e-05)
-        D_b = Constant(4.0e-05)
-        gamma = Constant(0.024)
-        kappa = Constant(0.06)
+        R = FunctionSpace(mesh_seq[index], "R", 0)
+        dt = Function(R).assign(mesh_seq.time_partition[index].timestep)
+        D_a = Function(R).assign(8.0e-05)
+        D_b = Function(R).assign(4.0e-05)
+        gamma = Function(R).assign(0.024)
+        kappa = Function(R).assign(0.06)
 
         # Write the two equations in variational form
         F = (

demos/gray_scott_split.py

@@ -57,11 +57,12 @@ def form(index, sols):
         psi_b = TestFunction(mesh_seq.function_spaces["b"][index])
 
         # Define constants
-        dt = Constant(mesh_seq.time_partition[index].timestep)
-        D_a = Constant(8.0e-05)
-        D_b = Constant(4.0e-05)
-        gamma = Constant(0.024)
-        kappa = Constant(0.06)
+        R = FunctionSpace(mesh_seq[index], "R", 0)
+        dt = Function(R).assign(mesh_seq.time_partition[index].timestep)
+        D_a = Function(R).assign(8.0e-05)
+        D_b = Function(R).assign(4.0e-05)
+        gamma = Function(R).assign(0.024)
+        kappa = Function(R).assign(0.06)
 
         # Write the two equations in variational form
         F_a = (

demos/point_discharge2d-goal_oriented.py

@@ -36,11 +36,16 @@ def get_form(mesh_seq):
     def form(index, sols):
         c, c_ = sols["c"]
         function_space = mesh_seq.function_spaces["c"][index]
-        D = Constant(0.1)
-        u = Constant(as_vector([1, 0]))
         h = CellSize(mesh_seq[index])
         S = source(mesh_seq[index])
 
+        # Define constants
+        R = FunctionSpace(mesh_seq[index], "R", 0)
+        D = Function(R).assign(0.1)
+        u_x = Function(R).assign(1.0)
+        u_y = Function(R).assign(0.0)
+        u = as_vector([u_x, u_y])
+
         # SUPG stabilisation parameter
         unorm = sqrt(dot(u, u))
         tau = 0.5 * h / unorm

demos/point_discharge2d-hessian.py

@@ -40,11 +40,16 @@ def get_form(mesh_seq):
     def form(index, sols):
         c, c_ = sols["c"]
         function_space = mesh_seq.function_spaces["c"][index]
-        D = Constant(0.1)
-        u = Constant(as_vector([1, 0]))
         h = CellSize(mesh_seq[index])
         S = source(mesh_seq[index])
 
+        # Define constants
+        R = FunctionSpace(mesh_seq[index], "R", 0)
+        D = Function(R).assign(0.1)
+        u_x = Function(R).assign(1.0)
+        u_y = Function(R).assign(0.0)
+        u = as_vector([u_x, u_y])
+
         # SUPG stabilisation parameter
         unorm = sqrt(dot(u, u))
         tau = 0.5 * h / unorm

demos/point_discharge2d.py

@@ -71,11 +71,16 @@ def get_form(mesh_seq):
     def form(index, sols):
         c, c_ = sols["c"]
         function_space = mesh_seq.function_spaces["c"][index]
-        D = Constant(0.1)
-        u = Constant(as_vector([1, 0]))
         h = CellSize(mesh_seq[index])
         S = source(mesh_seq[index])
 
+        # Define constants
+        R = FunctionSpace(mesh_seq[index], "R", 0)
+        D = Function(R).assign(0.1)
+        u_x = Function(R).assign(1.0)
+        u_y = Function(R).assign(0.0)
+        u = as_vector([u_x, u_y])
+
         # SUPG stabilisation parameter
         unorm = sqrt(dot(u, u))
         tau = 0.5 * h / unorm

demos/solid_body_rotation.py

@@ -151,10 +151,14 @@ def form(index, sols, field="c"):
         V = mesh_seq.function_spaces[field][index]
         mesh = mesh_seq[index]
 
+        # Define velocity field
         x, y = SpatialCoordinate(mesh)
         u = as_vector([-y, x])
-        dt = Constant(mesh_seq.time_partition[index].timestep)
-        theta = Constant(0.5)
+
+        # Define constants
+        R = FunctionSpace(mesh_seq[index], "R", 0)
+        dt = Function(R).assign(mesh_seq.time_partition[index].timestep)
+        theta = Function(R).assign(0.5)
 
         psi = TrialFunction(V)
         phi = TestFunction(V)