diff --git a/josie/FourEq/eos.py b/josie/FourEq/eos.py
index 400148ef..0164c681 100644
--- a/josie/FourEq/eos.py
+++ b/josie/FourEq/eos.py
@@ -3,7 +3,7 @@
 # SPDX-License-Identifier: BSD-3-Clause
 
 from josie.twofluid.state import PhasePair
-import abc
+from abc import ABC, abstractmethod
 import numpy as np
 
 
@@ -12,19 +12,19 @@
 ArrayAndScalar = Union[np.ndarray, float]
 
 
-class BarotropicEOS(metaclass=abc.ABCMeta):
+class BarotropicEOS(ABC):
     """An Abstract Base Class representing a barotropic EOS for an
     Euler System"""
 
-    @abc.abstractmethod
+    @abstractmethod
     def p(self, rho: ArrayAndScalar) -> ArrayAndScalar:
         raise NotImplementedError
 
-    @abc.abstractmethod
+    @abstractmethod
     def rho(self, p: ArrayAndScalar) -> ArrayAndScalar:
         raise NotImplementedError
 
-    @abc.abstractmethod
+    @abstractmethod
     def sound_velocity(self, rho: ArrayAndScalar) -> ArrayAndScalar:
         raise NotImplementedError
 
diff --git a/josie/bn/closure.py b/josie/bn/closure.py
index 25449822..b996f3a1 100644
--- a/josie/bn/closure.py
+++ b/josie/bn/closure.py
@@ -2,19 +2,19 @@
 #
 # SPDX-License-Identifier: BSD-3-Clause
 
-import abc
+from abc import ABC, abstractmethod
 import numpy as np
 
 from .state import Q
 
 
-class Closure(metaclass=abc.ABCMeta):
+class Closure(ABC):
     r"""A class representing the closure relation for :math:`p_I` and
     :math:`\vb{u}_I`. Use them as mixin with the Equation of State in order
     to provide full closure for the system
     """
 
-    @abc.abstractmethod
+    @abstractmethod
     def pI(self, state_array: Q) -> np.ndarray:
         r"""
         Parameters
@@ -32,7 +32,7 @@ def pI(self, state_array: Q) -> np.ndarray:
 
         raise NotImplementedError
 
-    @abc.abstractmethod
+    @abstractmethod
     def uI(self, state_array: Q) -> np.ndarray:
         r"""
         Parameters
diff --git a/josie/boundary/boundary.py b/josie/boundary/boundary.py
index ab4b6cb7..2ac40f08 100644
--- a/josie/boundary/boundary.py
+++ b/josie/boundary/boundary.py
@@ -7,7 +7,7 @@
 """
 from __future__ import annotations
 
-import abc
+from abc import ABC, abstractmethod
 import matplotlib.pyplot as plt
 import numpy as np
 
@@ -89,7 +89,7 @@ def apply_bc(self, cells: MeshCellSet, t: float):
         self.curve.bc(cells, self, t)
 
 
-class BoundaryCurve(metaclass=abc.ABCMeta):
+class BoundaryCurve(ABC):
     r"""A class representing a :class:`BoundaryCurve`. A
     :class:`BoundaryCurve` is parametrized with a single parameter. It
     implements a :meth:`__call__` method that returns the :math:`(x,y)` values
@@ -104,7 +104,7 @@ def bc(self):
     def bc(self, bc: BoundaryCondition):
         self._bc = bc
 
-    @abc.abstractmethod
+    @abstractmethod
     def __call__(self, xi):
         r"""The effective parametrization of the BoundaryCurve. Assume ``xi``
         (:math:`\xi`) to range into :math:`[0, 1]`
diff --git a/josie/euler/eos.py b/josie/euler/eos.py
index 83c93853..880e0236 100644
--- a/josie/euler/eos.py
+++ b/josie/euler/eos.py
@@ -5,7 +5,7 @@
 """ This module contains the different Equation of State (EOS) implementations
 """
 
-import abc
+from abc import ABC, abstractmethod
 import numpy as np
 
 from typing import Union
@@ -13,22 +13,22 @@
 ArrayAndScalar = Union[np.ndarray, float]
 
 
-class EOS(metaclass=abc.ABCMeta):
+class EOS(ABC):
     """An Abstract Base Class representing an EOS for an Euler System"""
 
-    @abc.abstractmethod
+    @abstractmethod
     def rhoe(self, rho: ArrayAndScalar, p: ArrayAndScalar) -> ArrayAndScalar:
         raise NotImplementedError
 
-    @abc.abstractmethod
+    @abstractmethod
     def p(self, rho: ArrayAndScalar, e: ArrayAndScalar) -> ArrayAndScalar:
         raise NotImplementedError
 
-    @abc.abstractmethod
+    @abstractmethod
     def rho(self, p: ArrayAndScalar, e: ArrayAndScalar) -> ArrayAndScalar:
         raise NotImplementedError
 
-    @abc.abstractmethod
+    @abstractmethod
     def sound_velocity(
         self, rho: ArrayAndScalar, p: ArrayAndScalar
     ) -> ArrayAndScalar:
diff --git a/josie/fields.py b/josie/fields.py
index e081fe44..d8ceeb9d 100644
--- a/josie/fields.py
+++ b/josie/fields.py
@@ -33,7 +33,7 @@ def __repr__(self):
 
 
 class FieldsMeta(type):
-    """This metaclass reproduces in a simpler form the behaviour of
+    """This metaclass reproduces in a simpler form the behavior of
     :class:`Enum`.  It tracks all the defined attributes of a class, it
     precomputes the number of fields and replaces the fields that have no int
     value with the corresponding int
diff --git a/josie/mesh/cell.py b/josie/mesh/cell.py
index 36349cb2..9b40ac0c 100644
--- a/josie/mesh/cell.py
+++ b/josie/mesh/cell.py
@@ -4,7 +4,7 @@
 
 from __future__ import annotations
 
-import abc
+from abc import ABC, abstractmethod
 
 import numpy as np
 
@@ -22,7 +22,7 @@
     from josie.mesh import Mesh
 
 
-class Cell(metaclass=abc.ABCMeta):
+class Cell(ABC):
     """This is a class interface representing a generic cell of a
     :class:`Mesh`.
 
@@ -44,35 +44,35 @@ class Cell(metaclass=abc.ABCMeta):
     """
 
     @property
-    @abc.abstractmethod
+    @abstractmethod
     def num_points(self) -> int:
         raise NotImplementedError
 
     @property
-    @abc.abstractmethod
+    @abstractmethod
     def num_dofs(self) -> int:
         raise NotImplementedError
 
     @property
-    @abc.abstractmethod
+    @abstractmethod
     def _meshio_cell_type(self) -> str:
         raise NotImplementedError
 
     @classmethod
-    @abc.abstractmethod
+    @abstractmethod
     def _generate_ghosts(cls, mesh: Mesh):
         """Generate ghost cells centroids"""
         raise NotImplementedError
 
     @classmethod
-    @abc.abstractmethod
+    @abstractmethod
     def centroid(
         cls, nw: PointType, sw: PointType, se: PointType, ne: PointType
     ) -> PointType:
         """Compute the centroid of the cell"""
 
     @classmethod
-    @abc.abstractmethod
+    @abstractmethod
     def volume(
         cls, nw: PointType, sw: PointType, se: PointType, ne: PointType
     ) -> np.ndarray:
@@ -80,19 +80,19 @@ def volume(
         raise NotImplementedError
 
     @classmethod
-    @abc.abstractmethod
+    @abstractmethod
     def face_surface(cls, p0: PointType, p1: PointType) -> np.ndarray:
         """Compute the surface of a face from its points."""
         raise NotImplementedError
 
     @classmethod
-    @abc.abstractmethod
+    @abstractmethod
     def face_normal(cls, p0: PointType, p1: PointType) -> np.ndarray:
         """Compute the normal vector to a face."""
         raise NotImplementedError
 
     @classmethod
-    @abc.abstractmethod
+    @abstractmethod
     def create_connectivity(cls, mesh: Mesh):
         r"""This method creates the connectivity from the given points of
         a mesh. It modifies attributes of the :class:`Mesh` instance.
@@ -123,7 +123,7 @@ def create_connectivity(cls, mesh: Mesh):
         mesh.cells.compute_min_length()
 
     @classmethod
-    @abc.abstractmethod
+    @abstractmethod
     def export_connectivity(cls, mesh: "Mesh") -> MeshIO:
         """This method exports the connectivity of the mesh in the format
         accepted by the :class:`~meshio.Mesh`.
diff --git a/josie/plot/backend.py b/josie/plot/backend.py
index 447c6aec..3edced5c 100644
--- a/josie/plot/backend.py
+++ b/josie/plot/backend.py
@@ -5,7 +5,7 @@
 """ Backends used to display mesh and mesh results """
 from __future__ import annotations
 
-import abc
+from abc import ABC, abstractmethod
 
 from typing import List, Union, TYPE_CHECKING
 
@@ -14,10 +14,10 @@
     from josie.solver import Solver
 
 
-class PlotBackend(metaclass=abc.ABCMeta):
+class PlotBackend(ABC):
     """An abstract interface representing a plot backend"""
 
-    @abc.abstractmethod
+    @abstractmethod
     def plot(self, mesh: Mesh):
         """Allocate a drawing instance in order to draw a single image plot.
         The drawing state is stored in :attr:`plot_state`.
@@ -30,7 +30,7 @@ def plot(self, mesh: Mesh):
 
         raise NotImplementedError
 
-    @abc.abstractmethod
+    @abstractmethod
     def append(self, solver: Solver, t):
         """Appends a new simulation time state
 
@@ -46,7 +46,7 @@ def append(self, solver: Solver, t):
 
         raise NotImplementedError
 
-    @abc.abstractmethod
+    @abstractmethod
     def update(self, solver: Solver):
         """Updates the :attr:`plot_state` with the state (i.e. the field data
         stored in the mesh, e.g. cell data) in the Solver.
@@ -63,7 +63,7 @@ def update(self, solver: Solver):
 
         raise NotImplementedError
 
-    @abc.abstractmethod
+    @abstractmethod
     def show(self, fields: Union[List, str]):
         """Show on screen a list of fields.
 
@@ -76,12 +76,12 @@ def show(self, fields: Union[List, str]):
 
         raise NotImplementedError
 
-    @abc.abstractmethod
+    @abstractmethod
     def show_grid(self):
         """Show the grid on screen"""
         raise NotImplementedError
 
-    @abc.abstractmethod
+    @abstractmethod
     def show_all(self):
         """Show on screen all the fields"""
 
diff --git a/josie/problem.py b/josie/problem.py
index 1834f2a8..8464b384 100644
--- a/josie/problem.py
+++ b/josie/problem.py
@@ -4,7 +4,7 @@
 
 from __future__ import annotations
 
-import abc
+from abc import ABC, abstractmethod
 import numpy as np
 
 from typing import Union
@@ -12,7 +12,7 @@
 from josie.mesh.cellset import CellSet, MeshCellSet
 
 
-class Problem(metaclass=abc.ABCMeta):
+class Problem(ABC):
     r"""A class representing a physical problem to be solved (governed by
     a PDE).
 
@@ -32,6 +32,7 @@ class Problem(metaclass=abc.ABCMeta):
     def __init__(self, **kwargs):
         super().__init__()
 
+    @abstractmethod
     def F(self, values) -> np.ndarray:
         r"""The convective flux operator :math:`\pdeConvective`
 
@@ -47,8 +48,9 @@ def F(self, values) -> np.ndarray:
             \times N_\text{dim}` containing the computed convective flux for
             the given problem
         """
-        pass
+        ...
 
+    @abstractmethod
     def B(self, values) -> np.ndarray:
         r"""This returns the tensor that pre-multiplies the non-conservative
         term of the problem.
@@ -76,6 +78,7 @@ def B(self, values) -> np.ndarray:
         """
         pass
 
+    @abstractmethod
     def K(self, cells: Union[CellSet, MeshCellSet]) -> np.ndarray:
         r"""This returns the tensor that pre-multiplies the gradient in the
         diffusive term of the problem.
@@ -103,6 +106,7 @@ def K(self, cells: Union[CellSet, MeshCellSet]) -> np.ndarray:
         """
         pass
 
+    @abstractmethod
     def s(self, values, t: float) -> np.ndarray:
         r"""This returns the values of the source terms