diff --git a/benchmarks/cookies-problem-propagation/README.md b/benchmarks/cookies-problem-propagation/README.md
index 4f353059..07d75336 100644
--- a/benchmarks/cookies-problem-propagation/README.md
+++ b/benchmarks/cookies-problem-propagation/README.md
@@ -2,15 +2,16 @@
 
 ## Overview
 
-This benchmark runs a forward uncertainty quantification problem for the [cookies model](https://github.com/UM-Bridge/benchmarks/tree/main/models/cookies-problem/README.md) using the [Sparse Grids Matlab Kit](https://github.com/lorenzo-tamellini/sparse-grids-matlab-kit) interface to UM-Bridge. See below for full description.
+This benchmark runs a forward uncertainty quantification problem for the [cookies model](https://github.com/UM-Bridge/benchmarks/tree/main/models/fenics-cookies-problem/README.md) using the [Sparse Grids Matlab Kit](https://github.com/lorenzo-tamellini/sparse-grids-matlab-kit) interface to UM-Bridge. See below for full description.
 
 ## Authors
+- [Benjamin Kent](kent@imati.cnr.it)
 - [Massimiliano Martinelli](mailto:martinelli@imati.cnr.it)
 - [Lorenzo Tamellini](mailto:tamellini@imati.cnr.it)
 
 ## Run
 ```
-docker run -it -p 4242:4242 linusseelinger/cookies-problem
+docker run -it -p 4242:4242 linusseelinger/cookiebenchmark
 ```
 
 ## Properties
@@ -24,7 +25,7 @@ benchmark | Sets the config options for the forward UQ benchmark (see below)
 Mapping | Dimensions   | Description
 ---     |---           |---
 input   | [8]          | These values modify the conductivity coefficient in the 8 cookies. They are i.i.d. uniform random variables in the range [-0.99, -0.2] (software does not check that inputs are within the bound) 
-output  | \[1\]        | The integral of the solution over the central subdomain (see definition of $\Psi$ at [cookies model](https://github.com/UM-Bridge/benchmarks/tree/main/models/cookies-problem/README.md) for info)
+output  | \[1\]        | The integral of the solution over the central subdomain (see definition of $$\Psi$$ at [cookies model](https://github.com/UM-Bridge/benchmarks/tree/main/models/fenics-cookies-problem/README.md) for info)
 
 Feature       | Supported
 ---           |---
@@ -33,11 +34,9 @@ Gradient      | False
 ApplyJacobian | False
 ApplyHessian  | False
 
-Config        | Type    | Default value   	| Can be changed in benchmark model	| Description
----           |---      |---      		|---		  			| ---	
-NumThreads    | integer | 1     		| yes					| Number of physical cores to be used by the solver 
-BasisDegree   | integer | 4       		| no					| Default degree of spline basis (must be a positive integer)
-Fidelity      | integer | 2       		| no					| Controls the number of mesh elements (must be a positive integer, see below for details)
+Config      | Type    | Default value | Description
+---         |---      |---      		| ---	
+None        |
 
 
 ## Mount directories
@@ -47,18 +46,25 @@ None            |
 
 ## Source code
 
-[Benchmark sources available at this folder.](https://github.com/UM-Bridge/benchmarks/tree/main/benchmarks/cookies-problem)
+[Benchmark sources available at this folder.](https://github.com/UM-Bridge/benchmarks/tree/main/benchmarks/cookies-problem-propagation)
 
 ## Description
 
-![cookies-problem](https://raw.githubusercontent.com/UM-Bridge/benchmarks/main/models/cookies-problem/cookies_domain.png "geometry of the cookies problem")
+![cookies-problem](https://raw.githubusercontent.com/UM-Bridge/benchmarks/main/models/fenics-cookies-problem/cookies_domain.png "geometry of the cookies problem")
 
-The benchmark implements a forward uncertainty quantification problem for the [cookies model](https://github.com/UM-Bridge/benchmarks/tree/main/models/cookies-problem/README.md). More specifically, we assume that the uncertain parameters $y_n$ appearing in the definition of the diffusion coefficient are uniform i.i.d. random variables on the range $[-0.99, -0.2]$ and we aim at computing the expected value of the quantity of interest (i.e., output of the model) $\Psi$, which is defined as the integral of the solution over $F$.
+The benchmark implements a forward uncertainty quantification problem for the elliptic version of the [cookies model](https://github.com/UM-Bridge/benchmarks/tree/main/models/fenics-cookies-problem/README.md). More specifically, we assume that the uncertain parameters $$y_n$$ appearing in the definition of the diffusion coefficient are uniform i.i.d. random variables on the range $$[-0.99, -0.2]$$ and we aim at computing the expected value of the quantity of interest (i.e., output of the model) $$\Psi$$, which is defined as the integral of the solution over $$F$$.
 
-The PDE is solved with an IGA solver that uses as basis splines of degree $p=4$ and maximal regularity, i.e. of continuity $3$, and the mesh has $200 \times 200$ elements (i.e., the fidelity config parameter is set to $2$). The structure of this benchmark is identical to the one discussed in [[Bäck et al.,2011]](https://doi.org/10.1007/978-3-642-15337-2_3); however, raw numbers are different since in [[Bäck et al.,2011]](https://doi.org/10.1007/978-3-642-15337-2_3) the PDE solver employed was different (standard FEM with piecewise linear basis) and the mesh was also different.
+The benchmark configuration of the docker uses all config options set to their default values, see againg the [cookies model page](https://github.com/UM-Bridge/benchmarks/tree/main/models/fenics-cookies-problem/README.md). The structure of this benchmark thus is identical to the one discussed in [[Bäck et al.,2011]](https://doi.org/10.1007/978-3-642-15337-2_3); however, raw numbers are different since in [[Bäck et al.,2011]](https://doi.org/10.1007/978-3-642-15337-2_3) a different mesh was used.
 
+As a reference value, we provide the approximation of the expected value computed with a standard Smolyak sparse grid, based on Clenshaw--Curtis points, for levels $$w=0,1,\ldots,5$$, see e.g. [[Piazzola et al.,2024]](https://doi.org/10.1145/3630023). 
 
-As a reference value, we provide the approximation of the expected value computed with a standard Smolyak sparse grid, based on Clenshaw--Curtis points, for level $w=5$, see e.g. [[Piazzola et al.,2023]](https://doi.org/10.48550/arXiv.2203.09314). The resulting sparse grid has 15713 points, and the corresponding approximation of the expected value is $0.064196096847169$.
+Sparse grid $$w$$ | number of collocation points    | Estimate of $$\Psi$$ 
+----------------  |-------------------------------- |-------------------
+0                 | 1                               | 0.062255257529767
+1                 | 17                              | 0.064176316082952
+2                 | 145                             | 0.064206407272061
+3                 | 849                             | 0.064202639076811
+4                 | 3937                            | 0.064202350667514
+5                 | 15713                           | 0.064202367186117
 
-
-The script available [here](https://github.com/UM-Bridge/benchmarks/tree/main/benchmarks/cookies-problem/run_forward_benchmark_in_matlab.m) generates the results, using the Sparse Grids Matlab Kit [[Piazzola et al.,2023]](https://doi.org/10.48550/arXiv.2203.09314) for generating sparse grids. The Grids Matlab Kit is available on Github [here](https://github.com/lorenzo-tamellini/sparse-grids-matlab-kit) and a dedicated website with full resources including user manual is available [here](https://sites.google.com/view/sparse-grids-kit).
+The script available [here](https://github.com/UM-Bridge/benchmarks/tree/main/benchmarks/cookies-problem-propagation/run_forward_benchmark_in_matlab.m) generates the results, using the Sparse Grids Matlab Kit [[Piazzola et al.,2024]](https://doi.org/10.1145/3630023) for generating sparse grids. The Grids Matlab Kit is available on Github [here](https://github.com/lorenzo-tamellini/sparse-grids-matlab-kit) and a dedicated website with full resources including user manual is available [here](https://sites.google.com/view/sparse-grids-kit).
diff --git a/benchmarks/cookies-problem-propagation/run_forward_benchmark_in_matlab.m b/benchmarks/cookies-problem-propagation/run_forward_benchmark_in_matlab.m
index 498351d6..348ea3a9 100644
--- a/benchmarks/cookies-problem-propagation/run_forward_benchmark_in_matlab.m
+++ b/benchmarks/cookies-problem-propagation/run_forward_benchmark_in_matlab.m
@@ -14,7 +14,8 @@
 % 2) add sparse grids matlab kit and um-bridge matlab client to path
 % addpath(genpath('~/GIT_projects/Github/sparse-grids-matlab-kit/')) % this is version 23.5 Robert
 % addpath(genpath('~/GIT_projects/Github/umbridge/matlab/'))
-
+addpath(genpath('../../../sparse-grids-matlab-kit/'))
+addpath(genpath('../../../umbridge'))
 
 % 3) to save results on file (sparse grids in .txt / .mat, results in .txt / .mat), set the flag below to true
 saving_stuff = true;
@@ -31,7 +32,7 @@
 model = HTTPModel(uri,'benchmark');
 
 % config cookie solver with num threades
-config = struct('NumThreads',4);
+config = struct();
 
 
 % wrap model in an @-function too
@@ -50,7 +51,7 @@
 knots = @(n) knots_CC(n,-0.99,-0.2);
 lev2knots = @lev2knots_doubling;
 idxset_rule = @(i) sum(i-1);
-idxset_level_max = 1; % <--- controls max size of sparse grid
+idxset_level_max = 5; % <--- controls max size of sparse grid
 
 
 % To recycle evaluations from one grid to the next, we need some containers
diff --git a/models/cookies-problem/Dockerfile b/models/cookies-problem/Dockerfile
deleted file mode 100644
index 8be924e9..00000000
--- a/models/cookies-problem/Dockerfile
+++ /dev/null
@@ -1,60 +0,0 @@
-FROM ubuntu:mantic
-
-RUN apt update
-RUN DEBIAN_FRONTEND=noninteractive apt install -y autoconf automake autoconf-archive bash python3-pip pbzip2 wget clang-17 git cmake curl zip pkg-config gfortran lld libomp5-17 libomp-17-dev
-RUN pip3 install --break-system-packages umbridge
-
-
-RUN wget -O- https://apt.repos.intel.com/intel-gpg-keys/GPG-PUB-KEY-INTEL-SW-PRODUCTS.PUB \ | gpg --dearmor | tee /usr/share/keyrings/oneapi-archive-keyring.gpg > /dev/null
-RUN echo "deb [signed-by=/usr/share/keyrings/oneapi-archive-keyring.gpg] https://apt.repos.intel.com/oneapi all main" | tee /etc/apt/sources.list.d/oneAPI.list
-
-RUN apt update
-RUN apt install -y intel-basekit
-# RUN . /opt/intel/oneapi/setvars.sh
-
-RUN git clone https://gitlab.com/max.martinelli/igatools.git && \
-  cd igatools && \
-  git submodule init && \
-  git submodule update && \
-  git checkout fc6135ca63ed9062eefedfec5fda19eb3e3f4769 && \
-  cd vcpkg && ./bootstrap-vcpkg.sh
-
-#  git checkout 576316cfcec2be3212cba7a95f2b2df60a3f64c2 && \
-#  git checkout a335daa3bf480cb6e57e2ea8522e3817e841a66d && \
-
-
-WORKDIR /build_igatools
-RUN . /opt/intel/oneapi/setvars.sh && cd /build_igatools && \
-  cmake /igatools \
-  -G"Unix Makefiles" \
-  -DIGATOOLS_COMPONENT_DOCUMENTATION=OFF \
-  -DIGATOOLS_WITH_GMP=ON \
-  -DIGATOOLS_WITH_MKL=ON \
-  -DMKL_DIR=/opt/intel/oneapi/mkl/latest \
-  -DIGATOOLS_WITH_QUADRUPLE_PRECISION=OFF \
-  -DIGATOOLS_WITH_SERIALIZATION=OFF \
-  -DIGATOOLS_WITH_SUITESPARSE=OFF \
-  -DIGATOOLS_WITH_SUPERLU=OFF \
-  -DIGATOOLS_WITH_TRILINOS=OFF \
-  -DTBB_DIR=/opt/intel/oneapi/tbb/latest \
-  -DTBB_LIBRARY_DIR=/opt/intel/oneapi/tbb/latest/lib/intel64/gcc4.8/ \
-  -DTRILINOS_DIR=/opt/trilinos/current \
-  -DIGATOOLS_WITH_THREADS=ON \
-  -Wno-dev \
-  -DCMAKE_INSTALL_PREFIX=igatools \
-  -DCMAKE_BUILD_TYPE=Release \
-  -DIGATOOLS_STATIC_EXECUTABLE=OFF \
-  -DCMAKE_C_COMPILER=clang-17 \
-  -DCMAKE_CXX_COMPILER=clang++-17 \
-  -DBUILD_SHARED_LIBS=ON
-
-RUN make -j4 && \
-  make setup_tests && \
-  cd tests/ && \
-  make poisson_lorenzo.release
-
-WORKDIR /
-
-COPY umbridge-server.py /
-
-CMD python3 umbridge-server.py
diff --git a/models/cookies-problem/README.md b/models/cookies-problem/README.md
deleted file mode 100644
index d40471ce..00000000
--- a/models/cookies-problem/README.md
+++ /dev/null
@@ -1,88 +0,0 @@
-# The cookies model
-
-## Overview
-This model implements the so-called 'cookies problem' or 'cookies in the oven problem' (see for reference [[Bäck et al.,2011]](https://doi.org/10.1007/978-3-642-15337-2_3), [[Ballani et al.,2015]](https://doi.org/10.1137/140960980), [[Kressner et al., 2011]](https://doi.org/10.1137/100799010)), i.e., a simplified thermal equation in which the conductivity coefficient is uncertain in 8 circular subdomains ('the cookies'), whereas it is known (and constant) in the remaining of the domain ('the oven'). The PDE is solved by an isogeometric solver with maximum continuity splines, whose degree can be set by the user. See below for full description. 
-
-
-## Authors
-- [Massimiliano Martinelli](mailto:martinelli@imati.cnr.it)
-- [Lorenzo Tamellini](mailto:tamellini@imati.cnr.it)
-
-## Run
-```
-docker run -it -p 4242:4242 linusseelinger/cookies-problem
-```
-
-## Properties
-
-Model     | Description
----       | ---
-forward   | Forward evaluation of the cookies model, all config options can be modified by the user (see below)
-benchmark | Sets the config options for the forward UQ benchmark [(see benchmark page)](https://github.com/UM-Bridge/benchmarks/tree/main/benchmarks/cookies-problem/README.md)
-
-### forward
-
-Mapping | Dimensions | Description
----     |---         |---
-input   | [8]        | These values modify the conductivity coefficient in the 8 cookies, each of them must be greater than -1 (software does not check that input values are valid)  
-output  | [1]        | The integral of the solution over the central subdomain (see definition of $\Psi$ below)
-
-Feature       | Supported
----           |---
-Evaluate      | True
-Gradient      | False
-ApplyJacobian | False
-ApplyHessian  | False
-
-Config        | Type    | Default | Description
----           |---      |---      |---
-NumThreads    | integer | 1       | Number of physical cores to be used by the solver
-BasisDegree   | integer | 4       | Default degree of spline basis (must be a positive integer)
-Fidelity      | integer | 2       | Controls the number of mesh elements (must be a positive integer, see below for details)
-
-
-## Mount directories
-Mount directory | Purpose
----             |---
-None            | 
-
-## Source code
-
-[Model sources here.](https://github.com/UM-Bridge/benchmarks/tree/main/models/cookies-problem)
-
-## Description
-
-![cookies-problem](https://raw.githubusercontent.com/UM-Bridge/benchmarks/main/models/cookies-problem/cookies_domain.png "geometry of the cookies problem")
-
-The model implements the version of the cookies problem in [[Bäck et al.,2011]](https://doi.org/10.1007/978-3-642-15337-2_3), see also e.g. [[Ballani et al.,2015]](https://doi.org/10.1137/140960980), [[Kressner et al., 2011]](https://doi.org/10.1137/100799010) for slightly different versions. With reference to the computational domain $D=[0,1]^2$ in the figure above, the cookies model consists in the thermal diffusion problem below, where $\mathbf{y}$ are the uncertain parameters discussed in the following and $\mathrm{x}$ are physical coordinates 
-
-$$-\mathrm{div}\Big[ a(\mathbf{x},\mathbf{y}) \nabla u(\mathbf{x},\mathbf{y}) \Big] = f(\mathrm{x}), \quad \mathbf{x}\in D$$
-
-with homogeneous Dirichlet boundary conditions and forcing term defined as
-
-$$f(\mathrm{x}) = \begin{cases} 
-100 &\text{if } \,  \mathrm{x} \in F \\
-0 &\text{otherwise} 
-\end{cases}$$
-
-where $F$ is the square $[0.4, 0.6]^2$. The 8 subdomains with uncertain diffusion coefficient (the cookies) are circles with radius 0.13 and the following center coordinates:
-
-cookie | 1   | 2   | 3   | 4   | 5   | 6   | 7   | 8   |
---     | --  | --  | --  | --  | --  | --  | --  | --  |
-x      | 0.2 | 0.5 | 0.8 | 0.2 | 0.8 | 0.2 | 0.5 | 0.8 |
-y      | 0.2 | 0.2 | 0.2 | 0.5 | 0.5 | 0.8 | 0.8 | 0.8 |
-
-The uncertain diffusion coefficient is defined as
-
-$$a = 1 + \sum_{n=1}^8 y_n \chi_n(\mathrm{x})$$ 
-
-
-where $y_n>-1$ and 
-
-$$\chi_n(\mathrm{x}) = \begin{cases} 1 &\text{inside the n-th cookie} \\ 0 &\text{otherwise} \end{cases}$$
-
-
-The output of the simulation is the integral of the solution over $F$, i.e. $\Psi = \int_F u(\mathrm{x}) d \mathrm{x}$
-
-
-The PDE is solved with an IGA solver (see e.g. [[Da Veiga et al., 2014]](https://doi.org/10.1017/S096249291400004X)) that uses as basis splines of degree $p$ (tunable by the user, default $p=4$) of maximal regularity, i.e. of continuity $p-1$. The computational mesh is an $N\times N$ quadrilateral mesh (cartesian product of knot lines) with square elements, with $N=100 \times \mathrm{Fidelity}$. The implementation is done using the C++ library IGATools [[Pauletti et al., 2015]](https://doi.org/10.1137/140955252), available at [gitlab.com/max.martinelli/igatools](https://gitlab.com/max.martinelli/igatools).
diff --git a/models/cookies-problem/draw_geometry.m b/models/cookies-problem/draw_geometry.m
deleted file mode 100644
index fdae0c90..00000000
--- a/models/cookies-problem/draw_geometry.m
+++ /dev/null
@@ -1,28 +0,0 @@
-clear
-
-% basic domain
-
-% circles
-xc = [0.2 0.5 0.8 0.2 0.8 0.2 0.5 0.8];
-yc = [0.2 0.2 0.2 0.5 0.5 0.8 0.8 0.8];
-N = length(xc);
-
-radius = 0.13;
-theta = linspace(0,2*pi,200);
-
-for n=1:N
-    plot(xc(n) + radius*cos(theta), yc(n) + radius*sin(theta),'-k','LineWidth',2);
-    hold on
-    text(xc(n)-0.05,yc(n)+0.01,strcat('# ',num2str(n)),'FontSize',15)
-end
-
-% central square
-xsquare = [0.4 0.6 0.6 0.4 0.4];
-ysquare = [0.4 0.4 0.6 0.6 0.4];
-plot(xsquare,ysquare,'-k','LineWidth',2);
-text(0.485,0.5,'F','FontSize',15)
-
-axis square
-
-saveas(gcf,'cookies_domain','png')
-
diff --git a/models/cookies-problem/umbridge-server.py b/models/cookies-problem/umbridge-server.py
deleted file mode 100644
index 927a665e..00000000
--- a/models/cookies-problem/umbridge-server.py
+++ /dev/null
@@ -1,81 +0,0 @@
-import umbridge
-import os
-
-class TestModel(umbridge.Model):
-
-    def __init__(self):
-        super().__init__("forward")
-
-    def get_input_sizes(self, config):
-        return [8]
-
-    def get_output_sizes(self, config):
-        return [1]
-
-    def __call__(self, parameters, config):
-        arguments = " ".join([str(x) for x in parameters[0]]);
-
-        num_threads = str(config.get("NumThreads",1));
-        basis_degree = str(config.get("BasisDegree",4));
-        fidelity = str(config.get("Fidelity",2));
-
-        arguments = arguments + " " + basis_degree + " " + fidelity
-
-        # System call, cd into working directory and call model binary
-        os.system('. /opt/intel/oneapi/setvars.sh && export IGATOOLS_NUM_THREADS=' + num_threads + ' && /build_igatools/tests/models/poisson_lorenzo/poisson_lorenzo.release ' + arguments)
-
-        # Read second line of output file
-        with open('/poisson_results.dat', 'r') as f:
-            line = f.readline() # Read first line
-
-        return [[float(line)]]
-
-
-    def supports_evaluate(self):
-        return True
-
-
-
-
-class TestBenchmark(umbridge.Model):
-
-    def __init__(self):
-        super().__init__("benchmark")
-
-    def get_input_sizes(self, config):
-        return [8]
-
-    def get_output_sizes(self, config):
-        return [1]
-
-    def __call__(self, parameters, config):
-        arguments = " ".join([str(x) for x in parameters[0]]);
-
-        num_threads = str(config.get("NumThreads",1));
-        basis_degree = str(4);
-        fidelity = str(2);
-
-        arguments = arguments + " " + basis_degree + " " + fidelity
-
-        # System call, cd into working directory and call model binary
-        os.system('. /opt/intel/oneapi/setvars.sh && export IGATOOLS_NUM_THREADS=' + num_threads + ' && /build_igatools/tests/models/poisson_lorenzo/poisson_lorenzo.release ' + arguments)
-
-        # Read second line of output file
-        with open('/poisson_results.dat', 'r') as f:
-            line = f.readline() # Read first line
-
-        return [[float(line)]]
-
-
-    def supports_evaluate(self):
-        return True
-
-
-
-
-
-testmodel = TestModel()
-testbenchmark = TestBenchmark()
-
-
-umbridge.serve_models([testmodel,testbenchmark], 4242)
diff --git a/models/fenics-cookies-problem/Dockerfile b/models/fenics-cookies-problem/Dockerfile
new file mode 100644
index 00000000..1839deec
--- /dev/null
+++ b/models/fenics-cookies-problem/Dockerfile
@@ -0,0 +1,19 @@
+FROM ubuntu:mantic
+
+# Install fenics and python packages
+RUN apt update
+RUN DEBIAN_FRONTEND=noninteractive apt install -y bash fenics python3-pip
+RUN pip3 install --break-system-packages umbridge
+
+# Expose port for UM-BRIDGE server
+EXPOSE 4242
+
+RUN apt update
+
+# Set up internal folders and files
+WORKDIR /
+COPY umbridge-server.py /
+COPY cookiepde.py /
+
+# Start Python UM-BRDIGE server
+CMD python3 umbridge-server.py
\ No newline at end of file
diff --git a/models/fenics-cookies-problem/README.md b/models/fenics-cookies-problem/README.md
new file mode 100644
index 00000000..e1c25091
--- /dev/null
+++ b/models/fenics-cookies-problem/README.md
@@ -0,0 +1,132 @@
+# The Cookies model
+
+## Overview
+This model implements the so-called 'cookies problem' or 'cookies in the oven problem' (see for reference [[Bäck et al.,2011]](https://doi.org/10.1007/978-3-642-15337-2_3), [[Ballani et al.,2015]](https://doi.org/10.1137/140960980), [[Kressner et al., 2011]](https://doi.org/10.1137/100799010)), i.e., a simplified thermal equation in which the conductivity coefficient is uncertain in 8 circular subdomains ('the cookies'), whereas it is known (and constant) in the remaining of the domain ('the oven'). The equation can be either stationary (elliptic PDE) or time-dependent (parabolic PDE). The spatial approximation is constructed using the legacy version of [FEniCs](https://fenicsproject.org/) via the Python interface, using quadrilateral meshes (FEM degree configurable by the user); the time-advancing scheme is adaptive with local error control, specifically an implementation of the TR-AB2 scheme, see, e.g., [[Iserles, 2008]](https://doi.org/10.1017/cbo9780511995569.009). See below for full description.
+
+## Authors
+- [Benjamin Kent](mailto:kent@imati.cnr.it)
+- [Massimiliano Martinelli](mailto:martinelli@imati.cnr.it)
+- [Lorenzo Tamellini](mailto:tamellini@imati.cnr.it)
+
+## Run 
+```
+docker run -p 4242:4242 -it linusseelinger/cookies-problem:latest
+```
+
+## Properties
+The Docker container serves three models, two for the elliptic version of the PDE and one for the parabolic one.
+
+Model     | Description
+---       | ---
+forward   | Forward evaluation of the elliptic version of the cookies model, all config options can be modified by the user (see below)
+forwardparabolic   | Forward evaluation of the parabolic version of the cookies model, all config options can be modified by the user (see below)
+benchmark | Sets the config options for the forward UQ benchmark [(see benchmark page)](https://github.com/UM-Bridge/benchmarks/tree/main/benchmarks/cookies-problem-propagation/README.md)
+
+### forward
+Mapping | Dimensions	| Description
+--------|-------------|------------
+input |	[8] |	These values modify the conductivity coefficient in the 8 cookies, each of them must be greater than -1 (software does not check that input values are valid)
+output |	[1] |	The integral of the solution over the central subdomain $F$ (see below for details)
+
+Feature	| Supported
+--------|---------
+Evaluate|	True
+Gradient|	False
+ApplyJacobian|	False
+ApplyHessian|	False
+
+Config          | Type      | Default      | Description
+--------------- |---------- | ----------| -------------
+ N              | integer   | 400       | The number of cells in each dimension (i.e. a mesh of $N^2$ elements).
+ Fidelity       | integer   | 4         | If provided, the mesh will be generated with N = 100 x Fidelity.  The 'fidelity' config key takes precedent over 'N' if both defined.
+BasisDegree     | Integer   | 1         | The degree of the piecewise polynomial FE approximation
+advection       | Integer   | 0         | A flag to add an an advection field to the problem. If set to 1 an advection term with advection field (see below) is added to the PDE problem.
+quad_degree     | Integer   | 8         | The quadrature degree used to evaluate integrals in the matrix assembly.
+diffzero        | Float     | 1.0       | Defines the background diffusion field $$a_0$$ (see below)
+directsolver    | Integer   | 1         | Uses a direct LU solve for linear system. Set to 0 to use GM-RES.
+pc              | String    | 'none'    | Preconditioning for the GM-RES solver, can be either 'ILU' or 'JACOBI'
+tol             | Float     | 1e-4      | Relative tolerance for the GM-RES solver.
+
+### forwardparabolic
+Mapping | Dimensions	| Description
+--------|-------------|------------
+input |	[8] |	These values modify the conductivity coefficient in the 8 cookies, each of them must be greater than -1 (software does not check that input values are valid)
+output |	[1] |	The integral of the solution over the central subdomain $F$ (see below for details)
+
+Feature	| Supported
+--------|---------
+Evaluate|	True
+Gradient|	False
+ApplyJacobian|	False
+ApplyHessian|	False
+
+
+Config          | Type      | Default      | Description
+--------------- |---------- | ----------| -------------
+ N              | integer   | 400       | The number of cells in each dimension (i.e. a mesh of $N^2$ elements).
+ Fidelity       | integer   | 4         | If provided, the mesh will be generated with N = 100 x Fidelity.  The 'fidelity' config key takes precedent over 'N' if both defined.
+BasisDegree     | Integer   | 1         | The degree of the piecewise polynomial FE approximation
+advection       | Integer   | 0         | A flag to add an an advection field to the problem. If set to 1 an advection term with advection field (see below) is added to the PDE problem.
+quad_degree     | Integer   | 8         | The quadrature degree used to evaluate integrals in the matrix assembly.
+diffzero        | Float     | 1.0       | Defines the background diffusion field $$a_0$$ (see below)
+letol           | Float     | 1e-4     | Local timestepping error tolerance for a simple implementation of TR-AB2 timestepping.
+T               | Float      | 10.0     | Final time for timestepping approximation. The QoI is evaluated and returned for time T.
+
+
+
+## Mount directories
+Mount directory | Purpose
+---             |---
+None            | 
+
+## Source code
+
+[Model sources here.](https://github.com/UM-Bridge/benchmarks/tree/main/models/fenics-cookies-problem)
+
+## Description
+
+![cookies-problem](https://raw.githubusercontent.com/UM-Bridge/benchmarks/main/models/fenics-cookies-problem/cookies_domain.png "geometry of the cookies problem")
+
+In its elliptic variant, the model implements the version of the cookies problem in [[Bäck et al.,2011]](https://doi.org/10.1007/978-3-642-15337-2_3), see also e.g. [[Ballani et al.,2015]](https://doi.org/10.1137/140960980), [[Kressner et al., 2011]](https://doi.org/10.1137/100799010) for slightly different versions. With reference to the computational domain $$D=[0,1]^2$$ in the figure above, the cookies model consists in the thermal diffusion problem below, where $$\mathbf{y}$$ are the uncertain parameters discussed in the following and $$\mathrm{x}$$ are physical coordinates.  
+
+$$-\mathrm{div}\Big[ a(\mathbf{x},\mathbf{y}) \nabla u(\mathbf{x},\mathbf{y}) \Big] = f(\mathrm{x}), \quad \mathbf{x}\in D$$
+
+with homogeneous Dirichlet boundary conditions and forcing term defined as
+
+$$f(\mathrm{x}) = \begin{cases} 
+100 &\text{if } \,  \mathrm{x} \in F \\
+0 &\text{otherwise} 
+\end{cases}$$
+
+where $$F$$ is the square $$[0.4, 0.6]^2$$. The 8 subdomains with uncertain diffusion coefficient (the cookies) are circles with radius 0.13 and the following center coordinates:
+
+cookie | 1   | 2   | 3   | 4   | 5   | 6   | 7   | 8   |
+--     | --  | --  | --  | --  | --  | --  | --  | --  |
+x      | 0.2 | 0.5 | 0.8 | 0.2 | 0.8 | 0.2 | 0.5 | 0.8 |
+y      | 0.2 | 0.2 | 0.2 | 0.5 | 0.5 | 0.8 | 0.8 | 0.8 |
+
+The uncertain diffusion coefficient is defined as
+
+$$a = a_0 + \sum_{n=1}^8 y_n \chi_n(\mathrm{x})$$ 
+
+
+where $$a_0=1$$ by default (can be changed by the user), $$y_n>-1$$ and 
+
+$$\chi_n(\mathrm{x}) = \begin{cases} 1 &\text{inside the n-th cookie} \\ 0 &\text{otherwise} \end{cases}$$
+
+
+The output of the simulation is the integral of the solution over $$F$$, i.e. $$\Psi = \int_F u(\mathrm{x}) d \mathrm{x}$$
+
+An advection term can be added to the equation, by suitably setting the config options. In this case, the equation becomes
+
+$$-\mathrm{div}\Big[ a(\mathbf{x},\mathbf{y}) \nabla u(\mathbf{x},\mathbf{y}) \Big] + \mathbf{w}(\mathbf{x}) \cdot \nabla u(\mathbf{x},\mathbf{y}) = f(\mathrm{x}), \quad \mathbf{x}\in D$$
+
+with $$\mathbf{w}(\mathbf{x})=  [4(x_2-0.5)(1-4(x_1-0.5)^2), \,\, -4(x_1-0.5)*(1-4(x_2-0.5)^2)].$$
+
+The PDE is solved by a classical Finite Eement Method (using the legacy version of [FEniCs](https://fenicsproject.org/) via the Python interface), with standard Lagrangian polynomial bases of degree $$p$$ ove quadrilateral meshes (FEM degree configurable by the user). The meshes contain $$100\times K$$ elements per direction. The corresponding linear system can be solvers either directly (LU solver) or iteratively by GMRES method, with three choices of preconditioning (none, ILU, Jacobi), up to the users. The relative tolerance of the iterative method can also be set by the user.
+
+A parabolic (time-dependent) variant of the same problem is provided, that can be selected by using the forwardparabolic model introduced above. In this case, the equation becomes
+
+$$\displaystyle \frac{\partial u(\mathbf{x},t,\mathbf{y})}{\partial t} -\mathrm{div}\Big[ a(\mathbf{x},\mathbf{y}) \nabla u(\mathbf{x},\mathbf{y}) \Big] + \mathbf{w}(\mathbf{x}) \cdot \nabla u(\mathbf{x},t,\mathbf{y}) = f(\mathrm{x}), \quad \mathbf{x}\in D, t \in (0,T]$$
+
+with $$T$$ specified by the user and initial condition $$u=0$$. Like in the elliptic case, the advection term is null by default. All config options of the elliptic variant can be used here. Concerning time-stepping, the time-advancing scheme is adaptive with local error control, specifically an implementation of the TR-AB2 scheme, see, e.g., [[Iserles, 2008]](https://doi.org/10.1017/cbo9780511995569.009). 
diff --git a/models/fenics-cookies-problem/cookiepde.py b/models/fenics-cookies-problem/cookiepde.py
new file mode 100644
index 00000000..1354ca9b
--- /dev/null
+++ b/models/fenics-cookies-problem/cookiepde.py
@@ -0,0 +1,503 @@
+from dolfin import *
+import logging
+import sys
+import gc
+import numpy as np
+from petsc4py import PETSc
+
+
+class CookiePDE:
+    """
+    Class representing an elliptic or parabolic partial differential equation with ``cookie'' diffusion coefficient.
+    """
+
+    def __init__(self, N, basis_degree):
+        """
+        Initialize the PDE problem.
+
+        Parameters:
+        N (int): Number of mesh divisions in each direction.
+        """
+        sh = logging.StreamHandler(sys.stdout)
+        logging.basicConfig(level=logging.DEBUG, handlers=[sh])
+        set_log_active(True)
+        set_log_level(LogLevel.DEBUG)
+
+        # Create FEniCS mesh and define function space
+        # Define the corners of the rectangular domain
+        x0, y0 = 0.0, 0.0  # Bottom-left corner
+        x1, y1 = 1.0, 1.0  # Top-right corner
+
+        # Define the number of elements in each direction
+        nx, ny = N, N
+
+        p0 = Point(x0,y0)
+        p1 = Point(x1,y1)
+        # Create a rectangular mesh with quadrilateral elements
+        mesh = RectangleMesh.create([p0,p1],[nx,ny],CellType.Type.quadrilateral)
+
+        # Can alternatively use a triangular mesh
+        # mesh = UnitSquareMesh(N, N)
+
+        # Define approximation space
+        V = FunctionSpace(mesh, "Lagrange", basis_degree)
+        self.mesh = mesh
+        self.V = V
+
+    def setupProblem(self, difftype, y, quad_degree=8, varcoeffs=None, advection=0):
+        """
+        Set up the variational problem for the PDE.
+
+        Parameters:
+        difftype (str): Type of diffusion field (only 'cookie' implemented).
+        y (array): Array of coefficients for the diffusion field.
+        quad_degree (int, optional): Quadrature degree for integration (default is 8).
+        varcoeffs (array, optional): Array of variable coefficients (default is None).
+        advection (integer, optional): Whether to include advection term (default is 0 (false)).
+        """
+        V = self.V
+
+        # Define boundary condition
+        u0 = Constant(0.0)
+
+        # Homogeneous Dirichlet BC for unit square domain
+        def boundary(x):
+            return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS or x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS
+
+        bc = DirichletBC(V, u0, boundary)
+
+        # Define variational problem
+        u = TrialFunction(V)
+        v = TestFunction(V)
+
+        # Forcing on subdomain F defined by f_indicator
+        # F is an indicator function for the square subdomain [0.4,0.6]^2
+        # NOTE: I think using SpatialCoordinate avoids interpolating the discontinuity into the function space.
+        x = SpatialCoordinate(self.mesh)
+        f_indicator = conditional(lt(abs(x[0] - 0.5),0.1), 1.0, 0.0)*conditional(lt(abs(x[1] - 0.5),0.1), 1.0, 0.0)
+        # Define f using the indicator function
+        f = 100 * f_indicator
+
+        # Diffusion field defined via indicators on appropriate subdomains
+        # See the definition of the benchmark problem (https://doi.org/10.48550/arXiv.2402.13768)
+        if difftype == 'cookie':
+            # Input varcoeffs allows the background diffusion field to be changed
+            if varcoeffs is None:
+                a0 = 1.0
+            else:
+                a0 = varcoeffs[0]
+
+            # Cookie centres
+            cxlist = [0.2, 0.5, 0.8, 0.2, 0.8, 0.2, 0.5, 0.8]
+            cylist = [0.2, 0.2, 0.2, 0.5, 0.5, 0.8, 0.8, 0.8]
+
+            # Cookie radius
+            r = 0.13
+            diff = a0
+            for ii in range(8):
+                cx_ii = cxlist[ii]
+                cy_ii = cylist[ii]
+                distance = sqrt(pow(x[0] - cx_ii, 2) + pow(x[1] - cy_ii, 2))
+                indicator_ii = conditional(lt(distance, r), 1.0, 0.0)
+                diff = diff + indicator_ii * Constant(1.0-a0 + y[ii])
+        else:
+            error("Unknown diffusion field")
+
+        # Diffusion bilinear form
+        dx_q = Measure('dx', domain=self.mesh, metadata={'quadrature_degree': quad_degree})
+        a = inner(diff*grad(u), grad(v))*dx_q
+
+        # Optional: additional advection term. This is not included in the benchmark
+        if advection==1:
+            w = as_vector((4*(x[1]-0.5)*(1-4*pow(x[0]-0.5,2)), -4*(x[0]-0.5)*(1-4*pow(x[1]-0.5,2))))
+            a = a + inner(w, grad(u)) * v * dx_q
+
+        # Define linear form on RHS
+        L = f*v*dx_q
+
+        # Bilinear form for mass matrix
+        m = inner(u, v)*dx_q
+
+        # Generic function in approx space V
+        u = Function(V)
+
+        # Write to object
+        # self.diffproject = project(diff, V)
+        # self.fproject = project(f, V)
+        self.f_indicator = f_indicator
+        self.a = a
+        self.L = L
+        self.bc = bc
+        self.m = m
+        self.u = u
+        self.y = y
+
+    def solve(self, directsolver, pctype, tol):
+        """
+        Solve the elliptic PDE problem using PETSc.
+
+        Parameters:
+        directsolver (int): Use direct LU solver if equal to 1, otherwise GM-RES iterative solver.
+        pctype (str): Type of preconditioner to use with iterative solver ('ILU', 'JACOBI', or 'none').
+        tol (float): Tolerance for the iterative solver (if used).
+
+        Returns:
+        numpy.ndarray: Solution vector.
+        """
+        # Assemble FEniCS matrices
+        A = assemble(self.a)
+        b = assemble(self.L)
+        self.bc.apply(A)
+        self.bc.apply(b)
+
+        print("Solving for y=" + str(self.y))
+
+        # Convert from fenics to petsc
+        A_petsc = as_backend_type(A).mat()
+        b_petsc = as_backend_type(b).vec()
+        u_petsc = as_backend_type(self.u.vector()).vec()
+
+        # Construct solver
+        ksp = PETSc.KSP().create()
+        pc = PETSc.PC().create()
+
+        if directsolver == 1:
+            # Direct solver is full LU preconditioner (and hence no iterations required)
+            pc.setType(PETSc.PC.Type.LU)  # Set the preconditioner type
+            ksp.setType(PETSc.KSP.Type.PREONLY)  # Choose a solver type
+            pc.setOperators(A_petsc)  # Attach the matrix to the preconditioner
+            ksp.setPC(pc)
+        else:
+            # Use GMRES to approximate solution of linear system
+            ksp.setType(PETSc.KSP.Type.GMRES)  # Choose a solver type
+            ksp.setTolerances(atol=tol)  # Set tolerance
+
+            # Preconditioning
+            if pctype == "ILU":
+                pc.setType(PETSc.PC.Type.ILU)  # Set the preconditioner type
+            elif pctype == "JACOBI":
+                pc.setType(PETSc.PC.Type.JACOBI)  # Set the preconditioner type
+            else:
+                pass
+            pc.setOperators(A_petsc)  # Attach the matrix to the 
+            ksp.setPC(pc)
+        # Set up solver operator
+        ksp.setOperators(A_petsc)  # Set the operator
+        ksp.setFromOptions()  # Set PETSc options for the solver
+
+
+        # Solve the linear system
+        ksp.solve(b_petsc, u_petsc)  # Solve A*x = b for x
+        iternum = ksp.getIterationNumber()
+        print("Solved in: "+str(iternum)+" iterations")
+
+        # Write vector to FEniCS vector
+        self.u.vector()[:] = u_petsc[:]
+
+        # Destroy petsc objects
+        A_petsc.destroy()
+        b_petsc.destroy()
+        u_petsc.destroy()
+        ksp.destroy()
+        pc.destroy()
+
+        # Return dof as vector
+        return self.u.vector().get_local()
+
+    def projectref(self, n_ref):
+        """
+        Interpolate the solution to a finer mesh.
+
+        Parameters:
+        n_ref (int): Number of mesh divisions in each direction for the finer mesh.
+
+        Returns:
+        numpy.ndarray: Interpolated solution vector on the finer mesh.
+        """
+        # Define the reference mesh
+        N = n_ref
+        mesh = UnitSquareMesh(N, N)
+        V = FunctionSpace(mesh, "Lagrange", 1)
+        # Interpolate on mesh
+        u_int = interpolate(self.u, V)
+        # Return values at nodes
+        return u_int.vector().get_local()
+
+    def computebenchmarkqoi(self):
+        """
+        Compute the benchmark quantity of interest.
+
+        Returns:
+        float: Integral of the solution over the region F.
+        """
+        integral = assemble(self.f_indicator*self.u*dx)
+        return integral
+
+    def computenorm(self, x, normtype):
+        """
+        Compute the norm of the solution.
+
+        Parameters:
+        x (array): Solution vector.
+        normtype (str): Type of norm to compute.
+
+        Returns:
+        float: Computed norm of the solution.
+        """
+        u = Function(self.V)
+        u.vector().set_local(x[:])
+        return norm(u, normtype)
+
+    def writeSln(self, filename):
+        """
+        Write the solution to a Paraview file.
+
+        Parameters:
+        filename (str): Base name for the output file.
+        """
+        fileout = File(filename + ".pvd")
+        fileout << self.u
+        # fileout << self.diffproject
+        # fileout << self.fproject
+
+
+    def solveTime(self, tol, finalTime):
+        """
+        Solves the time-dependent problem using PETSc's TS solver.
+
+        Args:
+            tol (float): The relative tolerance for the solver.
+            finalTime (float): The final time up to which the solution is computed.
+
+        Returns:
+            numpy.ndarray: The local part of the solution vector.
+        """
+        # Compute solution
+        A = assemble(self.a)
+        b = assemble(self.L)
+        M = assemble(self.m)
+        self.bc.apply(A)
+        self.bc.apply(b)
+        self.bc.apply(M)
+
+        print("Solving for y=" + str(self.y))
+
+        # Set up PETSc formulation
+        A_petsc = as_backend_type(A).mat()
+        M_petsc = as_backend_type(M).mat()
+        b_petsc = as_backend_type(b).vec()
+        u_petsc = as_backend_type(self.u.vector()).vec()
+
+        f = u_petsc.copy()
+        fm = u_petsc.copy()
+
+        ts = PETSc.TS().create()
+
+        def rhs_function(ts, t, u, F, A, b):
+            A.mult(u, F)
+            F.axpby(1.0, -1.0, b)
+            return
+
+        def rhs_function_specific(ts, t, u, F):
+            rhs_function(ts, t, u, F, A_petsc, b_petsc)
+
+        ts.setRHSFunction(rhs_function_specific, f)
+
+        def mass_matrix_multiply(ts, t, u, udot, F, M):
+            M.mult(udot, F)  # F = M * u
+
+        def mass_matrix_multiply_specific(ts, t, u, udot, F):
+            mass_matrix_multiply(ts, t, u, udot, F, M_petsc)
+
+        ts.setIFunction(mass_matrix_multiply_specific, fm)
+
+        vtkfile = File("output.pvd")  # Output file name (with .pvd extension)
+
+        def monitor(ts, step, t, u):
+            print(f"Time Step {step}: Time = {t}")
+            self.u.vector()[:] = u[:]
+            self.u.rename("u", "label")
+            vtkfile << (self.u, t)  # Write function u to VTK file
+
+        # Set the monitor function for the TS solver
+        ts.setMonitor(monitor)
+
+        # Create a PETSc options object
+        opts = PETSc.Options()
+
+        # Set the adaptive basis type (e.g., default is "basic")
+        opts.setValue('-ts_adapt_type', 'basic')
+        # Other adaptive basis options can be set here
+
+        # Set solver options for adaptive time-stepping
+
+        ts.setProblemType(ts.ProblemType.LINEAR)
+        # ts.setEquationType(ts.EquationType.IMPLICIT)
+        ts.setType(ts.Type.BEULER)
+
+        ts.setSolution(u_petsc)
+        ts.setMaxTime(float(finalTime))
+
+        ts.setTolerances(rtol=tol)
+        ts.setTimeStep(1e-8)
+        ts.setFromOptions()
+
+        # ts.setTimeStep(1e-8)
+        ts.solve(u_petsc)
+
+        PETScOptions.clear()
+        ts.reset()
+        ts.destroy()
+
+        self.u.vector()[:] = u_petsc[:]
+
+        # Destroy PETSC objects
+        A_petsc.destroy()
+        M_petsc.destroy()
+        b_petsc.destroy()
+        u_petsc.destroy()
+        f.destroy()
+        fm.destroy()
+
+        gc.collect()
+        # PETSc.Log.view()
+
+        # Return solution vector at time T
+        return self.u.vector().get_local()
+
+
+    def solveTimeSimple(self, tol, finalTime):
+        """
+        Solves the time-dependent problem using a simple time-stepping approach.
+
+        Args:
+            tol (float): The relative tolerance for the solver.
+            finalTime (float): The final time up to which the solution is computed.
+
+        Returns:
+            numpy.ndarray: The local part of the solution vector.
+        """
+        # PETSc.Options().setValue('-log_view', '')
+        # PETSc.Options().setValue('-malloc_view', '')
+        # PETSc.Log.begin()
+
+        # Compute solution
+        A = assemble(self.a)
+        b = assemble(self.L)
+        M = assemble(self.m)
+        self.bc.apply(A)
+        self.bc.apply(b)
+        self.bc.apply(M)
+
+        print("Solving for y=" + str(self.y))
+
+        A_petsc = as_backend_type(A).mat()
+        M_petsc = as_backend_type(M).mat()
+        b_petsc = as_backend_type(b).vec()
+        u_petsc = as_backend_type(self.u.vector()).vec()
+        A = None
+        b = None
+        M = None
+        gc.collect()
+
+        # Output file name (with .pvd extension)
+        vtkfile = File("output" + str(tol) + ".pvd")
+
+        def monitor(ts, step, t, u, u2, enorm):
+            print(f"Time = {t:.4g},\t dt = {ts:.4g},\t E = {enorm:.4g},\t acc = {acc:.4g}")
+            # self.u.vector()[:] = u[:]
+            # self.u.rename("u", "label")
+            # vtkfile << (self.u, t)  # Write function u to VTK file
+
+        ksp = PETSc.KSP().create()
+
+        # Initialise
+        t = 0.0
+        dt = 1e-5
+        dtold = 1e9
+
+        lhs = A_petsc.copy()
+        lhs_old = A_petsc.copy()
+        rhs = u_petsc.copy()
+        rhs_temp = u_petsc.copy()
+        abf2 = u_petsc.copy()
+        u_petsc_m1 = u_petsc.copy()
+        e_petsc = u_petsc.copy()
+        e_temp = u_petsc.copy()
+
+        ii = 0
+        # Timestepping loop.
+        # Timestep dt is adapted using TR-AB2 pair and the Milne device to estimate local error.
+        # Solution is returned at finalTime
+        while t < finalTime:
+            if t + dt > finalTime:
+                dt = finalTime - t
+
+            lhs.zeroEntries()
+            lhs_old.zeroEntries()
+            rhs.zeroEntries()
+
+            # lhs_old = M_petsc + 0.5*dt*A_petsc
+            lhs.axpy(0.5 * dt, A_petsc)
+            lhs.axpy(1.0, M_petsc)
+
+            # rhs = (M_petsc - 0.5*dt*A_petsc) * u_petsc + dt*b_petsc
+            M_petsc.mult(u_petsc, rhs)
+            A_petsc.mult(u_petsc, rhs_temp)
+            rhs.axpy(-0.5 * dt, rhs_temp)
+            rhs.axpy(dt, b_petsc)
+
+            ksp.setOperators(lhs)  # Set the operator
+
+            # abf2 = u_petsc + dt * (-1.5 * A_petsc * u_petsc + 0.5 * A_petsc * u_petsc_m1)
+            A_petsc.mult(u_petsc, abf2)
+            A_petsc.mult(u_petsc_m1, rhs_temp)
+            abf2.axpby(0.5 * (dt / dtold), -(1 + 0.5 * dt / dtold), rhs_temp)
+            abf2.axpby(1.0, dt, u_petsc)
+            abf2.axpy(dt, b_petsc)
+
+            u_petsc_m1.axpby(1.0, 0.0, u_petsc)
+
+            # Solve linear system for timestep
+            ksp.solve(rhs, u_petsc)
+
+            # print(u_petsc[:])
+
+            # Estimate local timestepping error
+            e_petsc.zeroEntries()
+            e_petsc.axpby(1.0, 0.0, abf2)
+            e_petsc.axpy(-1.0, u_petsc)
+            M_petsc.mult(e_petsc, e_temp)
+            enorm = (dt / (3 * (1 + dtold / dt))) * np.sqrt(e_petsc.dot(e_temp))
+            ii += 1
+            t += dt
+
+            # Compute timestep acceleration
+            acc = float(np.power(tol / enorm, 2 / 3))
+            if acc > 10:
+                acc = 10  # Limit to exponential growth
+            monitor(dt, ii, t, u_petsc, abf2, enorm)
+            dtold = dt
+            if ii != 1:
+                dt = dt * acc  # Don't adapt first step
+
+            ksp.reset()
+
+        # Clean up
+        ksp.destroy()
+
+        self.u.vector()[:] = u_petsc[:]
+
+        PETScOptions.clear()
+
+        A_petsc.destroy()
+        M_petsc.destroy()
+        b_petsc.destroy()
+        u_petsc.destroy()
+
+        lhs.destroy()
+        lhs_old.destroy()
+        rhs.destroy()
+
+        # Return approximation for finalTime T
+        return self.u.vector().get_local()
\ No newline at end of file
diff --git a/models/cookies-problem/cookies_domain.png b/models/fenics-cookies-problem/cookies_domain.png
similarity index 100%
rename from models/cookies-problem/cookies_domain.png
rename to models/fenics-cookies-problem/cookies_domain.png
diff --git a/models/fenics-cookies-problem/test_output.py b/models/fenics-cookies-problem/test_output.py
new file mode 100644
index 00000000..9ea3e610
--- /dev/null
+++ b/models/fenics-cookies-problem/test_output.py
@@ -0,0 +1,50 @@
+#!/usr/bin/env python3
+
+# first run the container as
+#
+# docker run -it -p 4242:4242 <image name>  
+#
+# then run this script as python3 test_output.py http://localhost:4242
+
+import argparse
+import umbridge
+import pytest
+
+parser = argparse.ArgumentParser(description='Minimal HTTP model demo.')
+parser.add_argument('url', metavar='url', type=str, help='the ULR on which the model is running, for example http://localhost:4242')
+args = parser.parse_args()
+print(f"Connecting to host URL {args.url}")
+
+# Set up a model by connecting to URL
+model = umbridge.HTTPModel(args.url, "forward")
+
+#test get methods
+output = model.get_input_sizes()
+print("get_input_sizes() returns "+str(output[0]))
+assert pytest.approx(output[0]) == 8, "get_input_sizes() returns wrong value"
+
+
+output = model.get_output_sizes()
+print("get_output_sizes() returns "+str(output[0]))
+assert pytest.approx(output[0]) == 1, "get input sizes returns wrong value"
+
+
+#test output for default config (p=1, fid=4)
+param = [[-2.3939786473373e-01, -8.6610045659126e-01, -2.1086275315687e-01, -9.2604304103162e-01, -6.0002531612112e-01, -5.5677423053456e-01, -7.7546408441658e-01, -7.6957620518706e-01]]
+output = model(param)
+print("model output (quantity of interest) for default config values = "+str(output[0][0]))
+assert  pytest.approx(output[0][0], abs=1e-12) == 0.06933717558839111, "Output not as expected"
+
+#test output for another config
+output = model(param,{"BasisDegree": 3, "Fidelity": 3})
+print("model output (quantity of interest) = "+str(output[0][0]))
+assert  pytest.approx(output[0][0], abs=1e-12) == 0.06937554567023371, "Output not as expected"
+
+
+#another test, this time for the benchmark version (i.e. p=1, fid=4 again, same as model default)
+model_B = umbridge.HTTPModel(args.url, "benchmark")
+
+param = [[-0.2,-0.2,-0.3,-0.4,-0.5,-0.6,-0.7,-0.8]]
+output = model_B(param)
+print("model output (quantity of interest) in benchmark configuration = "+str(output[0][0]))
+assert pytest.approx(output[0][0], abs=1e-12) == 0.05725576523730292, "Output not as expected"
diff --git a/models/cookies-problem/test_output.py b/models/fenics-cookies-problem/test_output_parabolic.py
similarity index 76%
rename from models/cookies-problem/test_output.py
rename to models/fenics-cookies-problem/test_output_parabolic.py
index 8fe6d9fa..ff681d00 100644
--- a/models/cookies-problem/test_output.py
+++ b/models/fenics-cookies-problem/test_output_parabolic.py
@@ -16,7 +16,7 @@
 print(f"Connecting to host URL {args.url}")
 
 # Set up a model by connecting to URL
-model = umbridge.HTTPModel(args.url, "forward")
+model = umbridge.HTTPModel(args.url, "forwardparabolic")
 
 #test get methods
 output = model.get_input_sizes()
@@ -33,19 +33,17 @@
 param = [[-2.3939786473373e-01, -8.6610045659126e-01, -2.1086275315687e-01, -9.2604304103162e-01, -6.0002531612112e-01, -5.5677423053456e-01, -7.7546408441658e-01, -7.6957620518706e-01]]
 output = model(param)
 print("model output (quantity of interest) for default config values = "+str(output[0][0]))
-assert pytest.approx(output[0][0]) == 0.0693282746248043, "Output not as expected"
+assert  pytest.approx(output[0][0], abs=1e-6) == 0.06933585905253055, "Output not as expected"
 
 #test output for another config
-output = model(param,{"NumThreads": 10, "BasisDegree": 3, "Fidelity": 3})
+output = model(param,{"BasisDegree": 3, "Fidelity": 2})
 print("model output (quantity of interest) = "+str(output[0][0]))
-assert pytest.approx(output[0][0]) == 0.06934748547844366, "Output not as expected"
-
+assert  pytest.approx(output[0][0], abs=1e-6) == 0.06936772917504516, "Output not as expected"
 
 #another test, this time for the benchmark version (i.e. p=4, fid=2 again, same as model default)
-model_B = umbridge.HTTPModel(args.url, "benchmark")
+model_B = umbridge.HTTPModel(args.url, "benchmarkparabolic")
 
 param = [[-0.2,-0.2,-0.3,-0.4,-0.5,-0.6,-0.7,-0.8]]
-output = model_B(param,{"NumThreads": 10})
+output = model_B(param)
 print("model output (quantity of interest) in benchmark configuration = "+str(output[0][0]))
-assert pytest.approx(output[0][0]) == 0.05725269745090122, "Output not as expected"
-
+assert pytest.approx(output[0][0], abs=1e-6) == 0.05725981992101194, "Output not as expected"
\ No newline at end of file
diff --git a/models/fenics-cookies-problem/umbridge-server.py b/models/fenics-cookies-problem/umbridge-server.py
new file mode 100644
index 00000000..41f912b1
--- /dev/null
+++ b/models/fenics-cookies-problem/umbridge-server.py
@@ -0,0 +1,308 @@
+import umbridge
+from cookiepde import CookiePDE
+from dolfin import *
+
+
+class CookieForward(umbridge.Model):
+    """
+    Class representing the forward cookie elliptic PDE model.
+    """
+
+    def __init__(self):
+        """
+        Initialize the CookieForward model.
+        """
+        super().__init__("forward")
+
+    def get_input_sizes(self, config):
+        """
+        Get the sizes of the input parameters.
+
+        Parameters:
+        config (dict): Configuration dictionary.
+
+        Returns:
+        list: List containing the size of the input parameter vector.
+        """
+        return [8]
+
+    def get_output_sizes(self, config):
+        """
+        Get the sizes of the output parameters.
+
+        Parameters:
+        config (dict): Configuration dictionary.
+
+        Returns:
+        list: List containing the size of the output parameter vector.
+        """
+        return [1]
+
+    def __call__(self, parameters, config):
+        """
+        Evaluate the CookieForward model.
+
+        Parameters:
+        parameters (list): List of input parameters for the model.
+        config (dict): Configuration dictionary.
+
+        Returns:
+        list: List containing the quantity of interest.
+        """
+        # Fill missing entries in config with default values
+        config = verifyConfig(config)
+        # Initialize PDE model
+        model = CookiePDE(config['N'], config['BasisDegree'])
+        # Set up cookie problem
+        model.setupProblem('cookie', parameters[0], config['quad_degree'], varcoeffs=config['diffzero'])
+        # Solve linear system with preconditioning pc and solver tolerance tol
+        model.solve(config['directsolver'], config['pc'], config['tol'])
+        # Compute quantity of interest (QoI) on solution
+        integral = model.computebenchmarkqoi()
+        # Return QoI
+        return [[integral]]
+
+    def supports_evaluate(self):
+        """
+        Check if the model supports evaluation.
+
+        Returns:
+        bool: True, indicating that the model supports evaluation.
+        """
+        return True
+
+    def supports_gradient(self):
+        """
+        Check if the model supports gradient computation.
+
+        Returns:
+        bool: False, indicating that the model does not support gradient computation.
+        """
+        return False
+
+
+class CookieBenchmark(umbridge.Model):
+    """
+    Class representing the Cookie Benchmark elliptic PDE model.
+    """
+
+    def __init__(self):
+        """
+        Initialize the CookieBenchmark model.
+        """
+        super().__init__("benchmark")
+
+    def get_input_sizes(self, config):
+        """
+        Get the sizes of the input parameters.
+
+        Parameters:
+        config (dict): Configuration dictionary.
+
+        Returns:
+        list: List containing the size of the input parameter vector.
+        """
+        return [8]
+
+    def get_output_sizes(self, config):
+        """
+        Get the sizes of the output parameters.
+
+        Parameters:
+        config (dict): Configuration dictionary.
+
+        Returns:
+        list: List containing the size of the output parameter vector.
+        """
+        return [1]
+
+    def __call__(self, parameters, config):
+        """
+        Evaluate the Cookie Benchmark model.
+
+        Parameters:
+        parameters (list): List of input parameters for the model.
+        config (dict): Configuration dictionary.
+
+        Returns:
+        list: List containing the quantity of interest.
+        """
+        # Use default values by ensuring dictionary is empty
+        config={}
+        config = verifyConfig(config)
+        # Initialize PDE model
+        model = CookiePDE(config['N'], config['BasisDegree'])
+        # Set up cookie problem
+        model.setupProblem(
+            'cookie', parameters[0], config['quad_degree'], varcoeffs=config['diffzero'])
+        # Solve linear system, possibly with preconditioning pc and solver tolerance tol
+        model.solve(config['directsolver'],config['pc'], config['tol'])
+        # Compute quantity of interest (QoI) on solution
+        integral = model.computebenchmarkqoi()
+        # Return QoI
+        return [[integral]]
+
+    def supports_evaluate(self):
+        """
+        Check if the model supports evaluation.
+
+        Returns:
+        bool: True, indicating that the model supports evaluation.
+        """
+        return True
+
+    def supports_gradient(self):
+        """
+        Check if the model supports gradient computation.
+
+        Returns:
+        bool: False, indicating that the model does not support gradient computation.
+        """
+        return False
+
+class CookieTime(umbridge.Model):
+    """
+    A model for parabolic formulation of the cookie model.
+
+    Inherits from umbridge.Model.
+
+    Attributes:
+        None
+    """
+
+    def __init__(self):
+        """
+        Initializes the CookieTime object.
+
+        Args:
+            None
+
+        Returns:
+            None
+        """
+        super().__init__("forwardparabolic")
+
+    def get_input_sizes(self, config):
+        """
+        Returns the input sizes for the model.
+
+        Args:
+            config (dict): A dictionary containing configuration parameters.
+
+        Returns:
+            list: A list containing the input sizes.
+        """
+        return [8]
+
+    def get_output_sizes(self, config):
+        """
+        Returns the output sizes for the model.
+
+        Args:
+            config (dict): A dictionary containing configuration parameters.
+
+        Returns:
+            list: A list containing the output sizes.
+        """
+        return [1]
+    
+    def __call__(self, parameters, config):
+        """
+        Evaluates the parabolic cookie model.
+
+        Args:
+            parameters (list): A list of parameters.
+            config (dict): A dictionary containing configuration parameters.
+
+        Returns:
+            list: A list containing the computed integral.
+        """
+        # Verfiy config and fills in empty keys.
+        config = verifyConfig(config)
+
+        # Set up discrete formulation.
+        model = CookiePDE(config['N'], config['BasisDegree'])
+        # Note that we have an additional advection term
+        model.setupProblem('cookie', parameters[0], varcoeffs=config['diffzero'], advection=config['advection'])
+        # Use the custom TR-AB2 solver. Optional solveTime function uses built in PETSC TS solver.
+        u = model.solveTimeSimple(config['letol'],config['T'])
+        # Compute QoI at finalTime T
+        integral = model.computebenchmarkqoi()
+        # model.writeSln("outputFinal")
+        return [[integral]]
+
+    def supports_evaluate(self):
+        """
+        Checks if the model supports evaluation.
+
+        Returns:
+            bool: True if the model supports evaluation, False otherwise.
+        """
+        return True
+
+    def supports_gradient(self):
+        """
+        Checks if the model supports gradient computation.
+
+        Returns:
+            bool: True if the model supports gradient computation, False otherwise.
+        """
+        return False
+
+
+
+def verifyConfig(config):
+    if config is None:
+        config = {}
+
+    # Use direct solver by default
+    if 'directsolver' not in config:
+        config['directsolver'] = 1
+
+    # Use 400x400 mesh by default
+    if 'Fidelity' not in config:
+        config['N'] = 400
+    else:
+        config['N'] = int(100 * config['Fidelity'])
+
+    # Use Q1 approximation by default
+    if 'BasisDegree' not in config:
+        config['BasisDegree'] = 1
+    
+    # Use degree 8 quadrature by default
+    if 'quad_degree' not in config:
+        config['quad_degree'] = 8
+    
+    # Use background diffusion 1.0 by default
+    if 'diffzero' not in config:
+        config['diffzero'] = [1.0]
+    
+    # Use no preconditioning by default
+    if 'pc' not in config:
+        config['pc'] = "none"
+
+    # Use GM-RES tol 1e-4 by default
+    if 'tol' not in config:
+        config['tol'] = 1e-4
+    
+    # Use local timestepping error tolerance 1e-4 by default
+    if 'letol' not in config:
+        config['letol'] = 1e-4
+    
+    # Use a finalTime T = 10.0 by default
+    if 'T' not in config:
+        config['T'] = 10.0
+
+    if 'advection' not in config:
+        config['advection'] = 0
+
+    print(config)
+    return config
+
+# Initialise UM-BRIDGE models
+cookieforward = CookieForward()
+cookiebenchmark= CookieBenchmark()
+cookietime = CookieTime()
+
+# Start UM-BRIDGE server
+umbridge.serve_models([cookieforward,cookiebenchmark,cookietime], 4242)