Repository: damar-wicaksono/uqtestfuns (language: python)

Issue:
Add the Two-Dimensional Cooling Coffee Cup Model from Tennøe et al. (2018)
The two-dimensional cooling coffee cup test function models the cooling of a coffee cup following Newton's law of cooling:

$$
\frac{d T(t)}{dt} = -k (T(t) - T_{\mathrm{amb}}),
$$

where $T$ is the temperature of the cup, $k$ is the cooling rate constant, and $T_{\mathrm{amb}}$ is the ambient temperature.

The function is two-dimensional with $k$ and $T_{\mathrm{amb}}$ are considered uncertain input modeled as uniform random variables
$\mathcal{U}(0.025, 0.075)$ and $\mathcal{U}(15, 25)$, respectively.

The solver `scipy.odeint` is used to solve the differential equation with initial temperature $T_0 = T(t = 0) = 95 [\mathrm{^{o}C}]$ evaluated
at $150$ time points between $0$ and $200$ $[\mathrm{s}]$.

The model is taken from one of the examples of Uncertainpy package (Tennøe et al. (2018)[^Tennøe]); the model also appears in the illustration of EasyVVUQ package (Richardson et al. (2020)[^Richardson]).

[^Tennøe]: S. Tennøe, G. Halnes, and G. T. Einevoll, “Uncertainpy: A Python Toolbox for Uncertainty Quantification and Sensitivity Analysis in Computational Neuroscience,” *Frontiers in Neuroinformatics*, vol. 12, p. 49, 2018, DOI: [10.3389/fninf.2018.00049](https://doi.org/10.3389/fninf.2018.00049).

[^Richardson]: R. A. Richardson, D. W. Wright, W. Edeling, V. Jancauskas, J. Lakhlili, and P. V. Coveney, “EasyVVUQ: A Library for Verification, Validation and Uncertainty Quantification in High Performance Computing,” *Journal of Open Research Software*, vol. 8, no. 1, p. 11, 2020, DOI: [10.5334/jors.303](https://doi.org/10.5334/jors.303).



Discussion / hints:

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4ed9b51..4a4a2dc 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,6 +9,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
+- The two-dimensional, time-dependent (vector-valued) cooling cup model
+  for metamodeling exercise.
 - The 8-dimensional robot arm function for metamodeling exercises.
 
 ## [0.5.0] - 2024-11-18
diff --git a/docs/_toc.yml b/docs/_toc.yml
index 028eb9a..342d705 100644
--- a/docs/_toc.yml
+++ b/docs/_toc.yml
@@ -58,6 +58,8 @@ parts:
             title: Cheng and Sandu (2010) 2D
           - file: test-functions/circular-pipe-crack
             title: Circular Pipe Crack
+          - file: test-functions/coffee-cup
+            title: Coffee Cup
           - file: test-functions/convex-fail-domain
             title: Convex Failure Domain
           - file: test-functions/damped-cosine
diff --git a/docs/fundamentals/metamodeling.md b/docs/fundamentals/metamodeling.md
index 605d21e..0a96b1a 100644
--- a/docs/fundamentals/metamodeling.md
+++ b/docs/fundamentals/metamodeling.md
@@ -18,43 +18,44 @@ kernelspec:
 The table below listed the available test functions typically used
 in the comparison of metamodeling approaches.
 
-|                                  Name                                  | Input Dimension |      Constructor       |
-|:----------------------------------------------------------------------:|:---------------:|:----------------------:|
-|                 {ref}`Ackley <test-functions:ackley>`                  |        M        |       `Ackley()`       |
-|  {ref}`Alemazkoor & Meidani (2018) 2D <test-functions:alemazkoor-2d>`  |        2        |    `Alemazkoor2D()`    |
-| {ref}`Alemazkoor & Meidani (2018) 20D <test-functions:alemazkoor-20d>` |       20        |   `Alemazkoor20D()`    |
-|               {ref}`Borehole <test-functions:borehole>`                |        8        |      `Borehole()`      |
-|       {ref}`Cheng and Sandu (2010) 2D <test-functions:cheng2d>`        |        2        |       `Cheng2D`        |
-|          {ref}`Damped Cosine <test-functions:damped-cosine>`           |        1        |    `DampedCosine()`    |
-|      {ref}`Damped Oscillator <test-functions:damped-oscillator>`       |        7        |  `DampedOscillator()`  |
-|                  {ref}`Flood <test-functions:flood>`                   |        8        |       `Flood()`        |
-|       {ref}`Forrester et al. (2008) <test-functions:forrester>`        |        1        |   `Forrester2008()`    |
-|             {ref}`(1st) Franke <test-functions:franke-1>`              |        2        |      `Franke1()`       |
-|             {ref}`(2nd) Franke <test-functions:franke-2>`              |        2        |      `Franke2()`       |
-|             {ref}`(3rd) Franke <test-functions:franke-3>`              |        2        |      `Franke3()`       |
-|             {ref}`(4th) Franke <test-functions:franke-4>`              |        2        |      `Franke4()`       |
-|             {ref}`(5th) Franke <test-functions:franke-5>`              |        2        |      `Franke5()`       |
-|             {ref}`(6th) Franke <test-functions:franke-6>`              |        2        |      `Franke6()`       |
-|           {ref}`Friedman (6D) <test-functions:friedman-6d>`            |        6        |     `Friedman6D()`     |
-|          {ref}`Friedman (10D) <test-functions:friedman-10d>`           |       10        |    `Friedman10D()`     |
-|     {ref}`Gramacy (2007) 1D Sine <test-functions:gramacy-1d-sine>`     |        1        |   `Gramacy1DSine()`    |
-| {ref}`Lim et al. (2002) Non-Polynomial <test-functions:lim-non-poly>`  |        2        |     `LimNonPoly()`     |
-|     {ref}`Lim et al. (2002) Polynomial <test-functions:lim-poly>`      |        2        |      `LimPoly()`       |
-|              {ref}`McLain S1 <test-functions:mclain-s1>`               |        2        |      `McLainS1()`      |
-|              {ref}`McLain S2 <test-functions:mclain-s2>`               |        2        |      `McLainS2()`      |
-|              {ref}`McLain S3 <test-functions:mclain-s3>`               |        2        |      `McLainS3()`      |
-|              {ref}`McLain S4 <test-functions:mclain-s4>`               |        2        |      `McLainS4()`      |
-|              {ref}`McLain S5 <test-functions:mclain-s5>`               |        2        |      `McLainS5()`      |
-|      {ref}`Oakley & O'Hagan (2002) 1D <test-functions:oakley-1d>`      |        1        |      `Oakley1D()`      |
-|            {ref}`OTL Circuit <test-functions:otl-circuit>`             |     6 / 20      |     `OTLCircuit()`     |
-|            {ref}`Piston Simulation <test-functions:piston>`            |     7 / 20      |       `Piston()`       |
-|              {ref}`Robot Arm <test-functions:robot-arm>`               |        8        |      `RobotArm()`      |
-|          {ref}`Solar Cell Model <test-functions:solar-cell>`           |        5        |     `SolarCell()`      |
-|                 {ref}`Sulfur <test-functions:sulfur>`                  |        9        |       `Sulfur()`       |
-|    {ref}`Undamped Oscillator <test-functions:undamped-oscillator>`     |        6        | `UndampedOscillator()` |
-|      {ref}`Webster et al. (1996) 2D <test-functions:webster-2d>`       |        2        |     `Webster2D()`      |
-|         {ref}`Welch et al. (1992) <test-functions:welch1992>`          |       20        |     `Welch1992()`      |
-|            {ref}`Wing Weight <test-functions:wing-weight>`             |       10        |     `WingWeight()`     |
+|                                  Name                                   | Input Dimension |       Constructor       |
+|:-----------------------------------------------------------------------:|:---------------:|:-----------------------:|
+|                  {ref}`Ackley <test-functions:ackley>`                  |        M        |       `Ackley()`        |
+|  {ref}`Alemazkoor & Meidani (2018) 2D <test-functions:alemazkoor-2d>`   |        2        |    `Alemazkoor2D()`     |
+| {ref}`Alemazkoor & Meidani (2018) 20D <test-functions:alemazkoor-20d>`  |       20        |    `Alemazkoor20D()`    |
+|                {ref}`Borehole <test-functions:borehole>`                |        8        |      `Borehole()`       |
+|        {ref}`Cheng and Sandu (2010) 2D <test-functions:cheng2d>`        |        2        |        `Cheng2D`        |
+|           {ref}`Coffee Cup Model <test-functions:coffee-cup>`           |        2        |      `CoffeeCup()`      |
+|           {ref}`Damped Cosine <test-functions:damped-cosine>`           |        1        |    `DampedCosine()`     |
+|       {ref}`Damped Oscillator <test-functions:damped-oscillator>`       |        7        |  `DampedOscillator()`   |
+|                   {ref}`Flood <test-functions:flood>`                   |        8        |        `Flood()`        |
+|        {ref}`Forrester et al. (2008) <test-functions:forrester>`        |        1        |    `Forrester2008()`    |
+|              {ref}`(1st) Franke <test-functions:franke-1>`              |        2        |       `Franke1()`       |
+|              {ref}`(2nd) Franke <test-functions:franke-2>`              |        2        |       `Franke2()`       |
+|              {ref}`(3rd) Franke <test-functions:franke-3>`              |        2        |       `Franke3()`       |
+|              {ref}`(4th) Franke <test-functions:franke-4>`              |        2        |       `Franke4()`       |
+|              {ref}`(5th) Franke <test-functions:franke-5>`              |        2        |       `Franke5()`       |
+|              {ref}`(6th) Franke <test-functions:franke-6>`              |        2        |       `Franke6()`       |
+|            {ref}`Friedman (6D) <test-functions:friedman-6d>`            |        6        |     `Friedman6D()`      |
+|           {ref}`Friedman (10D) <test-functions:friedman-10d>`           |       10        |     `Friedman10D()`     |
+|     {ref}`Gramacy (2007) 1D Sine <test-functions:gramacy-1d-sine>`      |        1        |    `Gramacy1DSine()`    |
+|  {ref}`Lim et al. (2002) Non-Polynomial <test-functions:lim-non-poly>`  |        2        |     `LimNonPoly()`      |
+|      {ref}`Lim et al. (2002) Polynomial <test-functions:lim-poly>`      |        2        |       `LimPoly()`       |
+|               {ref}`McLain S1 <test-functions:mclain-s1>`               |        2        |      `McLainS1()`       |
+|               {ref}`McLain S2 <test-functions:mclain-s2>`               |        2        |      `McLainS2()`       |
+|               {ref}`McLain S3 <test-functions:mclain-s3>`               |        2        |      `McLainS3()`       |
+|               {ref}`McLain S4 <test-functions:mclain-s4>`               |        2        |      `McLainS4()`       |
+|               {ref}`McLain S5 <test-functions:mclain-s5>`               |        2        |      `McLainS5()`       |
+|      {ref}`Oakley & O'Hagan (2002) 1D <test-functions:oakley-1d>`       |        1        |      `Oakley1D()`       |
+|             {ref}`OTL Circuit <test-functions:otl-circuit>`             |     6 / 20      |     `OTLCircuit()`      |
+|            {ref}`Piston Simulation <test-functions:piston>`             |     7 / 20      |       `Piston()`        |
+|               {ref}`Robot Arm <test-functions:robot-arm>`               |        8        |      `RobotArm()`       |
+|           {ref}`Solar Cell Model <test-functions:solar-cell>`           |        5        |      `SolarCell()`      |
+|                  {ref}`Sulfur <test-functions:sulfur>`                  |        9        |       `Sulfur()`        |
+|     {ref}`Undamped Oscillator <test-functions:undamped-oscillator>`     |        6        | `UndampedOscillator()`  |
+|       {ref}`Webster et al. (1996) 2D <test-functions:webster-2d>`       |        2        |      `Webster2D()`      |
+|          {ref}`Welch et al. (1992) <test-functions:welch1992>`          |       20        |      `Welch1992()`      |
+|             {ref}`Wing Weight <test-functions:wing-weight>`             |       10        |     `WingWeight()`      |
 
 In a Python terminal, you can list all the available functions relevant
 for metamodeling applications using ``list_functions()``
diff --git a/docs/references.bib b/docs/references.bib
index d97ef2e..33fea4b 100644
--- a/docs/references.bib
+++ b/docs/references.bib
@@ -953,4 +953,15 @@ @Article{An2001
   doi     = {10.1006/jcom.2001.0588},
 }
 
+@Article{Richardson2020,
+  author  = {Richardson, Robin A. and Wright, David W. and Edeling, Wouter and Jancauskas, Vytautas and Lakhlili, Jalal and Coveney, Peter V.},
+  journal = {Journal of Open Research Software},
+  title   = {{EasyVVUQ}: a library for verification, validation and uncertainty quantification in high performance computing},
+  year    = {2020},
+  number  = {1},
+  pages   = {11},
+  volume  = {8},
+  doi     = {10.5334/jors.303},
+}
+
 @Comment{jabref-meta: databaseType:bibtex;}


Test patch:
diff --git a/docs/test-functions/available.md b/docs/test-functions/available.md
index 7734c90..fe5f835 100644
--- a/docs/test-functions/available.md
+++ b/docs/test-functions/available.md
@@ -32,6 +32,7 @@ regardless of their typical applications.
 |           {ref}`Cantilever Beam (2D) <test-functions:cantilever-beam-2d>`           |        2        |       `CantileverBeam2D `       |
 |              {ref}`Cheng and Sandu (2010) 2D <test-functions:cheng2d>`              |        2        |           `Cheng2D `            |
 |           {ref}`Circular Pipe Crack <test-functions:circular-pipe-crack>`           |        2        |      `CircularPipeCrack()`      |
+|                 {ref}`Coffee Cup Model <test-functions:coffee-cup>`                 |        2        |          `CoffeeCup()`          |
 |          {ref}`Convex Failure Domain <test-functions:convex-fail-domain>`           |        2        |      `ConvexFailDomain()`       |
 |                 {ref}`Damped Cosine <test-functions:damped-cosine>`                 |        1        |        `DampedCosine()`         |
 |             {ref}`Damped Oscillator <test-functions:damped-oscillator>`             |        7        |      `DampedOscillator()`       |
@@ -63,8 +64,8 @@ regardless of their typical applications.
 |            {ref}`Oakley & O'Hagan (2002) 1D <test-functions:oakley-1d>`             |        1        |          `Oakley1D()`           |
 |                   {ref}`OTL Circuit <test-functions:otl-circuit>`                   |     6 / 20      |         `OTLCircuit()`          |
 |                  {ref}`Piston Simulation <test-functions:piston>`                   |     7 / 20      |           `Piston()`            |
-|                     {ref}`Robot Arm <test-functions:robot-arm>`                     |        8        |          `RobotArm()`           |
 |             {ref}`Simple Portfolio Model <test-functions:portfolio-3d>`             |        3        |         `Portfolio3D()`         |
+|                     {ref}`Robot Arm <test-functions:robot-arm>`                     |        8        |          `RobotArm()`           |
 |              {ref}`RS - Circular Bar <test-functions:rs-circular-bar>`              |        2        |        `RSCircularBar()`        |
 |                 {ref}`RS - Quadratic <test-functions:rs-quadratic>`                 |        2        |         `RSQuadratic()`         |
 |               {ref}`SaltelliLinear <test-functions:saltelli-linear>`                |        M        |       `SaltelliLinear()`        |
diff --git a/docs/test-functions/coffee-cup.md b/docs/test-functions/coffee-cup.md
new file mode 100644
index 0000000..600171a
--- /dev/null
+++ b/docs/test-functions/coffee-cup.md
@@ -0,0 +1,151 @@
+---
+jupytext:
+  formats: ipynb,md:myst
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.14.1
+kernelspec:
+  display_name: Python 3 (ipykernel)
+  language: python
+  name: python3
+---
+
+(test-functions:coffee-cup)=
+# Cooling Coffee Cup Model
+
+```{code-cell} ipython3
+import numpy as np
+import matplotlib.pyplot as plt
+import uqtestfuns as uqtf
+```
+
+The cooling coffee cup model simulates the temperature evolution of a coffee
+cup as it cools to an ambient temperature by solving an initial value problem
+As a UQ test function, the model is expressed as a two-dimensional,
+vector-valued function.
+
+The model appeared in {cite}`Tennoee2018, Richardson2020`
+as an introductory example for metamodeling.
+
+Some realizations of the temperature evolutions are shown in figure below.
+
+```{code-cell} ipython3
+:tags: [remove-input]
+
+my_fun = uqtf.CoffeeCup()
+my_fun.prob_input.reset_rng(237324)
+num_sample = 20
+xx = my_fun.prob_input.get_sample(num_sample)
+
+yy = my_fun(xx)
+t_e = my_fun.parameters["t_e"]
+n_ts = my_fun.parameters["n_ts"]
+tt = np.linspace(0, t_e, n_ts)
+
+# --- Create temperature evolution
+for i in range(len(xx)):
+    plt.plot(tt, yy[i, :], color="#8da0cb", alpha=0.5)
+plt.grid()
+plt.xlabel("Time [s]", fontsize=14)
+plt.ylabel("Temperature [degC]", fontsize=14)
+
+plt.gcf().tight_layout()
+plt.gcf().set_dpi(150);
+```
+
+## Test function instance
+
+To create a default instance of the test function:
+
+```{code-cell} ipython3
+my_testfun = uqtf.CoffeeCup()
+```
+
+Check if it has been correctly instantiated:
+
+```{code-cell} ipython3
+print(my_testfun)
+```
+
+## Description
+
+The temperature evolution of a coffee cup, $T(t)$, as it cools down
+to an ambient temperature, $T_{\text{amb}}$, is described by
+the following initial value problem (IVP):
+
+$$
+\frac{dT (t)}{dt} = - \kappa (T (t) - T_{\text{amb}}), T \in [ 0, t_e ],
+$$
+
+with an initial condition (IC):
+
+$$
+T (t = 0) = T_0,
+$$
+
+where $\kappa$ is the thermal conductivity of the cup.
+
+The test function is the solution to the IVP:
+
+$$
+\mathcal{M}(\boldsymbol{x}; \boldsymbol{p}) = \left( T(t_i; \boldsymbol{x}, \boldsymbol{p}) \right), \; i = 0, \ldots, n_{ts},
+$$
+
+where:
+
+- $\boldsymbol{x} = \left( \kappa, T_{\text{amb}} \right)$ is a two-dimensional
+  vector of uncertain input variables, defined further below.
+- $\boldsymbol{p} = \{ T_0, t_e, n_{ts} \}$ is a set of fixed parameters
+  of the problem, also defined further below.
+
+## Probabilistic input
+
+The probabilistic input model for the test function is shown below.
+
+```{code-cell} ipython3
+:tags: [hide-input]
+
+print(my_testfun.prob_input)
+```
+
+## Parameters
+
+The default values of the parameters are shown below.
+
+```{code-cell} ipython3
+:tags: [hide-input]
+
+print(my_testfun.parameters)
+```
+
+## Notes on numerical algorithms
+
+The IVP described above is solved numerically using
+[`solve_ivp()`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html)
+fro SciPy with its default method and parameter values.
+
+The default parameter values for these methods can be overridden by providing
+a dictionary with new parameter values.
+
+For example, to change the method used to solve the IVP:
+
+```python
+fun.parameters.add("solve_ivp", {"method": "RK23"})  # 'solve_ivp' as the parameter keyword
+```
+
+```{note}
+In this example, `method` is acceptable keyword-named argument for `solve_ivp()`;
+indeed, the key-value pairs specified for the parameters of `solve_ivp()`
+must be recognized by the method.
+```
+
+## References
+
+```{bibliography}
+:style: unsrtalpha
+:filter: docname in docnames
+```
+
+[^location]: See Eqs. (9-12), Section 3.1 in {cite}`Constantine2015`.
diff --git a/src/uqtestfuns/core/parameters.py b/src/uqtestfuns/core/parameters.py
index 369becf..a79e72c 100644
--- a/src/uqtestfuns/core/parameters.py
+++ b/src/uqtestfuns/core/parameters.py
@@ -312,7 +312,7 @@ def _get_values_as_list(
         if isinstance(value, np.ndarray):
             value = f"{value.shape} array"
         elif isinstance(value, int):
-            value = f"{value:1.0e}"
+            value = f"{value:1.5e}"
         elif isinstance(value, float):
             value = f"{value:1.5e}"
         else:
diff --git a/src/uqtestfuns/core/uqtestfun_abc.py b/src/uqtestfuns/core/uqtestfun_abc.py
index 25ec141..f255f1c 100644
--- a/src/uqtestfuns/core/uqtestfun_abc.py
+++ b/src/uqtestfuns/core/uqtestfun_abc.py
@@ -9,7 +9,7 @@
 
 from copy import deepcopy
 from inspect import signature
-from typing import Callable, cast, List, Optional, Type
+from typing import Callable, cast, List, Optional, Tuple, Type, Union
 
 from .prob_input.marginal import Marginal
 from .prob_input.probabilistic_input import ProbInput
@@ -76,7 +76,6 @@ def __init__(
         self.prob_input = prob_input
         self._parameters = parameters
         self._function_id = function_id
-        self._output_dimension: Optional[int] = None
 
     @property
     def prob_input(self) -> ProbInput:
@@ -119,19 +118,18 @@ def input_dimension(self) -> int:
         return self.prob_input.input_dimension
 
     @property
-    def output_dimension(self) -> int:
-        if self._output_dimension is None:
-            xx = self.prob_input.get_sample()
-            yy = self(xx)
-            if yy.ndim == 1:
-                self._output_dimension = 1
-            elif yy.ndim == 2:
-                self._output_dimension = yy.shape[1]
-            else:
-                self._output_dimension = yy.shape[1:]
-            return self._output_dimension
+    def output_dimension(self) -> Union[int, Tuple[int, ...]]:
+        """The output dimension of the UQ test function."""
+        xx = self.prob_input.get_sample()
+        yy = self(xx)
+        if yy.ndim == 1:
+            output_dim = 1
+        elif yy.ndim == 2:
+            output_dim = yy.shape[1]
+        else:
+            output_dim = yy.shape[1:]
 
-        return self._output_dimension
+        return output_dim
 
     def transform_sample(
         self,
diff --git a/src/uqtestfuns/test_functions/__init__.py b/src/uqtestfuns/test_functions/__init__.py
index 54adf8e..16015b8 100644
--- a/src/uqtestfuns/test_functions/__init__.py
+++ b/src/uqtestfuns/test_functions/__init__.py
@@ -9,6 +9,7 @@
 from .cantilever_beam_2d import CantileverBeam2D
 from .cheng2010 import Cheng2D
 from .circular_pipe_crack import CircularPipeCrack
+from .coffee_cup import CoffeeCup
 from .convex_fail_domain import ConvexFailDomain
 from .damped_cosine import DampedCosine
 from .damped_oscillator import DampedOscillator, DampedOscillatorReliability
@@ -59,6 +60,7 @@
     "CantileverBeam2D",
     "Cheng2D",
     "CircularPipeCrack",
+    "CoffeeCup",
     "ConvexFailDomain",
     "DampedCosine",
     "DampedOscillator",
diff --git a/src/uqtestfuns/test_functions/coffee_cup.py b/src/uqtestfuns/test_functions/coffee_cup.py
new file mode 100644
index 0000000..e753d19
--- /dev/null
+++ b/src/uqtestfuns/test_functions/coffee_cup.py
@@ -0,0 +1,178 @@
+"""
+This module implements the cooling coffee cup model described in Tennøe (2018).
+
+The cooling coffee cup model simulates the temperature evolution of a coffee
+cup as it cools to an ambient temperature by solving an initial value problem
+(using `solve_ivp()` from SciPy).
+As a UQ test function, the model is expressed as a two-dimensional,
+vector-valued function.
+
+The model was featured in [1] and [2] as an introductory example
+for metamodeling.
+
+References
+----------
+
+1. S. Tennøe, G. Halnes, and G. T. Einevoll, “Uncertainpy: A Python Toolbox
+   for Uncertainty Quantification and Sensitivity Analysis
+   in Computational Neuroscience,” Frontier in Neuroinformatics, vol. 12,
+   p. 49, 2018.
+   DOI: 10.3389/fninf.2018.00049
+2. R. A. Richardson, D. W. Wright, W. Edeling, V. Jancauskas, J. Lakhlili,
+   and P. V. Coveney, “EasyVVUQ: A Library for Verification, Validation
+   and Uncertainty Quantification in High Performance Computing,”
+   Journal of Open Research Software, vol. 8, no. 1, p. 11, 2020.
+   DOI: 10.5334/jors.303
+"""
+
+import numpy as np
+from scipy.integrate import solve_ivp
+
+from uqtestfuns.core.custom_typing import (
+    MarginalSpecs,
+    FunParamSpecs,
+    ProbInputSpecs,
+)
+from uqtestfuns.core.uqtestfun_abc import UQTestFunFixDimABC
+
+__all__ = ["CoffeeCup"]
+
+
+MARGINALS_TENNOEE2018: MarginalSpecs = [
+    {
+        "name": "kappa",
+        "distribution": "uniform",
+        "parameters": [0.025, 0.075],
+        "description": "Thermal conductivity of the cup",
+    },
+    {
+        "name": "temp_amb",
+        "distribution": "uniform",
+        "parameters": [15.0, 25.0],
+        "description": "Ambient temperature [degC]",
+    },
+]
+
+AVAILABLE_INPUTS: ProbInputSpecs = {
+    "Tennoee2018": {
+        "function_id": "CoffeeCup",
+        "description": (
+            "Probabilistic input model for the cooling coffee cup model "
+            "from Tennøe et al. (2018)"
+        ),
+        "marginals": MARGINALS_TENNOEE2018,
+        "copulas": None,
+    },
+}
+
+
+AVAILABLE_PARAMETERS: FunParamSpecs = {
+    "Tennoee2018": {
+        "function_id": "CoffeeCup",
+        "description": (
+            "Parameter set for the cooling cup coffee cup model "
+            "Tennøe et al. (2018)"
+        ),
+        "declared_parameters": [
+            {
+                "keyword": "temp_0",
+                "value": 95.0,
+                "type": float,
+                "description": "Initial temperature [degC]",
+            },
+            {
+                "keyword": "t_e",
+                "value": 200.0,
+                "type": float,
+                "description": "End of transient [s]",
+            },
+            {
+                "keyword": "n_ts",
+                "value": 150,
+                "type": int,
+                "description": "Number of time steps",
+            },
+        ],
+    },
+}
+
+
+def fun_ivp(t: float, temp: float, kappa: float, temp_amb: float) -> float:
+    """The right-hand side of the initial value problem.
+
+    Parameters
+    ----------
+    t : float
+        The current time.
+    temp :float
+        The current temperature.
+    kappa : float
+        Thermal conductivity of the cup.
+    temp_amb : float
+        Ambient temperature.
+    """
+    return float(-1 * kappa * (temp - temp_amb))
+
+
+def evaluate(
+    xx: np.ndarray, temp_0: float, t_e: float, n_ts: int, **kwargs
+) -> np.ndarray:
+    """Compute the temperature evolution of the cooling coffee cup.
+
+    Parameters
+    ----------
+    xx : np.ndarray
+        A two-dimensional input values given by an N-by-2 array
+        where N is the number of input values.
+    temp_0 : float
+        The initial temperature of the coffee cup in degC.
+    t_e : float
+        The end of the transient in seconds.
+    n_ts : int
+        The number of time steps in the IVP solution.
+    kwargs : dict
+        The additional parameters as dictionary to be passed to `solve_ivp()`.
+
+    Returns
+    -------
+    np.ndarray
+        The output of the test function evaluated on the input values.
+        The output is a 2-dimensional array of shape (N, n_ts).
+    """
+
+    # Initialize the output
+    yy = np.empty((len(xx), n_ts))
+
+    # Get solve_ivp kwargs
+    solve_ivp_kwargs = kwargs.get("solve_ivp", {})
+
+    for i in range(len(xx)):
+
+        # Get the realization of uncertain inputs
+        kappa = xx[i, 0]
+        temp_amb = xx[i, 1]
+
+        # Solve the IVP
+        sol = solve_ivp(
+            fun_ivp,
+            t_span=(0.0, t_e),
+            y0=[temp_0],
+            t_eval=np.linspace(0.0, t_e, n_ts),
+            args=(kappa, temp_amb),
+            **solve_ivp_kwargs,
+        )
+
+        yy[i, :] = sol.y[0]
+
+    return yy
+
+
+class CoffeeCup(UQTestFunFixDimABC):
+    """Concrete implementation of the cooling coffee cup model."""
+
+    _tags = ["metamodeling"]
+    _description = "Cooling coffee cup model from Tennøe et al. (2018)"
+    _available_inputs = AVAILABLE_INPUTS
+    _available_parameters = AVAILABLE_PARAMETERS
+
+    evaluate = staticmethod(evaluate)  # type: ignore
diff --git a/tests/builtin_test_functions/test_coffee_cup.py b/tests/builtin_test_functions/test_coffee_cup.py
new file mode 100644
index 0000000..fa5aecd
--- /dev/null
+++ b/tests/builtin_test_functions/test_coffee_cup.py
@@ -0,0 +1,70 @@
+"""
+Test module for the Coffee Cup model.
+
+Notes
+-----
+- The tests defined in this module deals with
+  the correctness of the evaluation.
+"""
+
+import numpy as np
+import pytest
+
+from conftest import assert_call
+from uqtestfuns.test_functions import CoffeeCup
+
+
+@pytest.mark.parametrize(
+    "solve_ivp_kwargs", [{"method": "RK23"}, {"dense_output": True}]
+)
+def test_solve_ivp_kwargs(solve_ivp_kwargs):
+    """Test passing additional kwargs for solve_ivp."""
+
+    # Create an instance
+    fun = CoffeeCup()
+
+    # Add the new parameter
+    fun.parameters.add("solve_ivp", solve_ivp_kwargs)
+
+    # Generate test points
+    xx = fun.prob_input.get_sample(10)
+
+    # Assertion
+    assert_call(fun, xx)
+
+
+@pytest.mark.parametrize("n_ts", [10, 20, 50])
+def test_parameter_nts(n_ts):
+    # Create an instance
+    fun = CoffeeCup()
+
+    # Modify the parameters
+    fun.parameters["n_ts"] = n_ts
+
+    # Generate test points
+    xx = fun.prob_input.get_sample(10)
+
+    # Evaluate the function
+    yy = fun(xx)
+
+    # Assertions
+    assert yy.shape[1] == n_ts
+    assert fun.output_dimension == n_ts
+
+
+@pytest.mark.parametrize("temp_0", [50, 60, 70])
+def test_parameter_temp0(temp_0):
+    # Create an instance
+    fun = CoffeeCup()
+
+    # Modify the parameter
+    fun.parameters["temp_0"] = temp_0
+
+    # Generate test points
+    xx = fun.prob_input.get_sample(10)
+
+    # Evaluate the function
+    yy = fun(xx)
+
+    # Assertions
+    assert np.all(yy[:, 0] == temp_0)