diff --git a/.DS_Store b/.DS_Store
new file mode 100644
index 0000000..5f7d455
Binary files /dev/null and b/.DS_Store differ
diff --git a/examples/.DS_Store b/examples/.DS_Store
new file mode 100644
index 0000000..0140f22
Binary files /dev/null and b/examples/.DS_Store differ
diff --git a/examples/burgers/meshgen.py b/examples/burgers/meshgen.py
deleted file mode 100644
index 3467cea..0000000
--- a/examples/burgers/meshgen.py
+++ /dev/null
@@ -1,2 +0,0 @@
-def generate_geo(config, reverse=False):
-    return
diff --git a/examples/compute_importance.py b/examples/compute_importance.py
index 3e2bf4c..2751ad4 100644
--- a/examples/compute_importance.py
+++ b/examples/compute_importance.py
@@ -20,7 +20,7 @@
     "model",
     help="The model",
     type=str,
-    choices=["steady_turbine"],
+    choices=["steady_turbine", "pyroteus_burgers"],
 )
 parser.add_argument(
     "num_training_cases",
diff --git a/examples/makefile b/examples/makefile
index 259ce67..1820ad1 100644
--- a/examples/makefile
+++ b/examples/makefile
@@ -3,7 +3,7 @@ all: setup network test
 # --- Configurable parameters
 
 APPROACHES			= anisotropic
-MODEL				= steady_turbine
+MODEL				= pyroteus_burgers
 NUM_TRAINING_CASES	= 100
 TESTING_CASES		= $(shell cat $(MODEL)/testing_cases.txt)
 PETSC_OPTIONS		= -dm_plex_metric_hausdorff_number 1
diff --git a/examples/models/.DS_Store b/examples/models/.DS_Store
new file mode 100644
index 0000000..5008ddf
Binary files /dev/null and b/examples/models/.DS_Store differ
diff --git a/examples/models/burgers.py b/examples/models/burgers.py
deleted file mode 100644
index f6fdea0..0000000
--- a/examples/models/burgers.py
+++ /dev/null
@@ -1,165 +0,0 @@
-from firedrake import *
-from firedrake.petsc import PETSc
-import nn_adapt.model
-import nn_adapt.solving
-
-
-class Parameters(nn_adapt.model.Parameters):
-    """
-    Class encapsulating all parameters required for a simple
-    Burgers equation test case.
-    """
-
-    qoi_name = "right boundary integral"
-    qoi_unit = r"m\,s^{-1}"
-
-    # Adaptation parameters
-    h_min = 1.0e-10  # Minimum metric magnitude
-    h_max = 1.0  # Maximum metric magnitude
-
-    # Physical parameters
-    viscosity_coefficient = 0.0001
-    initial_speed = 1.0
-
-    # Timestepping parameters
-    timestep = 0.05
-
-    solver_parameters = {}
-    adjoint_solver_parameters = {}
-
-    def bathymetry(self, mesh):
-        """
-        Compute the bathymetry field on the current `mesh`.
-
-        Note that there isn't really a concept of bathymetry
-        for Burgers equation. It is kept constant and should
-        be ignored by the network.
-        """
-        P0_2d = FunctionSpace(mesh, "DG", 0)
-        return Function(P0_2d).assign(1.0)
-
-    def drag(self, mesh):
-        """
-        Compute the bathymetry field on the current `mesh`.
-
-        Note that there isn't really a concept of bathymetry
-        for Burgers equation. It is kept constant and should
-        be ignored by the network.
-        """
-        P0_2d = FunctionSpace(mesh, "DG", 0)
-        return Function(P0_2d).assign(1.0)
-
-    def viscosity(self, mesh):
-        """
-        Compute the viscosity coefficient on the current `mesh`.
-        """
-        P0_2d = FunctionSpace(mesh, "DG", 0)
-        return Function(P0_2d).assign(self.viscosity_coefficient)
-
-    def ic(self, mesh):
-        """
-        Initial condition
-        """
-        x, y = SpatialCoordinate(mesh)
-        expr = self.initial_speed * sin(pi * x)
-        return as_vector([expr, 0])
-
-
-PETSc.Sys.popErrorHandler()
-parameters = Parameters()
-
-
-def get_function_space(mesh):
-    r"""
-    Construct the :math:`\mathbb P2` finite element space
-    used for the prognostic solution.
-    """
-    return VectorFunctionSpace(mesh, "CG", 2)
-
-
-class Solver(nn_adapt.solving.Solver):
-    """
-    Solver object based on current mesh and state.
-    """
-
-    def __init__(self, mesh, ic, **kwargs):
-        """
-        :arg mesh: the mesh to define the solver on
-        :arg ic: the current state / initial condition
-        """
-        self.mesh = mesh
-
-        # Collect parameters
-        dt = Constant(parameters.timestep)
-        nu = parameters.viscosity(mesh)
-
-        # Define variational formulation
-        V = self.function_space
-        u = Function(V)
-        u_ = Function(V)
-        v = TestFunction(V)
-        self._form = (
-            inner((u - u_) / dt, v) * dx
-            + inner(dot(u, nabla_grad(u)), v) * dx
-            + nu * inner(grad(u), grad(v)) * dx
-        )
-        problem = NonlinearVariationalProblem(self._form, u)
-
-        # Set initial condition
-        u_.project(parameters.ic(mesh))
-
-        # Create solver
-        self._solver = NonlinearVariationalSolver(problem)
-        self._solution = u
-
-    @property
-    def function_space(self):
-        r"""
-        The :math:`\mathbb P2` finite element space.
-        """
-        return get_function_space(self.mesh)
-
-    @property
-    def form(self):
-        """
-        The weak form of Burgers equation
-        """
-        return self._form
-
-    @property
-    def solution(self):
-        return self._solution
-
-    def iterate(self, **kwargs):
-        """
-        Take a single timestep of Burgers equation
-        """
-        self._solver.solve()
-
-
-def get_initial_condition(function_space):
-    """
-    Compute an initial condition based on the initial
-    speed parameter.
-    """
-    u = Function(function_space)
-    u.interpolate(parameters.ic(function_space.mesh()))
-    return u
-
-
-def get_qoi(mesh):
-    """
-    Extract the quantity of interest function from the :class:`Parameters`
-    object.
-
-    It should have one argument - the prognostic solution.
-    """
-
-    def qoi(sol):
-        return inner(sol, sol) * ds(2)
-
-    return qoi
-
-
-# Initial mesh for all test cases
-initial_mesh = UnitSquareMesh(30, 30)
diff --git a/examples/models/pyroteus_burgers.py b/examples/models/pyroteus_burgers.py
new file mode 100644
index 0000000..45f66c0
--- /dev/null
+++ b/examples/models/pyroteus_burgers.py
@@ -0,0 +1,196 @@
+from firedrake import *
+from pyroteus import *
+import pyroteus.go_mesh_seq
+from firedrake.petsc import PETSc
+import nn_adapt.model
+
+'''
+A memory hungry method solving time dependent PDE.
+'''
+class Parameters(nn_adapt.model.Parameters):
+    """
+    Class encapsulating all parameters required for a simple
+    Burgers equation test case.
+    """
+
+    qoi_name = "right boundary integral"
+    qoi_unit = r"m\,s^{-1}"
+
+    # Adaptation parameters
+    h_min = 1.0e-10  # Minimum metric magnitude
+    h_max = 1.0  # Maximum metric magnitude
+
+    # Physical parameters
+    viscosity_coefficient = 0.0001
+    initial_speed = 1.0
+    
+    # Offset for creating more initial conditions
+    x_offset = 0
+    y_offset = 0
+
+    # Timestepping parameters
+    timestep = 5
+
+    solver_parameters = {}
+    adjoint_solver_parameters = {}
+
+    def bathymetry(self, mesh):
+        """
+        Compute the bathymetry field on the current `mesh`.
+        Note that there isn't really a concept of bathymetry
+        for Burgers equation. It is kept constant and should
+        be ignored by the network.
+        """
+        P0_2d = FunctionSpace(mesh, "DG", 0)
+        return Function(P0_2d).assign(1.0)
+
+    def drag(self, mesh):
+        """
+        Compute the bathymetry field on the current `mesh`.
+        Note that there isn't really a concept of bathymetry
+        for Burgers equation. It is kept constant and should
+        be ignored by the network.
+        """
+        P0_2d = FunctionSpace(mesh, "DG", 0)
+        return Function(P0_2d).assign(1.0)
+
+    def viscosity(self, mesh):
+        """
+        Compute the viscosity coefficient on the current `mesh`.
+        """
+        P0_2d = FunctionSpace(mesh, "DG", 0)
+        return Function(P0_2d).assign(self.viscosity_coefficient)
+
+    def ic(self, mesh):
+        """
+        Initial condition
+        """
+        x, y = SpatialCoordinate(mesh)
+        # x_expr = self.initial_speed * sin(pi * x + self.x_offset)
+        # y_expr = self.initial_speed * sin(pi * y + self.y_offset)
+        x_expr = self.initial_speed * sin(pi * x)
+        y_expr = 0
+        return as_vector([x_expr, y_expr])
+    
+    
+def get_function_spaces(mesh):
+    return {"u": VectorFunctionSpace(mesh, "CG", 2)}
+
+
+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 = parameters.viscosity_coefficient
+
+        # Setup variational problem
+        v = TestFunction(u.function_space())
+        F = (
+            inner((u - u_) / dt, v) * dx
+            + inner(dot(u, nabla_grad(u)), v) * dx
+            + nu * inner(grad(u), grad(v)) * dx
+        )
+        return F
+
+    return form
+
+
+def get_solver(mesh_seq):
+    def solver(index, ic):
+        function_space = mesh_seq.function_spaces["u"][index]
+        u = Function(function_space)
+
+        # Initialise 'lagged' solution
+        u_ = Function(function_space, name="u_old")
+        u_.assign(ic["u"])
+
+        # Define form
+        F = mesh_seq.form(index, {"u": (u, u_)})
+
+        # Time integrate from t_start to t_end
+        P = mesh_seq.time_partition
+        t_start, t_end = P.subintervals[index]
+        dt = P.timesteps[index]
+        t = t_start
+        step = 0
+        sp = {'snes_max_it': 100}
+        while t < t_end - 1.0e-05:
+            step += 1
+            solve(F == 0, u, ad_block_tag="u", solver_parameters=sp)
+            u_.assign(u)
+            t += dt
+        return {"u": u}
+
+    return solver
+
+
+def get_initial_condition(mesh_seq):
+    fs = mesh_seq.function_spaces["u"][0]
+    return {"u": interpolate(parameters.ic(fs), fs)}
+
+def get_qoi(mesh_seq, solutions, index):
+    def end_time_qoi():
+        u = solutions["u"]
+        fs = mesh_seq.function_spaces["u"][0]
+        x, y = SpatialCoordinate(fs)
+        partial = conditional(And(And(x > 0.1, x < 0.9), And(y > 0.45, y < 0.55)), 1, 0)
+        return inner(u, u) * dx
+
+    def time_integrated_qoi(t):
+        fs = mesh_seq.function_spaces["u"][0]
+        x, y = SpatialCoordinate(fs)
+        partial = conditional(And(And(x > 0.7, x < 0.8), And(y > 0.45, y < 0.55)), 1, 0)
+        dt = Constant(mesh_seq.time_partition[index].timestep)
+        u = solutions["u"]
+        return dt * partial * inner(u, u) * dx
+    
+    if mesh_seq.qoi_type == "end_time":
+        return end_time_qoi
+    else:
+        return time_integrated_qoi
+
+
+PETSc.Sys.popErrorHandler()
+parameters = Parameters()
+
+def GoalOrientedMeshSeq(mesh, **kwargs):
+    fields = ["u"]
+    
+    try:
+        num_subintervals = len(mesh)
+    except:
+        num_subintervals = 1
+
+    # setup time steps and export steps
+    dt = 1
+    steps_subintervals = 10
+    end_time = num_subintervals * steps_subintervals * dt
+    timesteps_per_export = 1
+    
+    # setup pyroteus time_partition
+    time_partition = TimePartition(
+        end_time,
+        num_subintervals,
+        dt,
+        fields,
+        timesteps_per_export=timesteps_per_export,
+        )
+    
+    mesh_seq = pyroteus.go_mesh_seq.GoalOrientedMeshSeq(
+        time_partition,
+        mesh,
+        get_function_spaces=get_function_spaces,
+        get_initial_condition=get_initial_condition,
+        get_form=get_form,
+        get_solver=get_solver,
+        get_qoi=get_qoi,
+        qoi_type="end_time",
+    )
+    return mesh_seq
+    
+
+initial_mesh = lambda: UnitSquareMesh(30, 30)
+        
\ No newline at end of file
diff --git a/examples/models/steady_turbine.py b/examples/models/pyroteus_turbine.py
similarity index 55%
rename from examples/models/steady_turbine.py
rename to examples/models/pyroteus_turbine.py
index b2ade51..00260e6 100644
--- a/examples/models/steady_turbine.py
+++ b/examples/models/pyroteus_turbine.py
@@ -1,8 +1,7 @@
 from thetis import *
+from pyroteus import *
+import pyroteus.go_mesh_seq
 import nn_adapt.model
-import nn_adapt.solving
-import numpy as np
-
 
 class Parameters(nn_adapt.model.Parameters):
     """
@@ -18,6 +17,10 @@ class Parameters(nn_adapt.model.Parameters):
     # Adaptation parameters
     h_min = 1.0e-08
     h_max = 500.0
+    
+    # time dependent parameters
+    end_time = 100.
+    time_steps = 10.
 
     # Physical parameters
     viscosity_coefficient = 0.5
@@ -206,123 +209,145 @@ def viscosity(self, mesh):
 
 PETSc.Sys.popErrorHandler()
 parameters = Parameters()
+    
+def GoalOrientedMeshSeq(mesh):
+    
+    fields = ["q"]
+    time_partition = TimeInterval(parameters.end_time, 
+                                        parameters.time_steps, 
+                                        fields, timesteps_per_export=1)
+
+
+    def get_solver(mesh_seq):
+
+        def solver(index, ic):
+            V = mesh_seq.function_spaces["q"][index]
+            q = ic["q"]
+            mesh_seq.form(index, {"q": (q, q)})
+            u_init, eta_init = q.split()
+            mesh_seq._thetis_solver.assign_initial_conditions(uv=u_init, elev=eta_init)
+            mesh_seq._thetis_solver.iterate()
+            
+            return {"q": mesh_seq._thetis_solver.fields.solution_2d}
+
+        return solver
+
+
+    def get_form(mesh_seq):
+        def form(index, ic):
+            P = mesh_seq.time_partition
+            
+            bathymetry = parameters.bathymetry(mesh_seq[index])
+            Cd = parameters.drag_coefficient
+
+            # Create solver object
+            mesh_seq._thetis_solver = solver2d.FlowSolver2d(mesh_seq[index], bathymetry)
+            options = mesh_seq._thetis_solver.options
+            options.element_family = "dg-cg"
+            options.timestep = P.timestep
+            options.simulation_export_time = P.timestep * P.timesteps_per_export[index]
+            options.simulation_end_time = P.end_time
+            options.swe_timestepper_type = "SteadyState"
+            options.swe_timestepper_options.solver_parameters = parameters.solver_parameters
+            options.use_grad_div_viscosity_term = False
+            options.horizontal_viscosity = parameters.viscosity(mesh_seq[index])
+            options.quadratic_drag_coefficient = Cd
+            options.use_lax_friedrichs_velocity = True
+            options.lax_friedrichs_velocity_scaling_factor = Constant(1.0)
+            options.use_grad_depth_viscosity_term = False
+            options.no_exports = True
+
+            # Apply boundary conditions
+            mesh_seq._thetis_solver.create_function_spaces()
+            P1v_2d = mesh_seq._thetis_solver.function_spaces.P1v_2d
+            u_inflow = interpolate(parameters.u_inflow(mesh_seq[index]), P1v_2d)
+            mesh_seq._thetis_solver.bnd_functions["shallow_water"] = {
+                1: {"uv": u_inflow},  # inflow
+                2: {"elev": Constant(0.0)},  # outflow
+                3: {"un": Constant(0.0)},  # free-slip
+                4: {"uv": Constant(as_vector([0.0, 0.0]))},  # no-slip
+                5: {"elev": Constant(0.0), "un": Constant(0.0)}  # weakly reflective
+            }
+
+            # # Create tidal farm
+            options.tidal_turbine_farms = parameters.farm(mesh_seq[index])
+            mesh_seq._thetis_solver.create_timestepper()
+            
+            return mesh_seq._thetis_solver.timestepper.F
+
+        return form
+
+
+    def get_function_space(mesh):
+        """
+        Construct the (mixed) finite element space used for the
+        prognostic solution.
+        """
+        P1v_2d = get_functionspace(mesh, "DG", 1, vector=True)
+        P2_2d = get_functionspace(mesh, "CG", 2)
+        return {"q": P1v_2d * P2_2d}
+
+
+    def get_initial_condition(mesh_seq):
+        """
+        Compute an initial condition based on the inflow velocity
+        and zero free surface elevation.
+        """
+        V = mesh_seq.function_spaces["q"][0]
+        q = Function(V)
+        u, eta = q.split()
+        u.interpolate(parameters.ic(mesh_seq))
+        return {"q": q}
+
+
+    def get_qoi(mesh_seq, solutions, index):
+        """
+        Extract the quantity of interest function from the :class:`Parameters`
+        object.
+
+        It should have one argument - the prognostic solution.
+        """
+        dt = Constant(mesh_seq.time_partition[index].timestep)
+        rho = parameters.density
+        Ct = parameters.corrected_thrust_coefficient
+        At = parameters.swept_area
+        Cd = 0.5 * Ct * At * parameters.turbine_density(mesh_seq[index])
+        tags = parameters.turbine_ids
+        sol = solutions["q"]
+
+        def qoi():
+            u, eta = split(sol)
+            return sum([rho * Cd * pow(dot(u, u), 1.5) * dx(tag) for tag in tags])
+
+        return qoi
+    
+    
+    mesh_seq = pyroteus.go_mesh_seq.GoalOrientedMeshSeq(
+        time_partition,
+        mesh,
+        get_function_spaces=get_function_space,
+        get_initial_condition=get_initial_condition,
+        get_form=get_form,
+        get_solver=get_solver,
+        get_qoi=get_qoi,
+        qoi_type="end_time")
+    
+    return mesh_seq
+    
+        
+        
+
+# print(solutions["q"].keys())
+
+# fwd_file = File("./out/forward.pvd")
+# for i in range(len(solutions["u"]["forward"][0])):
+#     fwd_file.write(*solutions["q"]["forward"][0][i].split())
+    
+# adj_file = File("./out/adjoint.pvd")
+# for i in range(len(solutions["u"]["adjoint"][0])):
+#     adj_file.write(*solutions["q"]["adjoint"][0][i].split())
+
+# adj_next_file = File("./out/adjoint_next.pvd")
+# for i in range(len(solutions["u"]["adjoint_next"][0])):
+#     adj_next_file.write(*solutions["q"]["adjoint_next"][0][i].split())
 
-
-def get_function_space(mesh):
-    """
-    Construct the (mixed) finite element space used for the
-    prognostic solution.
-    """
-    P1v_2d = get_functionspace(mesh, "DG", 1, vector=True)
-    P2_2d = get_functionspace(mesh, "CG", 2)
-    return P1v_2d * P2_2d
-
-
-class Solver(nn_adapt.solving.Solver):
-    """
-    Set up a Thetis :class:`FlowSolver2d` object, based on
-    the current mesh and initial condition.
-    """
-
-    def __init__(self, mesh, ic, **kwargs):
-        """
-        :arg mesh: the mesh to define the solver on
-        :arg ic: the initial condition
-        """
-        bathymetry = parameters.bathymetry(mesh)
-        Cd = parameters.drag_coefficient
-        sp = kwargs.pop("solver_parameters", None)
-
-        # Create solver object
-        self._thetis_solver = solver2d.FlowSolver2d(mesh, bathymetry)
-        options = self._thetis_solver.options
-        options.element_family = "dg-cg"
-        options.timestep = 20.0
-        options.simulation_export_time = 20.0
-        options.simulation_end_time = 18.0
-        options.swe_timestepper_type = "SteadyState"
-        options.swe_timestepper_options.solver_parameters = (
-            sp or parameters.solver_parameters
-        )
-        options.use_grad_div_viscosity_term = False
-        options.horizontal_viscosity = parameters.viscosity(mesh)
-        options.quadratic_drag_coefficient = Cd
-        options.use_lax_friedrichs_velocity = True
-        options.lax_friedrichs_velocity_scaling_factor = Constant(1.0)
-        options.use_grad_depth_viscosity_term = False
-        options.no_exports = True
-        options.update(kwargs)
-        self._thetis_solver.create_equations()
-
-        # Apply boundary conditions
-        P1v_2d = self._thetis_solver.function_spaces.P1v_2d
-        u_inflow = interpolate(parameters.u_inflow(mesh), P1v_2d)
-        self._thetis_solver.bnd_functions["shallow_water"] = {
-            1: {"uv": u_inflow},  # inflow
-            2: {"elev": Constant(0.0)},  # outflow
-            3: {"un": Constant(0.0)},  # free-slip
-            4: {"uv": Constant(as_vector([0.0, 0.0]))},  # no-slip
-            5: {"elev": Constant(0.0), "un": Constant(0.0)}  # weakly reflective
-        }
-
-        # Create tidal farm
-        options.tidal_turbine_farms = parameters.farm(mesh)
-
-        # Apply initial guess
-        u_init, eta_init = ic.split()
-        self._thetis_solver.assign_initial_conditions(uv=u_init, elev=eta_init)
-
-    @property
-    def function_space(self):
-        r"""
-        The :math:`\mathbb P1_{DG}-\mathbb P2` function space.
-        """
-        return self._thetis_solver.function_spaces.V_2d
-
-    @property
-    def form(self):
-        """
-        The weak form of the shallow water equations.
-        """
-        return self._thetis_solver.timestepper.F
-
-    def iterate(self, **kwargs):
-        """
-        Solve the nonlinear shallow water equations.
-        """
-        self._thetis_solver.iterate(**kwargs)
-
-    @property
-    def solution(self):
-        return self._thetis_solver.fields.solution_2d
-
-
-def get_initial_condition(function_space):
-    """
-    Compute an initial condition based on the inflow velocity
-    and zero free surface elevation.
-    """
-    q = Function(function_space)
-    u, eta = q.split()
-    u.interpolate(parameters.ic(function_space.mesh()))
-    return q
-
-
-def get_qoi(mesh):
-    """
-    Extract the quantity of interest function from the :class:`Parameters`
-    object.
-
-    It should have one argument - the prognostic solution.
-    """
-    rho = parameters.density
-    Ct = parameters.corrected_thrust_coefficient
-    At = parameters.swept_area
-    Cd = 0.5 * Ct * At * parameters.turbine_density(mesh)
-    tags = parameters.turbine_ids
-
-    def qoi(sol):
-        u, eta = split(sol)
-        return sum([rho * Cd * pow(dot(u, u), 1.5) * dx(tag) for tag in tags])
-
-    return qoi
diff --git a/examples/plot_convergence.py b/examples/plot_convergence.py
index d8e2255..1df914f 100644
--- a/examples/plot_convergence.py
+++ b/examples/plot_convergence.py
@@ -71,6 +71,7 @@
 if len(approaches.keys()) == 0:
     print("Nothing to plot.")
     sys.exit(0)
+print(qois)
 
 # Drop first iteration because timings include compilation   # FIXME: Why?
 dofs["uniform"] = dofs["uniform"][1:]
diff --git a/examples/plot_importance.py b/examples/plot_importance.py
index 596d80f..d605068 100644
--- a/examples/plot_importance.py
+++ b/examples/plot_importance.py
@@ -19,7 +19,7 @@
     "model",
     help="The model",
     type=str,
-    choices=["steady_turbine"],
+    choices=["steady_turbine", "pyroteus_burgers"],
 )
 parser.add_argument(
     "num_training_cases",
diff --git a/examples/plot_progress.py b/examples/plot_progress.py
index 09482a3..ed8c5d2 100644
--- a/examples/plot_progress.py
+++ b/examples/plot_progress.py
@@ -18,7 +18,7 @@
     "model",
     help="The model",
     type=str,
-    choices=["steady_turbine", "burgers"],
+    choices=["steady_turbine", "pyroteus_burgers"],
 )
 parser.add_argument(
     "--tag",
diff --git a/examples/pyroteus_burgers/.DS_Store b/examples/pyroteus_burgers/.DS_Store
new file mode 100644
index 0000000..416cd27
Binary files /dev/null and b/examples/pyroteus_burgers/.DS_Store differ
diff --git a/examples/burgers/config.py b/examples/pyroteus_burgers/config.py
similarity index 58%
rename from examples/burgers/config.py
rename to examples/pyroteus_burgers/config.py
index 0bd51e3..62c7004 100644
--- a/examples/burgers/config.py
+++ b/examples/pyroteus_burgers/config.py
@@ -1,4 +1,4 @@
-from models.burgers import *
+from models.pyroteus_burgers import *
 from nn_adapt.ann import sample_uniform
 import numpy as np
 
@@ -25,12 +25,18 @@ def initialise(case, discrete=False):
         # Random initial speed from 0.01 m/s to 6 m/s
         parameters.initial_speed = sample_uniform(0.01, 6.0)
 
-        # Random viscosity from 0.00001 m^2/s to 1 m^2/s
-        parameters.viscosity_coefficient = sample_uniform(0.1, 1.0) * 10 ** np.random.randint(-3, 1)
+        # # Random viscosity from 0.00001 m^2/s to 1 m^2/s
+        # parameters.viscosity_coefficient = sample_uniform(0.1, 1.0) * 10 ** np.random.randint(-3, 1)
+        # Random viscosity from 0.001 m^2/s to 1 m^2/s
+        parameters.viscosity_coefficient = sample_uniform(0.01, 1.0) * 10 ** np.random.randint(-1, 1)
+        
+        # Random offset for initial conditions
+        parameters.x_offset = sample_uniform(0, 1/2*pi)
+        parameters.y_offset = sample_uniform(0, 1/2*pi)
         return
     elif "demo" in case:
-        parameters.viscosity_coefficient = 0.0001
-        parameters.initial_speed = 1.0
+        parameters.viscosity_coefficient = 0.001
+        parameters.initial_speed = 1
     else:
         raise ValueError(f"Test case {test_case} not recognised")
 
diff --git a/examples/steady_turbine/meshgen.py b/examples/pyroteus_burgers/meshgen.py
similarity index 100%
rename from examples/steady_turbine/meshgen.py
rename to examples/pyroteus_burgers/meshgen.py
diff --git a/examples/burgers/network.py b/examples/pyroteus_burgers/network.py
similarity index 94%
rename from examples/burgers/network.py
rename to examples/pyroteus_burgers/network.py
index b7f9907..5f437fb 100644
--- a/examples/burgers/network.py
+++ b/examples/pyroteus_burgers/network.py
@@ -16,7 +16,7 @@ class NetLayout(NetLayoutBase):
           + [boundary element?]
           + [12 forward DoFs per element]
           + [12 adjoint DoFs per element]
-          = 30
+          = 29
 
     Hidden layer:
     -------------
@@ -30,7 +30,7 @@ class NetLayout(NetLayoutBase):
     """
 
     inputs = (
-        "estimator_coarse",
+        # "estimator_coarse",
         "physics_viscosity",
         "mesh_d",
         "mesh_h1",
diff --git a/examples/pyroteus_burgers/plotting.py b/examples/pyroteus_burgers/plotting.py
new file mode 100644
index 0000000..b5b629a
--- /dev/null
+++ b/examples/pyroteus_burgers/plotting.py
@@ -0,0 +1,83 @@
+from firedrake import *
+import matplotlib
+import numpy as np
+
+
+matplotlib.rcParams["font.size"] = 12
+
+
+def plot_config(config, mesh, axes):
+    """
+    Plot a given configuration of a problem on a given
+    mesh and axes.
+    """
+    tags = config.parameters.turbine_ids
+    P0 = FunctionSpace(mesh, "DG", 0)
+    footprints = assemble(sum(TestFunction(P0) * dx(tag) for tag in tags))
+    footprints.interpolate(conditional(footprints > 0, 0, 1))
+    triplot(
+        mesh,
+        axes=axes,
+        boundary_kw={"color": "dodgerblue"},
+        interior_kw={"edgecolor": "w"},
+    )
+    tricontourf(footprints, axes=axes, cmap="Blues", levels=[0, 1])
+
+    # Bounding box
+    xmin = 0
+    xmax = 1200
+    ymin = 0
+    ymax = 500
+    eps = 5
+
+    # Adjust axes
+    W = assemble(Constant(1.0, domain=mesh) * ds(1))
+    L = 0.5 * assemble(
+        Constant(1.0, domain=mesh) * ds(3)
+    )  # NOTE: both top and bottom are tagged as 3
+    dL = 0.5 * (xmax - L)
+    dW = 0.5 * (ymax - W)
+    axes.axis(False)
+    axes.set_xlim([xmin - dL - eps, xmax - dL + eps])
+    axes.set_ylim([ymin - dW - eps, ymax - dW + eps])
+
+    # Annotate with viscosity coefficient and bathymetry
+    nu = config.parameters.viscosity_coefficient
+    b = config.parameters.depth
+    u_in = config.parameters.inflow_speed
+    txt = r"$\nu$ = %.3f, $b$ = %.2f, $u_{\mathrm{in}}$ = %.2f" % (nu, b, u_in)
+    axes.annotate(
+        txt, xy=(0.025 * L, -0.25 * W), bbox={"fc": "w"}, annotation_clip=False
+    )
+
+
+def process_sensitivities(data, layout):
+    """
+    Separate sensitivity experiment data by variable.
+
+    :arg data: the output of `compute_importance.py`
+    :arg layout: the :class:`NetLayout` instance
+    """
+    i = 0
+    sensitivities = {}
+    dofs = {"u": 3, "v": 3, r"\eta": 6}
+    for label in ("physics", "mesh", "forward", "adjoint"):
+        n = layout.count_inputs(label)
+        if n == 0:
+            continue
+        if label in ("forward", "adjoint"):
+            assert n == sum(dofs.values())
+            for key, dof in dofs.items():
+                S = np.zeros(6)
+                for j in range(dof):
+                    S[j] = data[i + j]
+                l = (r"$%s$" if label == "forward" else r"$%s^*$") % key
+                sensitivities[l] = S
+                i += dof
+        else:
+            S = np.zeros(6)
+            for j in range(n):
+                S[j] = data[i + j]
+            i += n
+            sensitivities[label.capitalize()] = S
+    return sensitivities
diff --git a/examples/pyroteus_burgers/testing_cases.txt b/examples/pyroteus_burgers/testing_cases.txt
new file mode 100644
index 0000000..1549b67
--- /dev/null
+++ b/examples/pyroteus_burgers/testing_cases.txt
@@ -0,0 +1 @@
+demo
diff --git a/examples/pyroteus_turbine/.DS_Store b/examples/pyroteus_turbine/.DS_Store
new file mode 100644
index 0000000..a90647e
Binary files /dev/null and b/examples/pyroteus_turbine/.DS_Store differ
diff --git a/examples/steady_turbine/config.py b/examples/pyroteus_turbine/config.py
similarity index 99%
rename from examples/steady_turbine/config.py
rename to examples/pyroteus_turbine/config.py
index 551d614..d21cd0a 100644
--- a/examples/steady_turbine/config.py
+++ b/examples/pyroteus_turbine/config.py
@@ -1,4 +1,4 @@
-from models.steady_turbine import *
+from models.pyroteus_turbine import *
 from nn_adapt.ann import sample_uniform
 
 
diff --git a/examples/pyroteus_turbine/meshgen.py b/examples/pyroteus_turbine/meshgen.py
new file mode 100644
index 0000000..d25c405
--- /dev/null
+++ b/examples/pyroteus_turbine/meshgen.py
@@ -0,0 +1,99 @@
+def generate_geo(config, reverse=False):
+    """
+    Given a configuration object for a given training
+    or testing case, generate a gmsh geometry file
+    defining the initial mesh, with all turbines
+    meshed explicitly.
+
+    :arg config: the configuration file
+    :kwarg reverse: should the flow direction be reversed?
+    """
+    if config.parameters.case in ("pipe", "headland"):
+        return
+    tc = config.parameters.turbine_coords
+    num_turbines = len(tc)
+    f = """// Domain and turbine specification
+L = 1200.0;
+W = 500.0;
+D = 18.0;
+dx_outer = 20.0;
+dx_inner = 20.0;
+"""
+    for i, xy in enumerate(tc):
+        f += "xt%d = %f;  // x-location of turbine %d\n" % (i, xy[0], i)
+        f += "yt%d = %f;  // y-location of turbine %d\n" % (i, xy[1], i)
+    f += """
+// Domain and turbine footprints
+Point(1) = {0, 0, 0, dx_outer};
+Point(2) = {L, 0, 0, dx_outer};
+Point(3) = {L, W, 0, dx_outer};
+Point(4) = {0, W, 0, dx_outer};
+Line(1) = {1, 2};
+Line(2) = {2, 3};
+Line(3) = {3, 4};
+Line(4) = {4, 1};
+Physical Line(1) = {%d};   // Left boundary
+Physical Line(2) = {%d};   // Right boundary
+Physical Line(3) = {1, 3}; // Sides
+Line Loop(1) = {1, 2, 3, 4};  // outside loop
+""" % (
+        2 if reverse else 4,
+        4 if reverse else 2,
+    )
+    i = 5
+    j = 2
+    for k in range(num_turbines):
+        f += """
+Point(%d) = {xt%d-D/2, yt%d-D/2, 0., dx_inner};
+Point(%d) = {xt%d+D/2, yt%d-D/2, 0., dx_inner};
+Point(%d) = {xt%d+D/2, yt%d+D/2, 0., dx_inner};
+Point(%d) = {xt%d-D/2, yt%d+D/2, 0., dx_inner};
+Line(%d) = {%d, %d};
+Line(%d) = {%d, %d};
+Line(%d) = {%d, %d};
+Line(%d) = {%d, %d};
+Line Loop(%d) = {%d, %d, %d, %d};
+""" % (
+            i,
+            k,
+            k,
+            i + 1,
+            k,
+            k,
+            i + 2,
+            k,
+            k,
+            i + 3,
+            k,
+            k,
+            i,
+            i,
+            i + 1,
+            i + 1,
+            i + 1,
+            i + 2,
+            i + 2,
+            i + 2,
+            i + 3,
+            i + 3,
+            i + 3,
+            i,
+            j,
+            i,
+            i + 1,
+            i + 2,
+            i + 3,
+        )
+        i += 4
+        j += 1
+    f += """
+// Surfaces
+Plane Surface(1) = %s;
+Physical Surface(1) = {1};  // outside turbines
+""" % set(
+        range(1, num_turbines + 2)
+    )
+    for i in range(1, num_turbines + 1):
+        f += "Plane Surface(%d) = {%d};\n" % (i + 1, i + 1)
+        f += "Physical Surface(%d) = {%d};  // inside turbine %d\n" % (i + 1, i + 1, i)
+    return f[:-1]
diff --git a/examples/steady_turbine/network.py b/examples/pyroteus_turbine/network.py
similarity index 100%
rename from examples/steady_turbine/network.py
rename to examples/pyroteus_turbine/network.py
diff --git a/examples/steady_turbine/plot_pipe.py b/examples/pyroteus_turbine/plot_pipe.py
similarity index 100%
rename from examples/steady_turbine/plot_pipe.py
rename to examples/pyroteus_turbine/plot_pipe.py
diff --git a/examples/steady_turbine/plotting.py b/examples/pyroteus_turbine/plotting.py
similarity index 100%
rename from examples/steady_turbine/plotting.py
rename to examples/pyroteus_turbine/plotting.py
diff --git a/examples/steady_turbine/testing_cases.txt b/examples/pyroteus_turbine/testing_cases.txt
similarity index 100%
rename from examples/steady_turbine/testing_cases.txt
rename to examples/pyroteus_turbine/testing_cases.txt
diff --git a/examples/run_adapt.py b/examples/run_adapt.py
index 69d9cd4..252bba5 100644
--- a/examples/run_adapt.py
+++ b/examples/run_adapt.py
@@ -47,9 +47,14 @@
 setup.initialise(test_case)
 unit = setup.parameters.qoi_unit
 if hasattr(setup, "initial_mesh"):
-    mesh = setup.initial_mesh
+    mesh = setup.initial_mesh()
 else:
     mesh = Mesh(f"{model}/meshes/{test_case}.msh")
+    
+try:
+    num_subinterval = len(mesh)
+except:
+    num_subinterval = 1
 
 # Run adaptation loop
 kwargs = {
@@ -62,7 +67,10 @@
     "h_max": setup.parameters.h_max,
     "a_max": 1.0e5,
 }
-ct = ConvergenceTracker(mesh, parsed_args)
+if num_subinterval == 1:
+    ct = ConvergenceTracker(mesh, parsed_args)
+else:
+    ct = ConvergenceTracker(mesh[0], parsed_args)
 if not no_outputs:
     output_dir = f"{model}/outputs/{test_case}/GO/{approach}"
     fwd_file = File(f"{output_dir}/forward.pvd")
@@ -87,7 +95,8 @@
     out = go_metric(mesh, setup, convergence_checker=ct, **kwargs)
     qoi, fwd_sol = out["qoi"], out["forward"]
     print(f"    Quantity of Interest = {qoi} {unit}")
-    dof = sum(np.array([fwd_sol.function_space().dof_count]).flatten())
+    spaces = [sol[0][0].function_space() for sol in fwd_sol.values()]
+    dof = sum(np.array([fs.dof_count for fs in spaces]).flatten())
     print(f"    DoF count            = {dof}")
     if "adjoint" not in out:
         break
@@ -97,32 +106,22 @@
         break
     adj_sol, dwr, metric = out["adjoint"], out["dwr"], out["metric"]
     if not no_outputs:
-        fwd_file.write(*fwd_sol.split())
-        adj_file.write(*adj_sol.split())
-        ee_file.write(dwr)
+        fields = ()
+        for sol in fwd_sol.values():
+            fields += sol[0][0].split()  # FIXME: Only uses 0th
+        fwd_file.write(*fields)
+        fields = ()
+        for sol in adj_sol.values():
+            fields += sol[0][0].split()  # FIXME: Only uses 0th
+        adj_file.write(*fields)
+        ee_file.write(dwr[-1])  # FIXME: Only uses 0th
         metric_file.write(metric.function)
 
-    def proj(V):
-        """
-        After the first iteration, project the previous
-        solution as the initial guess.
-        """
-        ic = Function(V)
-        try:
-            ic.project(fwd_sol)
-        except NotImplementedError:
-            for c_init, c in zip(ic.split(), fwd_sol.split()):
-                c_init.project(c)
-        return ic
-
-    # Use previous solution for initial guess
-    if parsed_args.transfer:
-        kwargs["init"] = proj
-
     # Extract features
     if not optimise:
-        features = extract_features(setup, fwd_sol, adj_sol)
-        target = dwr.dat.data.flatten()
+        field = list(fwd_sol.keys())[0]  # FIXME: Only uses 0th field
+        features = extract_features(setup, fwd_sol[field][0][-1], adj_sol[field][0][-1])  # FIXME
+        target = dwr[-1].dat.data.flatten()  # FIXME: Only uses 0th
         assert not np.isnan(target).any()
         for key, value in features.items():
             np.save(f"{data_dir}/feature_{key}_{suffix}", value)
@@ -130,12 +129,23 @@ def proj(V):
 
     # Adapt the mesh and check for element count convergence
     with PETSc.Log.Event("Mesh adaptation"):
-        mesh = adapt(mesh, metric)
+        if num_subinterval == 1:
+            mesh = adapt(mesh, metric)
+        else:
+            for id in range(num_subinterval):
+                mesh[id] = adapt(mesh[id], metric[id])
     if not no_outputs:
         mesh_file.write(mesh.coordinates)
-    elements = mesh.num_cells()
-    print(f"  Mesh {ct.fp_iteration+1}")
-    print(f"    Element count        = {elements}")
+        
+    if num_subinterval == 1:
+        elements = mesh.num_cells()
+        print(f"  Mesh {ct.fp_iteration+1}")
+        print(f"    Element count        = {elements}")
+    else:
+        elements_list = np.array([mesh_i.num_cells() for mesh_i in mesh])
+        elements = elements_list.mean()
+        print(f"  Mesh {ct.fp_iteration+1}")
+        print(f"    Element list        = {elements_list}")
     if ct.check_elements(elements):
         break
     ct.check_maxiter()
diff --git a/examples/run_adapt_ml.py b/examples/run_adapt_ml.py
index 9c79d29..f15c166 100644
--- a/examples/run_adapt_ml.py
+++ b/examples/run_adapt_ml.py
@@ -43,7 +43,7 @@
 setup.initialise(test_case)
 unit = setup.parameters.qoi_unit
 if hasattr(setup, "initial_mesh"):
-    mesh = setup.initial_mesh
+    mesh = setup.initial_mesh()
 else:
     mesh = Mesh(f"{model}/meshes/{test_case}.msh")
 
@@ -76,37 +76,29 @@
     out = get_solutions(mesh, setup, convergence_checker=ct, **kwargs)
     qoi, fwd_sol = out["qoi"], out["forward"]
     print(f"    Quantity of Interest = {qoi} {unit}")
-    dof = sum(np.array([fwd_sol.function_space().dof_count]).flatten())
+    spaces = [sol[0][0].function_space() for sol in fwd_sol.values()]
+    dof = sum(np.array([fs.dof_count for fs in spaces]).flatten())
     print(f"    DoF count            = {dof}")
     if "adjoint" not in out:
         break
     adj_sol = out["adjoint"]
     if not optimise:
-        fwd_file.write(*fwd_sol.split())
-        adj_file.write(*adj_sol.split())
+        fields = ()
+        for sol in fwd_sol.values():
+            fields += sol[0][0].split()  # FIXME: Only uses 0th
+        fwd_file.write(*fields)
+        fields = ()
+        for sol in adj_sol.values():
+            fields += sol[0][0].split()  # FIXME: Only uses 0th
+        adj_file.write(*fields)
     P0 = FunctionSpace(mesh, "DG", 0)
     P1_ten = TensorFunctionSpace(mesh, "CG", 1)
 
-    def proj(V):
-        """
-        After the first iteration, project the previous
-        solution as the initial guess.
-        """
-        ic = Function(V)
-        try:
-            ic.project(fwd_sol)
-        except NotImplementedError:
-            for c_init, c in zip(ic.split(), fwd_sol.split()):
-                c_init.project(c)
-        return ic
-
-    # Use previous solution for initial guess
-    if parsed_args.transfer:
-        kwargs["init"] = proj
-
     # Extract features
     with PETSc.Log.Event("Network"):
-        features = collect_features(extract_features(setup, fwd_sol, adj_sol), layout)
+        field = list(fwd_sol.keys())[0]  # FIXME: Only uses 0th field
+        features = extract_features(setup, fwd_sol[field][0][0], adj_sol[field][0][0])  # FIXME
+        features = collect_features(features, layout)
 
         # Run model
         with PETSc.Log.Event("Propagate"):
@@ -120,6 +112,7 @@ def proj(V):
                     )
             dwr = Function(P0)
             dwr.dat.data[:] = np.abs(test_targets)
+            # FIXME: Only produces one error indicator
 
     # Check for error estimator convergence
     with PETSc.Log.Event("Error estimation"):
@@ -133,7 +126,10 @@ def proj(V):
     # Construct metric
     with PETSc.Log.Event("Metric construction"):
         if approach == "anisotropic":
-            hessian = combine_metrics(*get_hessians(fwd_sol), average=True)
+            field = list(fwd_sol.keys())[0]
+            fwd = fwd_sol[field][0]  # FIXME: Only uses 0th
+            hessians = sum([get_hessians(sol) for sol in fwd], start=())
+            hessian = combine_metrics(*hessians, average=True)
         else:
             hessian = None
         M = anisotropic_metric(
diff --git a/examples/run_adaptation_loop.py b/examples/run_adaptation_loop.py
index 126616a..91c508e 100644
--- a/examples/run_adaptation_loop.py
+++ b/examples/run_adaptation_loop.py
@@ -49,8 +49,7 @@
 # Run adaptation loop
 qois, dofs, elements, estimators, niter = [], [], [], [], []
 components = ("forward", "adjoint", "estimator", "metric", "adapt")
-times = {c: [] for c in components}
-times["all"] = []
+times = []
 print(f"Test case {test_case}")
 for i in range(num_refinements + 1):
     try:
@@ -66,15 +65,13 @@
             "a_max": 1.0e5,
         }
         if hasattr(setup, "initial_mesh"):
-            mesh = setup.initial_mesh
+            mesh = setup.initial_mesh()
         else:
             mesh = Mesh(f"{model}/meshes/{test_case}.msh")
         ct = ConvergenceTracker(mesh, parsed_args)
         print(f"  Target {target_complexity}\n    Mesh 0")
         print(f"      Element count        = {ct.elements_old}")
-        times["all"].append(-perf_counter())
-        for c in components:
-            times[c].append(0.0)
+        times.append(-perf_counter())
         for ct.fp_iteration in range(ct.maxiter + 1):
 
             # Ramp up the target complexity
@@ -85,47 +82,24 @@
             # Compute goal-oriented metric
             out = go_metric(mesh, setup, convergence_checker=ct, **kwargs)
             qoi = out["qoi"]
-            times["forward"][-1] += out["times"]["forward"]
             print(f"      Quantity of Interest = {qoi} {unit}")
             if "adjoint" not in out:
                 break
             estimator = out["estimator"]
-            times["adjoint"][-1] += out["times"]["adjoint"]
-            times["estimator"][-1] += out["times"]["estimator"]
             print(f"      Error estimator      = {estimator}")
             if "metric" not in out:
                 break
-            times["metric"][-1] += out["times"]["metric"]
             fwd_sol, adj_sol = (
                 out["forward"],
                 out["adjoint"],
             )
             dwr, metric = out["dwr"], out["metric"]
-            dof = sum(np.array([fwd_sol.function_space().dof_count]).flatten())
+            spaces = [sol[0][0].function_space() for sol in fwd_sol.values()]
+            dof = sum(np.array([fs.dof_count for fs in spaces]).flatten())
             print(f"      DoF count            = {dof}")
 
-            def proj(V):
-                """
-                After the first iteration, project the previous
-                solution as the initial guess.
-                """
-                ic = Function(V)
-                try:
-                    ic.project(fwd_sol)
-                except NotImplementedError:
-                    for c_init, c in zip(ic.split(), fwd_sol.split()):
-                        c_init.project(c)
-                return ic
-
-            # Use previous solution for initial guess
-            if parsed_args.transfer:
-                kwargs["init"] = proj
-
             # Adapt the mesh
-            out["times"]["adapt"] = -perf_counter()
             mesh = adapt(mesh, metric)
-            out["times"]["adapt"] += perf_counter()
-            times["adapt"][-1] += out["times"]["adapt"]
             print(f"    Mesh {ct.fp_iteration+1}")
             cells = mesh.num_cells()
             print(f"      Element count        = {cells}")
@@ -135,7 +109,7 @@ def proj(V):
         print(
             f"    Terminated after {ct.fp_iteration+1} iterations due to {ct.converged_reason}"
         )
-        times["all"][-1] += perf_counter()
+        times[-1] += perf_counter()
         qois.append(qoi)
         dofs.append(dof)
         elements.append(cells)
@@ -146,9 +120,7 @@ def proj(V):
         np.save(f"{model}/data/elements_GO{approach}_{test_case}", elements)
         np.save(f"{model}/data/estimators_GO{approach}_{test_case}", estimators)
         np.save(f"{model}/data/niter_GO{approach}_{test_case}", niter)
-        np.save(f"{model}/data/times_all_GO{approach}_{test_case}", times["all"])
-        for c in components:
-            np.save(f"{model}/data/times_{c}_GO{approach}_{test_case}", times[c])
+        np.save(f"{model}/data/times_all_GO{approach}_{test_case}", times)
     except ConvergenceError:
         print("Skipping due to convergence error")
         continue
diff --git a/examples/run_adaptation_loop_ml.py b/examples/run_adaptation_loop_ml.py
index 018329b..0011264 100644
--- a/examples/run_adaptation_loop_ml.py
+++ b/examples/run_adaptation_loop_ml.py
@@ -60,23 +60,20 @@
 # Run adaptation loop
 qois, dofs, elements, estimators, niter = [], [], [], [], []
 components = ("forward", "adjoint", "estimator", "metric", "adapt")
-times = {c: [] for c in components}
-times["all"] = []
+times = []
 print(f"Test case {test_case}")
 for i in range(num_refinements + 1):
     try:
         target_complexity = 100.0 * 2 ** (f * i)
         if hasattr(setup, "initial_mesh"):
-            mesh = setup.initial_mesh
+            mesh = setup.initial_mesh()
         else:
             mesh = Mesh(f"{model}/meshes/{test_case}.msh")
         ct = ConvergenceTracker(mesh, parsed_args)
         kwargs = {}
         print(f"  Target {target_complexity}\n    Mesh 0")
         print(f"      Element count        = {ct.elements_old}")
-        times["all"].append(-perf_counter())
-        for c in components:
-            times[c].append(0.0)
+        times.append(-perf_counter())
         for ct.fp_iteration in range(ct.maxiter + 1):
 
             # Ramp up the target complexity
@@ -87,35 +84,17 @@
             # Solve forward and adjoint and compute Hessians
             out = get_solutions(mesh, setup, convergence_checker=ct, **kwargs)
             qoi = out["qoi"]
-            times["forward"][-1] += out["times"]["forward"]
             print(f"      Quantity of Interest = {qoi} {unit}")
             if "adjoint" not in out:
                 break
-            times["adjoint"][-1] += out["times"]["adjoint"]
             fwd_sol, adj_sol = out["forward"], out["adjoint"]
-            dof = sum(np.array([fwd_sol.function_space().dof_count]).flatten())
+            spaces = [sol[0][0].function_space() for sol in fwd_sol.values()]
+            dof = sum(np.array([fs.dof_count for fs in spaces]).flatten())
             print(f"      DoF count            = {dof}")
 
-            def proj(V):
-                """
-                After the first iteration, project the previous
-                solution as the initial guess.
-                """
-                ic = Function(V)
-                try:
-                    ic.project(fwd_sol)
-                except NotImplementedError:
-                    for c_init, c in zip(ic.split(), fwd_sol.split()):
-                        c_init.project(c)
-                return ic
-
-            # Use previous solution for initial guess
-            if parsed_args.transfer:
-                kwargs["init"] = proj
-
             # Extract features
-            out["times"]["estimator"] = -perf_counter()
-            features = extract_features(setup, fwd_sol, adj_sol)
+            field = list(fwd_sol.keys())[0]  # FIXME: Only uses 0th field
+            features = extract_features(setup, fwd_sol[field][0][0], adj_sol[field][0][0])  # FIXME
             features = collect_features(features, layout)
 
             # Run model
@@ -130,19 +109,20 @@ def proj(V):
             P0 = FunctionSpace(mesh, "DG", 0)
             dwr = Function(P0)
             dwr.dat.data[:] = np.abs(test_targets)
+            # FIXME: Only produces one error indicator
 
             # Check for error estimator convergence
             estimator = dwr.vector().gather().sum()
-            out["times"]["estimator"] += perf_counter()
-            times["estimator"][-1] += out["times"]["estimator"]
             print(f"      Error estimator      = {estimator}")
             if ct.check_estimator(estimator):
                 break
 
             # Construct metric
-            out["times"]["metric"] = -perf_counter()
             if approach == "anisotropic":
-                hessian = combine_metrics(*get_hessians(fwd_sol), average=True)
+                field = list(fwd_sol.keys())[0]
+                fwd = fwd_sol[field][0]  # FIXME: Only uses 0th
+                hessians = sum([get_hessians(sol) for sol in fwd], start=())
+                hessian = combine_metrics(*hessians, average=True)
             else:
                 hessian = None
             P1_ten = TensorFunctionSpace(mesh, "CG", 1)
@@ -159,14 +139,9 @@ def proj(V):
             )
             metric = RiemannianMetric(mesh)
             metric.assign(M)
-            out["times"]["metric"] += perf_counter()
-            times["metric"][-1] += out["times"]["metric"]
 
             # Adapt the mesh and check for element count convergence
-            out["times"]["adapt"] = -perf_counter()
             mesh = adapt(mesh, metric)
-            out["times"]["adapt"] += perf_counter()
-            times["adapt"][-1] += out["times"]["adapt"]
             print(f"    Mesh {ct.fp_iteration+1}")
             cells = mesh.num_cells()
             print(f"      Element count        = {cells}")
@@ -176,7 +151,7 @@ def proj(V):
         print(
             f"    Terminated after {ct.fp_iteration+1} iterations due to {ct.converged_reason}"
         )
-        times["all"][-1] += perf_counter()
+        times[-1] += perf_counter()
         qois.append(qoi)
         dofs.append(dof)
         elements.append(cells)
@@ -187,9 +162,7 @@ def proj(V):
         np.save(f"{model}/data/elements_ML{approach}_{test_case}_{tag}", elements)
         np.save(f"{model}/data/estimators_ML{approach}_{test_case}_{tag}", estimators)
         np.save(f"{model}/data/niter_ML{approach}_{test_case}_{tag}", niter)
-        np.save(f"{model}/data/times_all_ML{approach}_{test_case}_{tag}", times["all"])
-        for c in components:
-            np.save(f"{model}/data/times_{c}_ML{approach}_{test_case}_{tag}", times[c])
+        np.save(f"{model}/data/times_all_ML{approach}_{test_case}_{tag}", times)
     except ConvergenceError:
         print("Skipping due to convergence error")
         continue
diff --git a/examples/run_fixed_mesh.py b/examples/run_fixed_mesh.py
index f8ad698..65a64fd 100644
--- a/examples/run_fixed_mesh.py
+++ b/examples/run_fixed_mesh.py
@@ -39,10 +39,12 @@
 qoi = out["qoi"]
 print_output(f"QoI for test case {test_case} = {qoi:.8f} {unit}")
 if not parsed_args.optimise:
-    File(f"{model}/outputs/{test_case}/fixed/forward.pvd").write(
-        *out["forward"].split()
-    )
-    File(f"{model}/outputs/{test_case}/fixed/adjoint.pvd").write(
-        *out["adjoint"].split()
-    )
+    fields = ()
+    for f, sol in out["forward"].items():
+        fields += sol[0][0].split()  # FIXME: only 0th timestep
+    File(f"{model}/outputs/{test_case}/fixed/forward.pvd").write(*fields)
+    fields = ()
+    for f, sol in out["adjoint"].items():
+        fields += sol[0][0].split()  # FIXME: only 0th timestep
+    File(f"{model}/outputs/{test_case}/fixed/adjoint.pvd").write(*fields)
 print_output(f"  Total time taken: {perf_counter() - start_time:.2f} seconds")
diff --git a/examples/run_uniform_refinement.py b/examples/run_uniform_refinement.py
index be3b3a8..cc58ca6 100644
--- a/examples/run_uniform_refinement.py
+++ b/examples/run_uniform_refinement.py
@@ -31,7 +31,10 @@
 setup = importlib.import_module(f"{model}.config")
 setup.initialise(test_case)
 unit = setup.parameters.qoi_unit
-mesh = Mesh(f"{model}/meshes/{test_case}.msh")
+if hasattr(setup, "initial_mesh"):
+    mesh = setup.initial_mesh()
+else:
+    mesh = Mesh(f"{model}/meshes/{test_case}.msh")
 mh = [mesh] + list(MeshHierarchy(mesh, num_refinements))
 tm = TransferManager()
 kwargs = {}
@@ -59,12 +62,13 @@ def prolong(V):
 
     if parsed_args.prolong:
         kwargs["init"] = prolong
-    fs = fwd_sol.function_space()
+    spaces = [sol[0][0].function_space() for sol in fwd_sol.values()]
     time = perf_counter() - start_time
     print_output(f"    Quantity of Interest = {qoi} {unit}")
     print_output(f"    Runtime: {time:.2f} seconds")
     qois.append(qoi)
-    dofs.append(sum(fs.dof_count))
+    dof = sum(np.array([fs.dof_count for fs in spaces]).flatten())
+    dofs.append(dof)
     times.append(time)
     elements.append(mesh.num_cells())
     niter.append(1)
diff --git a/examples/steady_turbine/meshes/headland.geo b/examples/steady_turbine/meshes/headland.geo
deleted file mode 100644
index 293e4a4..0000000
--- a/examples/steady_turbine/meshes/headland.geo
+++ /dev/null
@@ -1,56 +0,0 @@
-L = 1200.0;
-W = 500.0;
-D = 18.0;
-dx_outer = 20.0;
-dx_inner = 20.0;
-
-headland_x_scale = 0.2;
-headland_y = 150;
-
-site_x = 100;
-site_y = 80;
-site_x_start = L/2 - site_x/2;
-site_x_end = site_x_start + site_x;
-
-site_y_start = W/2 - site_y/2;
-site_y_end = site_y_start + site_y;
-
-Point(1) = {0, 0, 0, dx_outer};
-Point(2) = {L, 0, 0, dx_outer};
-
-
-// Headland
-res = 100;
-b = 10;
-For k In {0:res:1}
-    x = L/res*k;
-    y = W - headland_y*Exp(-0.5*((headland_x_scale*(x-L/2))/b)^2);
-	Point(10+k) = {x, y, 0, dx_outer};
-EndFor
-BSpline(100) = { 10 : res+10 };
-
-// Domain boundary
-Line(101) = {10, 1};
-Line(102) = {1, 2};
-Line(103) = {2, res+10};
-Line Loop(104) = {100, -103, -102, -101};
-
-// Generate site nodes
-Point(111) = {site_x_start, site_y_start, 0, dx_inner};
-Point(112) = {site_x_end, site_y_start, 0, dx_inner};
-Point(113) = {site_x_end, site_y_end, 0, dx_inner};
-Point(114) = {site_x_start, site_y_end, 0, dx_inner};
-Line(105) = {111, 112};
-Line(106) = {112, 113};
-Line(107) = {113, 114};
-Line(108) = {114, 111};
-Line Loop(110) = {105, 106, 107, 108};
-Plane Surface(111) = {104, 110};
-Plane Surface(112) = {110};
-Physical Line(1) = {101};  // inflow
-Physical Line(2) = {103};  // outflow
-Physical Line(3) = {100, 102};  // free-slip
-Physical Line(4) = {};  // no-slip
-Physical Line(5) = {};  // weakly reflective
-Physical Surface(1) = {111};
-Physical Surface(2) = {112};
diff --git a/examples/steady_turbine/meshes/pipe.geo b/examples/steady_turbine/meshes/pipe.geo
deleted file mode 100644
index da211d0..0000000
--- a/examples/steady_turbine/meshes/pipe.geo
+++ /dev/null
@@ -1,72 +0,0 @@
-D = 18.0;
-dx_outer = 20.0;
-dx_inner = 20.0;
-xt0 = 550.0;  // x-location of turbine 0
-yt0 = 300.0;  // y-location of turbine 0
-xt1 = 620.0;  // x-location of turbine 1
-yt1 = 390.0;  // y-location of turbine 1
-W = 700;  // Width of domain
-w = 200;  // Width of channel
-
-sqrt2 = 1.414213562373095;
-
-// Lower points
-Point(1) = {0,    0,   0, dx_outer};
-Point(2) = {200,  0,   0, dx_outer}; // First spline point
-Point(3) = {400,  0,   0, dx_outer};
-Point(4) = {600,  W/2-0.75*w, 0, dx_outer};
-Point(5) = {800,  W-w, 0, dx_outer};
-Point(6) = {1000, W-w, 0, dx_outer}; // Last spline point
-Point(7) = {1200, W-w, 0, dx_outer};
-// Upper points
-Point( 8) = {1200, W, 0, dx_outer};
-Point( 9) = {1000, W, 0, dx_outer}; // First spline point
-Point(10) = {800,  W, 0, dx_outer};
-Point(11) = {600,  W/2+0.75*w, 0, dx_outer};
-Point(12) = {400,  w, 0, dx_outer};
-Point(13) = {200,  w, 0, dx_outer}; // Last spline point
-Point(14) = {0,    w, 0, dx_outer};
-
-// Edges
-Line(1) = {1, 2};
-BSpline(2) = {2, 3, 4, 5, 6};
-Line(3) = {6, 7};
-Line(4) = {7, 8};
-Line(5) = {8, 9};
-BSpline(6) = {9, 10, 11, 12, 13};
-Line(7) = {13, 14};
-Line(8) = {14,  1};
-
-// Boundary and physical curves
-Curve Loop(9) = {1, 2, 3, 4, 5, 6, 7, 8};
-Physical Curve("Inflow", 1) = {8};
-Physical Curve("Outflow", 2) = {4};
-Physical Curve("Sides", 4) = {1, 2, 3, 5, 6, 7};
-
-Point(15) = {xt0-D/sqrt2, yt0, 0., dx_inner};
-Point(16) = {xt0, yt0-D/sqrt2, 0., dx_inner};
-Point(17) = {xt0+D/sqrt2, yt0, 0., dx_inner};
-Point(18) = {xt0, yt0+D/sqrt2, 0., dx_inner};
-Line(15) = {15, 16};
-Line(16) = {16, 17};
-Line(17) = {17, 18};
-Line(18) = {18, 15};
-Line Loop(10) = {15, 16, 17, 18};
-
-Point(19) = {xt1-D/sqrt2, yt1, 0., dx_inner};
-Point(20) = {xt1, yt1-D/sqrt2, 0., dx_inner};
-Point(21) = {xt1+D/sqrt2, yt1, 0., dx_inner};
-Point(22) = {xt1, yt1+D/sqrt2, 0., dx_inner};
-Line(19) = {19, 20};
-Line(20) = {20, 21};
-Line(21) = {21, 22};
-Line(22) = {22, 19};
-Line Loop(11) = {19, 20, 21, 22};
-
-// Domain and physical surface
-Plane Surface(1) = {9, 10, 11};
-Physical Surface("Pipe", 1) = {1};
-Plane Surface(2) = {10};
-Physical Surface(2) = {2};  // inside turbine 1
-Plane Surface(3) = {11};
-Physical Surface(3) = {3};  // inside turbine 2
diff --git a/examples/test_and_train.py b/examples/test_and_train.py
index e3cfd9b..992bdcc 100644
--- a/examples/test_and_train.py
+++ b/examples/test_and_train.py
@@ -197,7 +197,7 @@
                 key: np.load(f"{data_dir}/feature_{key}_{suffix}.npy")
                 for key in layout.inputs
             }
-            features = concat(features, collect_features(feature))
+            features = concat(features, collect_features(feature, layout))
             target = np.load(f"{data_dir}/target_{suffix}.npy")
             targets = concat(targets, target)
 print(f"Total number of features: {len(features.flatten())}")
diff --git a/nn_adapt/.DS_Store b/nn_adapt/.DS_Store
new file mode 100644
index 0000000..5008ddf
Binary files /dev/null and b/nn_adapt/.DS_Store differ
diff --git a/nn_adapt/features.py b/nn_adapt/features.py
index f185f89..2f0f23d 100644
--- a/nn_adapt/features.py
+++ b/nn_adapt/features.py
@@ -8,7 +8,7 @@
 import numpy as np
 from pyroteus.metric import *
 import ufl
-from nn_adapt.solving import dwr_indicator
+# from nn_adapt.solving import dwr_indicator
 from collections import Iterable
 
 
@@ -215,9 +215,9 @@ def extract_features(config, fwd_sol, adj_sol):
     """
     mesh = fwd_sol.function_space().mesh()
 
-    # Coarse-grained DWR estimator
-    with PETSc.Log.Event("Extract estimator"):
-        dwr = dwr_indicator(config, mesh, fwd_sol, adj_sol)
+    # # Coarse-grained DWR estimator  # TODO: Reintroduce
+    # with PETSc.Log.Event("Extract estimator"):
+    #     dwr = dwr_indicator(config, mesh, fwd_sol, adj_sol)
 
     # Features describing the mesh element
     with PETSc.Log.Event("Analyse element"):
@@ -229,12 +229,12 @@ def extract_features(config, fwd_sol, adj_sol):
         d, h1, h2 = (extract_array(p) for p in extract_components(JTJ))
 
         # Is the element on the boundary?
-        p0test = firedrake.TestFunction(dwr.function_space())
+        p0test = firedrake.TestFunction(firedrake.FunctionSpace(mesh, "DG", 0))
         bnd = firedrake.assemble(p0test * ufl.ds).dat.data
 
     # Combine the features together
     features = {
-        "estimator_coarse": extract_array(dwr),
+        # "estimator_coarse": extract_array(dwr),  # TODO: Reintroduce
         "physics_drag": extract_array(config.parameters.drag(mesh)),
         "physics_viscosity": extract_array(config.parameters.viscosity(mesh), project=True),
         "physics_bathymetry": extract_array(config.parameters.bathymetry(mesh), project=True),
diff --git a/nn_adapt/layout.py b/nn_adapt/layout.py
index 060502f..15e0651 100644
--- a/nn_adapt/layout.py
+++ b/nn_adapt/layout.py
@@ -13,10 +13,12 @@ class NetLayoutBase(object):
     for each of these parameters.
     """
 
+    inputs = None
+    num_hidden_neurons = None
     # TODO: Allow more general networks
 
     colours = {
-        "estimator": "b",
+        # "estimator": "b",
         "physics": "C0",
         "mesh": "deepskyblue",
         "forward": "mediumturquoise",
diff --git a/nn_adapt/metric.py b/nn_adapt/metric.py
index 22ed97e..634e147 100644
--- a/nn_adapt/metric.py
+++ b/nn_adapt/metric.py
@@ -6,7 +6,6 @@
 from nn_adapt.features import split_into_scalars
 from nn_adapt.solving import *
 from firedrake.meshadapt import RiemannianMetric
-from time import perf_counter
 
 
 def get_hessians(f, **kwargs):
@@ -22,11 +21,11 @@ def get_hessians(f, **kwargs):
         component
     """
     kwargs.setdefault("method", "Clement")
-    return [
+    return tuple(
         space_normalise(hessian_metric(recover_hessian(fij, **kwargs)), 4000.0, "inf")
         for i, fi in split_into_scalars(f).items()
         for fij in fi
-    ]
+    )
 
 
 def go_metric(
@@ -68,6 +67,11 @@ def go_metric(
     :kwarg convergence_checker: :class:`ConvergenceTracer`
         instance
     """
+    try:
+        num_subintervals = len(mesh)
+    except:
+        num_subintervals = 1
+        mesh = [mesh]
     h_min = kwargs.pop("h_min", 1.0e-30)
     h_max = kwargs.pop("h_max", 1.0e+30)
     a_max = kwargs.pop("a_max", 1.0e+30)
@@ -81,27 +85,67 @@ def go_metric(
     )
     if retall and "adjoint" not in out:
         return out
-    out["estimator"] = out["dwr"].vector().gather().sum()
-    if convergence_checker is not None:
-        if convergence_checker.check_estimator(out["estimator"]):
-            return out
+    
+    # single mesh for whole time interval
+    if num_subintervals == 1:
+        out["estimator"] = out["dwr"][0].vector().gather().sum()
+        if convergence_checker is not None:
+            if convergence_checker.check_estimator(out["estimator"]):
+                return out
+            
+    # multiple meshes for whole time interval
+    else:
+        out["estimator"] = [0 for _ in range(num_subintervals)]
+        for id in range(num_subintervals):
+            out["estimator"][id] = out["dwr"][id].vector().gather().sum()
+        if convergence_checker is not None:
+            max_estimator = np.array(out["estimator"]).mean()
+            if convergence_checker.check_estimator(max_estimator):
+                return out
 
-    out["times"]["metric"] = -perf_counter()
     with PETSc.Log.Event("Metric construction"):
-        if anisotropic:
-            hessian = combine_metrics(*get_hessians(out["forward"]), average=average)
-        else:
-            hessian = None
-        metric = anisotropic_metric(
-            out["dwr"],
-            hessian=hessian,
-            target_complexity=target_complexity,
-            target_space=TensorFunctionSpace(mesh, "CG", 1),
-            interpolant=interpolant,
-        )
-        space_normalise(metric, target_complexity, "inf")
-        enforce_element_constraints(metric, h_min, h_max, a_max)
-        out["metric"] = RiemannianMetric(mesh)
-        out["metric"].assign(metric)
-    out["times"]["metric"] += perf_counter()
+        # single mesh for whole time interval
+        if num_subintervals == 1:
+            if anisotropic:
+                field = list(out["forward"].keys())[0]
+                fwd = out["forward"][field][0]
+                hessians = sum([get_hessians(sol) for sol in fwd], start=())
+                hessian = combine_metrics(*hessians, average=average)
+            else:
+                hessian = None
+            metric = anisotropic_metric(
+                out["dwr"][0],
+                hessian=hessian,
+                target_complexity=target_complexity,
+                target_space=TensorFunctionSpace(mesh[0], "CG", 1),
+                interpolant=interpolant,
+            )
+            space_normalise(metric, target_complexity, "inf")
+            enforce_element_constraints(metric, h_min, h_max, a_max)
+            out["metric"] = RiemannianMetric(mesh[0])
+            out["metric"].assign(metric)
+        
+        # multiple meshes for whole time interval
+        else: 
+            out["metric"] = []
+            for id in range(num_subintervals):
+                if anisotropic:
+                    field = list(out["forward"].keys())[0]
+                    fwd = out["forward"][field][id]
+                    hessians = sum([get_hessians(sol) for sol in fwd], start=())
+                    hessian = combine_metrics(*hessians, average=average)
+                else:
+                    hessian = None
+                metric = anisotropic_metric(
+                    out["dwr"][id],
+                    hessian=hessian,
+                    target_complexity=target_complexity,
+                    target_space=TensorFunctionSpace(mesh[id], "CG", 1),
+                    interpolant=interpolant,
+                )
+                space_normalise(metric, target_complexity, "inf")
+                enforce_element_constraints(metric, h_min, h_max, a_max)
+                out["metric"].append(RiemannianMetric(mesh[id]))
+                out["metric"][-1].assign(metric)
+                
     return out if retall else out["metric"]
diff --git a/nn_adapt/parse.py b/nn_adapt/parse.py
index 7565690..12da786 100644
--- a/nn_adapt/parse.py
+++ b/nn_adapt/parse.py
@@ -105,11 +105,6 @@ def parse_approach(self):
             choices=["isotropic", "anisotropic"],
             default="anisotropic",
         )
-        self.add_argument(
-            "--transfer",
-            help="Transfer the solution from the previous mesh as initial guess",
-            action="store_true",
-        )
 
     def parse_target_complexity(self):
         self.add_argument(
diff --git a/nn_adapt/solving.py b/nn_adapt/solving.py
index b7a950c..ccb9876 100644
--- a/nn_adapt/solving.py
+++ b/nn_adapt/solving.py
@@ -3,69 +3,13 @@
 files and performing goal-oriented error estimation.
 """
 from firedrake import *
-from firedrake.petsc import PETSc
-from firedrake.mg.embedded import TransferManager
-from pyroteus.error_estimation import get_dwr_indicator
-import abc
-from time import perf_counter
-
-
-tm = TransferManager()
-
-
-class Solver(abc.ABC):
-    """
-    Base class that defines the API for solver objects.
-    """
-
-    @abc.abstractmethod
-    def __init__(self, mesh, ic, **kwargs):
-        """
-        Setup the solver.
-
-        :arg mesh: the mesh to define the solver on
-        :arg ic: the initial condition
-        """
-        pass
-
-    @property
-    @abc.abstractmethod
-    def function_space(self):
-        """
-        The function space that the PDE is solved in.
-        """
-        pass
-
-    @property
-    @abc.abstractmethod
-    def form(self):
-        """
-        Return the weak form.
-        """
-        pass
-
-    @abc.abstractmethod
-    def iterate(self, **kwargs):
-        """
-        Solve the PDE.
-        """
-        pass
-
-    @property
-    @abc.abstractmethod
-    def solution(self):
-        """
-        Return the solution field.
-        """
-        pass
+from firedrake_adjoint import *
 
 
 def get_solutions(
     mesh,
     config,
     solve_adjoint=True,
-    refined_mesh=None,
-    init=None,
     convergence_checker=None,
     **kwargs,
 ):
@@ -82,67 +26,32 @@ def get_solutions(
         adjoint problem?
     :kwarg refined_mesh: refined mesh to compute
         enriched adjoint solution on
-    :kwarg init: custom initial condition function
     :kwarg convergence_checker: :class:`ConvergenceTracer`
         instance
     :return: forward solution, adjoint solution
         and enriched adjoint solution (if requested)
     """
+    out = {}
+    solve_adjoint = True  # TODO: Account for the false case!
+    # NOTE: None of the timings will work!
 
     # Solve forward problem in base space
-    V = config.get_function_space(mesh)
-    out = {"times": {"forward": -perf_counter()}}
-    with PETSc.Log.Event("Forward solve"):
-        if init is None:
-            ic = config.get_initial_condition(V)
-        else:
-            ic = init(V)
-        solver_obj = config.Solver(mesh, ic, **kwargs)
-        solver_obj.iterate()
-        q = solver_obj.solution
-        J = config.get_qoi(mesh)(q)
-        qoi = assemble(J)
-    out["times"]["forward"] += perf_counter()
+    mesh_seq = config.GoalOrientedMeshSeq(mesh, **kwargs)
+    solutions = mesh_seq.solve_adjoint()
+    fields = mesh_seq.fields
+    qoi = mesh_seq.J
     out["qoi"] = qoi
-    out["forward"] = q
+    out["forward"] = {f: solutions[f]["forward"] for f in fields}
     if convergence_checker is not None:
         if convergence_checker.check_qoi(qoi):
             return out
+    else:
+        print("No convergence checker")
     if not solve_adjoint:
         return out
 
     # Solve adjoint problem in base space
-    out["times"]["adjoint"] = -perf_counter()
-    with PETSc.Log.Event("Adjoint solve"):
-        sp = config.parameters.adjoint_solver_parameters
-        q_star = Function(V)
-        F = solver_obj.form
-        dFdq = derivative(F, q, TrialFunction(V))
-        dFdq_transpose = adjoint(dFdq)
-        dJdq = derivative(J, q, TestFunction(V))
-        solve(dFdq_transpose == dJdq, q_star, solver_parameters=sp)
-        out["adjoint"] = q_star
-    out["times"]["adjoint"] += perf_counter()
-    if refined_mesh is None:
-        return out
-
-    # Solve adjoint problem in enriched space
-    out["times"]["estimator"] = -perf_counter()
-    with PETSc.Log.Event("Enrichment"):
-        V = config.get_function_space(refined_mesh)
-        q_plus = Function(V)
-        solver_obj = config.Solver(refined_mesh, q_plus, **kwargs)
-        q_plus = solver_obj.solution
-        J = config.get_qoi(refined_mesh)(q_plus)
-        F = solver_obj.form
-        tm.prolong(q, q_plus)
-        q_star_plus = Function(V)
-        dFdq = derivative(F, q_plus, TrialFunction(V))
-        dFdq_transpose = adjoint(dFdq)
-        dJdq = derivative(J, q_plus, TestFunction(V))
-        solve(dFdq_transpose == dJdq, q_star_plus, solver_parameters=sp)
-        out["enriched_adjoint"] = q_star_plus
-    out["times"]["estimator"] += perf_counter()
+    out["adjoint"] = {f: solutions[f]["adjoint"] for f in fields}
     return out
 
 
@@ -154,10 +63,10 @@ def split_into_components(f):
     return [f] if f.function_space().value_size == 1 else f.split()
 
 
-def indicate_errors(mesh, config, enrichment_method="h", retall=False, **kwargs):
+def indicate_errors(mesh, config, enrichment_method="h", retall=False, convergence_checker = None, **kwargs):
     """
-    Indicate errors according to ``dwr_indicator``,
-    using the solver given in the configuration file.
+    Indicate errors according to the ``GoalOrientedMeshSeq``
+    given in the configuration file.
 
     :arg mesh: input mesh
     :arg config: configuration file, which
@@ -168,59 +77,28 @@ def indicate_errors(mesh, config, enrichment_method="h", retall=False, **kwargs)
         solution and adjoint solution in addition
         to the dual-weighted residual error indicator
     """
+    out = {}
     if not enrichment_method == "h":
         raise NotImplementedError  # TODO
-    with PETSc.Log.Event("Enrichment"):
-        mesh, ref_mesh = MeshHierarchy(mesh, 1)
-
-    # Solve the forward and adjoint problems
-    out = get_solutions(mesh, config, refined_mesh=ref_mesh, **kwargs)
-    if retall and "adjoint" not in out:
-        return out
-
-    out["times"]["estimator"] -= perf_counter()
-    with PETSc.Log.Event("Enrichment"):
-        adj_sol_plus = out["enriched_adjoint"]
-
-        # Prolong
-        V_plus = adj_sol_plus.function_space()
-        fwd_sol_plg = Function(V_plus)
-        tm.prolong(out["forward"], fwd_sol_plg)
-        adj_sol_plg = Function(V_plus)
-        tm.prolong(out["adjoint"], adj_sol_plg)
-
-        # Subtract prolonged adjoint solution from enriched version
-        adj_error = Function(V_plus)
-        adj_sols_plus = split_into_components(adj_sol_plus)
-        adj_sols_plg = split_into_components(adj_sol_plg)
-        for i, err in enumerate(split_into_components(adj_error)):
-            err += adj_sols_plus[i] - adj_sols_plg[i]
-
-        # Evaluate errors
-        out["dwr"] = dwr_indicator(config, mesh, fwd_sol_plg, adj_error)
-    out["times"]["estimator"] += perf_counter()
-
+    mesh_seq = config.GoalOrientedMeshSeq(mesh, **kwargs)
+    fields = mesh_seq.fields
+    kw = {"enrichment_method": enrichment_method}
+    solutions, indicators = mesh_seq.indicate_errors(enrichment_kwargs=kw)
+    integrated = []
+    for i, indicator in enumerate(indicators):
+        dt = mesh_seq.time_partition.timesteps[i]
+        contrib = Function(indicator[0].function_space())
+        for indi in indicator:
+            contrib += dt * indi
+        integrated.append(contrib)
+    qoi = mesh_seq.J
+    out["qoi"] = qoi
+    out["forward"] = {f: solutions[f]["forward"] for f in fields}
+    out["adjoint"] = {f: solutions[f]["adjoint"] for f in fields}
+    out["dwr"] = integrated
+    if convergence_checker is not None:
+        if convergence_checker.check_qoi(qoi):
+            return out
+    else:
+        print("No convergence checker")
     return out if retall else out["dwr"]
-
-
-def dwr_indicator(config, mesh, q, q_star):
-    r"""
-    Evaluate the DWR error indicator as a :math:`\mathbb P0` field.
-
-    :arg mesh: the current mesh
-    :arg q: the forward solution, transferred into enriched space
-    :arg q_star: the adjoint solution in enriched space
-    """
-    mesh_plus = q.function_space().mesh()
-
-    # Extract indicator in enriched space
-    solver_obj = config.Solver(mesh_plus, q)
-    F = solver_obj.form
-    V = solver_obj.function_space
-    dwr_plus = get_dwr_indicator(F, q_star, test_space=V)
-
-    # Project down to base space
-    P0 = FunctionSpace(mesh, "DG", 0)
-    dwr = project(dwr_plus, P0)
-    dwr.interpolate(abs(dwr))
-    return dwr
diff --git a/nn_adapt/utility.py b/nn_adapt/utility.py
index 229de16..3668333 100644
--- a/nn_adapt/utility.py
+++ b/nn_adapt/utility.py
@@ -37,6 +37,15 @@ def _chk(self, val, old, rtol, reason):
                 self.converged_reason = reason
                 converged = True
         return converged
+    
+    def qoi_abs_chk(self, val, old, reason):
+        print("checking")
+        converged = False
+        if old is not None and self.fp_iteration >= self.miniter:
+            if abs(val - old) < 1e-8:
+                self.converged_reason = reason+"absolute error"
+                converged = True
+        return converged
 
     def check_qoi(self, val):
         """
@@ -44,6 +53,7 @@ def check_qoi(self, val):
         """
         r = "QoI convergence"
         converged = self._chk(val, self.qoi_old, self.qoi_rtol, r)
+        # converged = self.qoi_abs_chk(val, self.qoi_old, r)
         self.qoi_old = val
         return converged