Repository: mathLab/PINA

Issue:
Wrong link for tutorials .html files
**Describe the bug**
The provided link to the tutorials' .html files does not work properly.

Discussion:

Solution patch:
diff --git a/pina/operator.py b/pina/operator.py
index 0b3a326..cb2a2b6 100644
--- a/pina/operator.py
+++ b/pina/operator.py
@@ -275,44 +275,42 @@ def fast_laplacian(output_, input_, components, d, method="std"):
 def fast_advection(output_, input_, velocity_field, components, d):
     """
     Perform the advection operation on the ``output_`` with respect to the
-    ``input``. This operator support vector-valued functions with multiple input
-    coordinates.
+    ``input``. This operator supports vector-valued functions with multiple
+    input coordinates.
 
     Unlike ``advection``, this function performs no internal checks on input and
     output tensors. The user is required to specify both ``components`` and
     ``d`` as lists of strings. It is designed to enhance computation speed.
 
     :param LabelTensor output_: The output tensor on which the advection is
-        computed.
+        computed. It includes both the velocity and the quantity to be advected.
     :param LabelTensor input_: the input tensor with respect to which advection
         is computed.
-    :param str velocity_field: The name of the output variable used as velocity
-        field. It must be chosen among the output labels.
+    :param list[str] velocity_field: The name of the output variables used as
+        velocity field. It must be chosen among the output labels.
     :param list[str] components: The names of the output variables for which to
         compute the advection. It must be a subset of the output labels.
     :param list[str] d: The names of the input variables with respect to which
         the advection is computed. It must be a subset of the input labels.
     :return: The computed advection tensor.
-    :rtype: torch.Tensor
+    :rtype: LabelTensor
     """
     # Add a dimension to the velocity field for following operations
     velocity = output_.extract(velocity_field).unsqueeze(-1)
 
-    # Remove the velocity field from the components
-    filter_components = [c for c in components if c != velocity_field]
-
     # Compute the gradient
     grads = fast_grad(
-        output_=output_, input_=input_, components=filter_components, d=d
+        output_=output_, input_=input_, components=components, d=d
     )
 
     # Reshape into [..., len(filter_components), len(d)]
-    tmp = grads.reshape(*output_.shape[:-1], len(filter_components), len(d))
+    tmp = grads.reshape(*output_.shape[:-1], len(components), len(d))
 
     # Transpose to [..., len(d), len(filter_components)]
     tmp = tmp.transpose(-1, -2)
 
-    return (tmp * velocity).sum(dim=tmp.tensor.ndim - 2)
+    adv = (tmp * velocity).sum(dim=tmp.tensor.ndim - 2)
+    return LabelTensor(adv, labels=[f"adv_{c}" for c in components])
 
 
 def grad(output_, input_, components=None, d=None):
@@ -425,15 +423,16 @@ def laplacian(output_, input_, components=None, d=None, method="std"):
 def advection(output_, input_, velocity_field, components=None, d=None):
     """
     Perform the advection operation on the ``output_`` with respect to the
-    ``input``. This operator support vector-valued functions with multiple input
-    coordinates.
+    ``input``. This operator supports vector-valued functions with multiple
+    input coordinates.
 
     :param LabelTensor output_: The output tensor on which the advection is
-        computed.
+        computed. It includes both the velocity and the quantity to be advected.
     :param LabelTensor input_: the input tensor with respect to which advection
         is computed.
-    :param str velocity_field: The name of the output variable used as velocity
+    :param velocity_field: The name of the output variables used as velocity
         field. It must be chosen among the output labels.
+    :type velocity_field: str | list[str]
     :param components: The names of the output variables for which to compute
         the advection. It must be a subset of the output labels.
         If ``None``, all output variables are considered. Default is ``None``.
@@ -444,18 +443,29 @@ def advection(output_, input_, velocity_field, components=None, d=None):
     :type d: str | list[str]
     :raises TypeError: If the input tensor is not a LabelTensor.
     :raises TypeError: If the output tensor is not a LabelTensor.
-    :raises RuntimeError: If the velocity field is not in the output labels.
+    :raises RuntimeError: If the velocity field is not a subset of the output
+        labels.
+    :raises RuntimeError: If the dimensionality of the velocity field does not
+        match that of the input tensor.
     :return: The computed advection tensor.
-    :rtype: torch.Tensor
+    :rtype: LabelTensor
     """
     components, d = _check_values(
         output_=output_, input_=input_, components=components, d=d
     )
 
-    # Check if velocity field is present in the output labels
-    if velocity_field not in output_.labels:
+    # Map velocity_field to a list if it is a string
+    if isinstance(velocity_field, str):
+        velocity_field = [velocity_field]
+
+    # Check if all the velocity_field labels are present in the output labels
+    if not all(vi in output_.labels for vi in velocity_field):
+        raise RuntimeError("Velocity labels missing from output tensor.")
+
+    # Check if the velocity has the same dimensionality as the input tensor
+    if len(velocity_field) != len(d):
         raise RuntimeError(
-            f"Velocity {velocity_field} is not present in the output labels."
+            "Velocity dimensionality does not match input dimensionality."
         )
 
     return fast_advection(
diff --git a/tutorials/README.md b/tutorials/README.md
index f2a6322..e432e4a 100644
--- a/tutorials/README.md
+++ b/tutorials/README.md
@@ -9,39 +9,39 @@ The table below provides an overview of each tutorial. All tutorials are also av
 
 | Description   | Tutorial |
 |---------------|-----------|
-Introductory Tutorial: A Beginner’s Guide to PINA|[[.ipynb](tutorial17/tutorial.ipynb),[.py](tutorial17/tutorial.py),[.html](http://mathlab.github.io/PINA/_rst/tutorials/tutorial17/tutorial.html)]|
-How to build a `Problem` in PINA|[[.ipynb](tutorial16/tutorial.ipynb),[.py](tutorial16/tutorial.py),[.html](http://mathlab.github.io/PINA/_rst/tutorials/tutorial16/tutorial.html)]|
-Introduction to Solver classes|[[.ipynb](tutorial18/tutorial.ipynb),[.py](tutorial18/tutorial.py),[.html](http://mathlab.github.io/PINA/_rst/tutorials/tutorial18/tutorial.html)]|
-Introduction to `Trainer` class|[[.ipynb](tutorial11/tutorial.ipynb),[.py](tutorial11/tutorial.py),[.html](http://mathlab.github.io/PINA/_rst/tutorials/tutorial11/tutorial.html)]|
-Data structure for SciML: `Tensor`, `LabelTensor`, `Data` and `Graph` |[[.ipynb](tutorial19/tutorial.ipynb),[.py](tutorial19/tutorial.py),[.html](http://mathlab.github.io/PINA/_rst/tutorials/tutorial19/tutorial.html)]|
-Building geometries with `DomainInterface` class|[[.ipynb](tutorial6/tutorial.ipynb),[.py](tutorial6/tutorial.py),[.html](http://mathlab.github.io/PINA/_rst/tutorials/tutorial6/tutorial.html)]|
-Introduction to PINA `Equation` class|[[.ipynb](tutorial12/tutorial.ipynb),[.py](tutorial12/tutorial.py),[.html](http://mathlab.github.io/PINA/_rst/tutorials/tutorial12/tutorial.html)]|
+Introductory Tutorial: A Beginner’s Guide to PINA|[[.ipynb](tutorial17/tutorial.ipynb),[.py](tutorial17/tutorial.py),[.html](http://mathlab.github.io/PINA/tutorial17/tutorial.html)]|
+How to build a `Problem` in PINA|[[.ipynb](tutorial16/tutorial.ipynb),[.py](tutorial16/tutorial.py),[.html](http://mathlab.github.io/PINA/tutorial16/tutorial.html)]|
+Introduction to Solver classes|[[.ipynb](tutorial18/tutorial.ipynb),[.py](tutorial18/tutorial.py),[.html](http://mathlab.github.io/PINA/tutorial18/tutorial.html)]|
+Introduction to `Trainer` class|[[.ipynb](tutorial11/tutorial.ipynb),[.py](tutorial11/tutorial.py),[.html](http://mathlab.github.io/PINA/tutorial11/tutorial.html)]|
+Data structure for SciML: `Tensor`, `LabelTensor`, `Data` and `Graph` |[[.ipynb](tutorial19/tutorial.ipynb),[.py](tutorial19/tutorial.py),[.html](http://mathlab.github.io/PINA/tutorial19/tutorial.html)]|
+Building geometries with `DomainInterface` class|[[.ipynb](tutorial6/tutorial.ipynb),[.py](tutorial6/tutorial.py),[.html](http://mathlab.github.io/PINA/tutorial6/tutorial.html)]|
+Introduction to PINA `Equation` class|[[.ipynb](tutorial12/tutorial.ipynb),[.py](tutorial12/tutorial.py),[.html](http://mathlab.github.io/PINA/tutorial12/tutorial.html)]|
 
 
 ## Physics Informed Neural Networks
 | Description   | Tutorial  |
 |---------------|-----------|
-Introductory Tutorial:  Physics Informed Neural Networks with PINA |[[.ipynb](tutorial1/tutorial.ipynb),[.py](tutorial1/tutorial.py),[.html](http://mathlab.github.io/PINA/_rst/tutorials/tutorial1/tutorial.html)]|
-Enhancing PINNs with Extra Features to solve the Poisson Problem |[[.ipynb](tutorial2/tutorial.ipynb),[.py](tutorial2/tutorial.py),[.html](http://mathlab.github.io/PINA/_rst/tutorials/tutorial2/tutorial.html)]|
-Applying Hard Constraints in PINNs to solve the Wave Problem |[[.ipynb](tutorial3/tutorial.ipynb),[.py](tutorial3/tutorial.py),[.html](http://mathlab.github.io/PINA/_rst/tutorials/tutorial3/tutorial.html)]|
-Applying Periodic Boundary Conditions in PINNs to solve the Helmotz Problem |[[.ipynb](tutorial9/tutorial.ipynb),[.py](tutorial9/tutorial.py),[.html](http://mathlab.github.io/PINA/_rst/tutorials/tutorial9/tutorial.html)]|
-Inverse Problem Solving with Physics-Informed Neural Network |[[.ipynb](tutorial7/tutorial.ipynb),[.py](tutorial7/tutorial.py),[.html](http://mathlab.github.io/PINA/_rst/tutorials/tutorial7/tutorial.html)]|
-Learning Multiscale PDEs Using Fourier Feature Networks|[[.ipynb](tutorial13/tutorial.ipynb),[.py](tutorial13/tutorial.py),[.html](http://mathlab.github.io/PINA/_rst/tutorials/tutorial13/tutorial.html)]|
-Learning Bifurcating PDE Solutions with Physics-Informed Deep Ensembles|[[.ipynb](tutorial14/tutorial.ipynb),[.py](tutorial14/tutorial.py),[.html](http://mathlab.github.io/PINA/_rst/tutorials/tutorial14/tutorial.html)]|
+Introductory Tutorial:  Physics Informed Neural Networks with PINA |[[.ipynb](tutorial1/tutorial.ipynb),[.py](tutorial1/tutorial.py),[.html](http://mathlab.github.io/PINA/tutorial1/tutorial.html)]|
+Enhancing PINNs with Extra Features to solve the Poisson Problem |[[.ipynb](tutorial2/tutorial.ipynb),[.py](tutorial2/tutorial.py),[.html](http://mathlab.github.io/PINA/tutorial2/tutorial.html)]|
+Applying Hard Constraints in PINNs to solve the Wave Problem |[[.ipynb](tutorial3/tutorial.ipynb),[.py](tutorial3/tutorial.py),[.html](http://mathlab.github.io/PINA/tutorial3/tutorial.html)]|
+Applying Periodic Boundary Conditions in PINNs to solve the Helmholtz Problem |[[.ipynb](tutorial9/tutorial.ipynb),[.py](tutorial9/tutorial.py),[.html](http://mathlab.github.io/PINA/tutorial9/tutorial.html)]|
+Inverse Problem Solving with Physics-Informed Neural Network |[[.ipynb](tutorial7/tutorial.ipynb),[.py](tutorial7/tutorial.py),[.html](http://mathlab.github.io/PINA/tutorial7/tutorial.html)]|
+Learning Multiscale PDEs Using Fourier Feature Networks|[[.ipynb](tutorial13/tutorial.ipynb),[.py](tutorial13/tutorial.py),[.html](http://mathlab.github.io/PINA/tutorial13/tutorial.html)]|
+Learning Bifurcating PDE Solutions with Physics-Informed Deep Ensembles|[[.ipynb](tutorial14/tutorial.ipynb),[.py](tutorial14/tutorial.py),[.html](http://mathlab.github.io/PINA/tutorial14/tutorial.html)]|
 
 
 ## Neural Operator Learning
 | Description   | Tutorial  |
 |---------------|-----------|
-Introductory Tutorial: Neural Operator Learning with PINA |[[.ipynb](tutorial21/tutorial.ipynb),[.py](tutorial21/tutorial.py),[.html](http://mathlab.github.io/PINA/_rst/tutorials/tutorial21/tutorial.html)]|
-Modeling 2D Darcy Flow with the Fourier Neural Operator |[[.ipynb](tutorial5/tutorial.ipynb),[.py](tutorial5/tutorial.py),[.html](http://mathlab.github.io/PINA/_rst/tutorials/tutorial5/tutorial.html)]|
-Solving the Kuramoto–Sivashinsky Equation with Averaging Neural Operator |[[.ipynb](tutorial10/tutorial.ipynb),[.py](tutorial10/tutorial.py),[.html](http://mathlab.github.io/PINA/_rst/tutorials/tutorial10/tutorial.html)]|
+Introductory Tutorial: Neural Operator Learning with PINA |[[.ipynb](tutorial21/tutorial.ipynb),[.py](tutorial21/tutorial.py),[.html](http://mathlab.github.io/PINA/tutorial21/tutorial.html)]|
+Modeling 2D Darcy Flow with the Fourier Neural Operator |[[.ipynb](tutorial5/tutorial.ipynb),[.py](tutorial5/tutorial.py),[.html](http://mathlab.github.io/PINA/tutorial5/tutorial.html)]|
+Solving the Kuramoto–Sivashinsky Equation with Averaging Neural Operator |[[.ipynb](tutorial10/tutorial.ipynb),[.py](tutorial10/tutorial.py),[.html](http://mathlab.github.io/PINA/tutorial10/tutorial.html)]|
 
 ## Supervised Learning
 | Description   | Tutorial  |
 |---------------|-----------|
-Introductory Tutorial: Supervised Learning with PINA |[[.ipynb](tutorial20/tutorial.ipynb),[.py](tutorial20/tutorial.py),[.html](http://mathlab.github.io/PINA/_rst/tutorials/tutorial20/tutorial.html)]|
-Chemical Properties Prediction with Graph Neural Networks |[[.ipynb](tutorial15/tutorial.ipynb),[.py](tutorial15/tutorial.py),[.html](http://mathlab.github.io/PINA/_rst/tutorials/tutorial15/tutorial.html)]|
-Unstructured Convolutional Autoencoders with Continuous Convolution |[[.ipynb](tutorial4/tutorial.ipynb),[.py](tutorial4/tutorial.py),[.html](http://mathlab.github.io/PINA/_rst/tutorials/tutorial4/tutorial.html)]|
-Reduced Order Modeling with POD-RBF and POD-NN Approaches for Fluid Dynamics| [[.ipynb](tutorial8/tutorial.ipynb),[.py](tutorial8/tutorial.py),[.html](http://mathlab.github.io/PINA/_rst/tutorials/tutorial8/tutorial.html)]|
+Introductory Tutorial: Supervised Learning with PINA |[[.ipynb](tutorial20/tutorial.ipynb),[.py](tutorial20/tutorial.py),[.html](http://mathlab.github.io/PINA/tutorial20/tutorial.html)]|
+Chemical Properties Prediction with Graph Neural Networks |[[.ipynb](tutorial15/tutorial.ipynb),[.py](tutorial15/tutorial.py),[.html](http://mathlab.github.io/PINA/tutorial15/tutorial.html)]|
+Unstructured Convolutional Autoencoders with Continuous Convolution |[[.ipynb](tutorial4/tutorial.ipynb),[.py](tutorial4/tutorial.py),[.html](http://mathlab.github.io/PINA/tutorial4/tutorial.html)]|
+Reduced Order Modeling with POD-RBF and POD-NN Approaches for Fluid Dynamics| [[.ipynb](tutorial8/tutorial.ipynb),[.py](tutorial8/tutorial.py),[.html](http://mathlab.github.io/PINA/tutorial8/tutorial.html)]|
 
diff --git a/tutorials/tutorial11/tutorial.ipynb b/tutorials/tutorial11/tutorial.ipynb
index 6ac6a66..07a60c1 100644
--- a/tutorials/tutorial11/tutorial.ipynb
+++ b/tutorials/tutorial11/tutorial.ipynb
@@ -287,7 +287,7 @@
    "metadata": {},
    "source": [
     "<p align=\\\"center\\\">\n",
-    "<img src=\"../static/logging.png\" alt=\\\"Logging API\\\" width=\\\"400\\\"/>\n",
+    "    <img src=\"http://raw.githubusercontent.com/mathLab/PINA/master/tutorials/static/logging.png\" alt=\\\"Logging API\\\" width=\\\"400\\\"/>\n",
     "</p>"
    ]
   },
diff --git a/tutorials/tutorial14/tutorial.ipynb b/tutorials/tutorial14/tutorial.ipynb
index 1312d6a..afafdcb 100644
--- a/tutorials/tutorial14/tutorial.ipynb
+++ b/tutorials/tutorial14/tutorial.ipynb
@@ -58,7 +58,7 @@
     "This approach allows the ensemble to capture different perspectives of the problem, leading to more accurate and reliable predictions.\n",
     "\n",
     "<p align=\"center\">\n",
-    "    <img src=\"../static/deep_ensemble.png\" alt=\"PINA Workflow\" width=\"600\"/>\n",
+    "    <img src=\"http://raw.githubusercontent.com/mathLab/PINA/master/tutorials/static/deep_ensemble.png\" alt=\"Deep ensemble\" width=\"600\"/>\n",
     "</p>\n",
     "\n",
     "The image above illustrates a Deep Ensemble setup, where multiple models attempt to predict the text from an image. While individual models may make errors (e.g., predicting \"PONY\" instead of \"PINA\"), combining their outputs—such as taking the majority vote—often leads to the correct result. This ensemble effect improves reliability by mitigating the impact of individual model biases.\n",
diff --git a/tutorials/tutorial17/tutorial.ipynb b/tutorials/tutorial17/tutorial.ipynb
index e909b99..6d36a25 100644
--- a/tutorials/tutorial17/tutorial.ipynb
+++ b/tutorials/tutorial17/tutorial.ipynb
@@ -11,9 +11,10 @@
     "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/mathLab/PINA/blob/master/tutorials/tutorial17/tutorial.ipynb)\n",
     "\n",
     "<p align=\"left\">\n",
-    "    <img src=\"../static/pina_logo.png\" alt=\"PINA Logo\" width=\"90\"/>\n",
+    "    <img src=\"https://raw.githubusercontent.com/mathLab/PINA/master/readme/pina_logo.png\" alt=\"PINA logo\" width=\"90\"/>\n",
     "</p>\n",
     "\n",
+    "\n",
     "Welcome to **PINA**!\n",
     "\n",
     "PINA [1] is an open-source Python library designed for **Scientific Machine Learning (SciML)** tasks, particularly involving:\n",
@@ -39,7 +40,7 @@
     "## The PINA Workflow \n",
     "\n",
     "<p align=\"center\">\n",
-    "    <img src=\"../static/pina_wokflow.png\" alt=\"PINA Workflow\" width=\"1000\"/>\n",
+    "    <img src=\"http://raw.githubusercontent.com/mathLab/PINA/master/tutorials/static/pina_wokflow.png\" alt=\"PINA Workflow\" width=\"1000\"/>\n",
     "</p>\n",
     "\n",
     "Solving a differential problem in **PINA** involves four main steps:\n",
diff --git a/tutorials/tutorial21/tutorial.ipynb b/tutorials/tutorial21/tutorial.ipynb
index 9bc39e7..69a5f1b 100644
--- a/tutorials/tutorial21/tutorial.ipynb
+++ b/tutorials/tutorial21/tutorial.ipynb
@@ -163,10 +163,9 @@
     "At their core, **Neural Operators** transform an input function $a$ into an output function $u$. The general structure of a Neural Operator consists of three key components:\n",
     "\n",
     "<p align=\"center\">\n",
-    "    <img src=\"../static/neural_operator.png\" alt=\"Neural Operators\" width=\"800\"/>\n",
+    "    <img src=\"http://raw.githubusercontent.com/mathLab/PINA/master/tutorials/static/neural_operator.png\" alt=\"Neural Operators\" width=\"800\"/>\n",
     "</p>\n",
     "\n",
-    "\n",
     "1. **Encoder**: The encoder maps the input into a specific embedding space.\n",
     "\n",
     "2. **Processor**: The processor consists of multiple layers performing **function convolutions**, which is the core computational unit in a Neural Operator. \n",
diff --git a/tutorials/tutorial9/tutorial.ipynb b/tutorials/tutorial9/tutorial.ipynb
index 43fe8db..a0062e2 100644
--- a/tutorials/tutorial9/tutorial.ipynb
+++ b/tutorials/tutorial9/tutorial.ipynb
@@ -4,7 +4,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Tutorial: Applying Periodic Boundary Conditions in PINNs to solve the Helmotz Problem\n",
+    "# Tutorial: Applying Periodic Boundary Conditions in PINNs to solve the Helmholtz Problem\n",
     "\n",
     "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/mathLab/PINA/blob/master/tutorials/tutorial9/tutorial.ipynb)\n",
     "\n",


Test patch:
diff --git a/tests/test_operator.py b/tests/test_operator.py
index 20eb675..aade700 100644
--- a/tests/test_operator.py
+++ b/tests/test_operator.py
@@ -296,22 +296,183 @@ def test_laplacian(f):
         laplacian(output_=output_, input_=input_, components=["a", "b", "c"])
 
 
-def test_advection():
+def test_advection_scalar():
 
-    # Define input and output
+    # Define 3-dimensional input
     input_ = torch.rand((20, 3), requires_grad=True)
     input_ = LabelTensor(input_, ["x", "y", "z"])
-    output_ = LabelTensor(input_**2, ["u", "v", "c"])
 
-    # Define the velocity field
-    velocity = output_.extract(["c"])
+    # Define 3-dimensional velocity field and quantity to be advected
+    velocity = torch.rand((20, 3), requires_grad=True)
+    field = torch.sum(input_**2, dim=-1, keepdim=True)
+
+    # Combine velocity and field into a LabelTensor
+    labels = ["ux", "uy", "uz", "c"]
+    output_ = LabelTensor(torch.cat((velocity, field), dim=1), labels)
+
+    # Compute the pina advection
+    components = ["c"]
+    pina_adv = advection(
+        output_=output_,
+        input_=input_,
+        velocity_field=["ux", "uy", "uz"],
+        components=components,
+        d=["x", "y", "z"],
+    )
+
+    # Compute the true advection
+    grads = 2 * input_
+    true_adv = torch.sum(grads * velocity, dim=grads.ndim - 1, keepdim=True)
+
+    # Check the shape, labels, and value of the advection
+    assert pina_adv.shape == (*output_.shape[:-1], len(components))
+    assert pina_adv.labels == ["adv_c"]
+    assert torch.allclose(pina_adv, true_adv)
+
+    # Should fail if input not a LabelTensor
+    with pytest.raises(TypeError):
+        advection(
+            output_=output_,
+            input_=input_.tensor,
+            velocity_field=["ux", "uy", "uz"],
+        )
+
+    # Should fail if output not a LabelTensor
+    with pytest.raises(TypeError):
+        advection(
+            output_=output_.tensor,
+            input_=input_,
+            velocity_field=["ux", "uy", "uz"],
+        )
+
+    # Should fail for non-existent input labels
+    with pytest.raises(RuntimeError):
+        advection(
+            output_=output_,
+            input_=input_,
+            d=["x", "a"],
+            velocity_field=["ux", "uy", "uz"],
+        )
+
+    # Should fail for non-existent output labels
+    with pytest.raises(RuntimeError):
+        advection(
+            output_=output_,
+            input_=input_,
+            components=["a", "b", "c"],
+            velocity_field=["ux", "uy", "uz"],
+        )
+
+    # Should fail if velocity_field labels are not present in the output labels
+    with pytest.raises(RuntimeError):
+        advection(
+            output_=output_,
+            input_=input_,
+            velocity_field=["ux", "uy", "nonexistent"],
+            components=["c"],
+        )
+
+    # Should fail if velocity_field dimensionality does not match input tensor
+    with pytest.raises(RuntimeError):
+        advection(
+            output_=output_,
+            input_=input_,
+            velocity_field=["ux", "uy"],
+            components=["c"],
+        )
+
+
+def test_advection_vector():
 
-    # Compute the true advection and the pina advection
-    pina_advection = advection(
-        output_=output_, input_=input_, velocity_field="c"
+    # Define 3-dimensional input
+    input_ = torch.rand((20, 3), requires_grad=True)
+    input_ = LabelTensor(input_, ["x", "y", "z"])
+
+    # Define 3-dimensional velocity field
+    velocity = torch.rand((20, 3), requires_grad=True)
+
+    # Define 2-dimensional field to be advected
+    field_1 = torch.sum(input_**2, dim=-1, keepdim=True)
+    field_2 = torch.sum(input_**3, dim=-1, keepdim=True)
+
+    # Combine velocity and field into a LabelTensor
+    labels = ["ux", "uy", "uz", "c1", "c2"]
+    output_ = LabelTensor(
+        torch.cat((velocity, field_1, field_2), dim=1), labels
+    )
+
+    # Compute the pina advection
+    components = ["c1", "c2"]
+    pina_adv = advection(
+        output_=output_,
+        input_=input_,
+        velocity_field=["ux", "uy", "uz"],
+        components=components,
+        d=["x", "y", "z"],
     )
-    true_advection = velocity * 2 * input_.extract(["x", "y"])
 
-    # Check the shape of the advection
-    assert pina_advection.shape == (*output_.shape[:-1], output_.shape[-1] - 1)
-    assert torch.allclose(pina_advection, true_advection)
+    # Compute the true gradients of the fields "c1", "c2"
+    grads1 = 2 * input_
+    grads2 = 3 * input_**2
+
+    # Compute the true advection for each field
+    true_adv1 = torch.sum(grads1 * velocity, dim=grads1.ndim - 1, keepdim=True)
+    true_adv2 = torch.sum(grads2 * velocity, dim=grads2.ndim - 1, keepdim=True)
+    true_adv = torch.cat((true_adv1, true_adv2), dim=-1)
+
+    # Check the shape, labels, and value of the advection
+    assert pina_adv.shape == (*output_.shape[:-1], len(components))
+    assert pina_adv.labels == ["adv_c1", "adv_c2"]
+    assert torch.allclose(pina_adv, true_adv)
+
+    # Should fail if input not a LabelTensor
+    with pytest.raises(TypeError):
+        advection(
+            output_=output_,
+            input_=input_.tensor,
+            velocity_field=["ux", "uy", "uz"],
+        )
+
+    # Should fail if output not a LabelTensor
+    with pytest.raises(TypeError):
+        advection(
+            output_=output_.tensor,
+            input_=input_,
+            velocity_field=["ux", "uy", "uz"],
+        )
+
+    # Should fail for non-existent input labels
+    with pytest.raises(RuntimeError):
+        advection(
+            output_=output_,
+            input_=input_,
+            d=["x", "a"],
+            velocity_field=["ux", "uy", "uz"],
+        )
+
+    # Should fail for non-existent output labels
+    with pytest.raises(RuntimeError):
+        advection(
+            output_=output_,
+            input_=input_,
+            components=["a", "b", "c"],
+            velocity_field=["ux", "uy", "uz"],
+        )
+
+    # Should fail if velocity_field labels are not present in the output labels
+    with pytest.raises(RuntimeError):
+        advection(
+            output_=output_,
+            input_=input_,
+            velocity_field=["ux", "uy", "nonexistent"],
+            components=["c"],
+        )
+
+    # Should fail if velocity_field dimensionality does not match input tensor
+    with pytest.raises(RuntimeError):
+        advection(
+            output_=output_,
+            input_=input_,
+            velocity_field=["ux", "uy"],
+            components=["c"],
+        )