diff --git a/.github/workflows/Fortran.yml b/.github/workflows/Fortran.yml
index 2df5ecb1a..837b90a87 100644
--- a/.github/workflows/Fortran.yml
+++ b/.github/workflows/Fortran.yml
@@ -2,7 +2,7 @@ name: Grid and Mesh Solver
 on: [push]
 
 env:
-  PETSC_VERSION: '3.22.0'
+  PETSC_VERSION: '3.22.2'
   HOMEBREW_NO_ANALYTICS: 'ON'     # Make Homebrew installation a little quicker
   HOMEBREW_NO_AUTO_UPDATE: 'ON'
   HOMEBREW_NO_BOTTLE_SOURCE_FALLBACK: 'ON'
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index f39c2bdba..1a78ec2d1 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -34,22 +34,14 @@ variables:
   # ===============================================================================================
   # Names of module files to load
   # ===============================================================================================
-  # ++++++++++++ Compiler +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-  COMPILER_GNU: "Compiler/GNU/10"
-  COMPILER_INTELLLVM: "Compiler/oneAPI/2022.0.1 Libraries/IMKL/2022.0.1"
+  # ++++++++++++ Compiler/MPI/PETSc +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   COMPILER_INTEL: "Compiler/Intel/2022.0.1 Libraries/IMKL/2022.0.1"
-  # ++++++++++++ MPI ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-  MPI_GNU: "MPI/GNU/10/OpenMPI/4.1.2"
-  MPI_INTELLLVM: "MPI/oneAPI/2022.0.1/IntelMPI/2021.5.0"
   MPI_INTEL: "MPI/Intel/2022.0.1/IntelMPI/2021.5.0"
-  # ++++++++++++ PETSc ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-  PETSC_GNU: "Libraries/PETSc/3.16.4/GNU-10-OpenMPI-4.1.2"
-  PETSC_INTELLLVM: "Libraries/PETSc/3.16.3/oneAPI-2022.0.1-IntelMPI-2021.5.0"
   PETSC_INTEL: "Libraries/PETSc/3.16.5/Intel-2022.0.1-IntelMPI-2021.5.0"
   # ++++++++++++ MSC Marc +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   MSC: "FEM/MSC/2024.1"
-  IntelMarc: "Compiler/Intel/19.1.2 Libraries/IMKL/2020"
-  HDF5Marc: "HDF5/1.12.2/Intel-19.1.2"
+  IntelMarc: "Compiler/Intel/2022.0.1 Libraries/IMKL/2022.0.1"
+  HDF5Marc: "HDF5/1.12.1/Intel-2022.0.1"
 
 
 ###################################################################################################
@@ -330,7 +322,7 @@ grid_performance:
 
 
 ###################################################################################################
-update_website:
+update_website_master:
   stage: finalize
   trigger:
     project: damask/website
@@ -338,6 +330,14 @@ update_website:
   only:
     - development
 
+update_website_3.0:
+  stage: finalize
+  trigger:
+    project: damask/website
+    branch: "3.0"
+  only:
+    - "3.0"
+
 update_revision:
   stage: finalize
   before_script:
diff --git a/python/damask/_colormap.py b/python/damask/_colormap.py
index c6c245c72..0e09c4bf3 100644
--- a/python/damask/_colormap.py
+++ b/python/damask/_colormap.py
@@ -170,11 +170,11 @@ def from_range(low: FloatSequence,
         model : {'rgb', 'hsv', 'hsl', 'xyz', 'lab', 'msh'}
             Color model used for input color definitions. Defaults to 'rgb'.
             The available color models are:
-            - 'rgb': Red Green Blue.
-            - 'hsv': Hue Saturation Value.
-            - 'hsl': Hue Saturation Luminance.
-            - 'xyz': CIE Xyz.
-            - 'lab': CIE Lab.
+            - 'rgb': Red Green Blue (RGB).
+            - 'hsv': Hue Saturation Value (HSV).
+            - 'hsl': Hue Saturation Luminance (HSL).
+            - 'xyz': CIEXYZ.
+            - 'lab': CIELAB.
             - 'msh': Msh (for perceptually uniform interpolation).
 
         Returns
@@ -182,6 +182,12 @@ def from_range(low: FloatSequence,
         new : damask.Colormap
             Colormap spanning given bounds.
 
+        Notes
+        -----
+        RGB values are in the range [[0,1],  [0,1],[0,1]].
+        HSV values are in the range [[0,360],[0,1],[0,1]].
+        HSL values are in the range [[0,360],[0,1],[0,1]].
+
         Examples
         --------
         >>> import damask
@@ -538,12 +544,12 @@ def _hsv2rgb(hsv: np.ndarray) -> np.ndarray:
         Parameters
         ----------
         hsv : numpy.ndarray, shape (3)
-            HSV values.
+            HSV values in the range [[0,360],[0,1],[0,1]].
 
         Returns
         -------
         rgb : numpy.ndarray, shape (3)
-            RGB values.
+            RGB values in the range [[0,1],[0,1],[0,1]].
 
         """
         return np.array(colorsys.hsv_to_rgb(hsv[0]/360.,hsv[1],hsv[2]))
@@ -556,12 +562,12 @@ def _rgb2hsv(rgb: np.ndarray) -> np.ndarray:
         Parameters
         ----------
         rgb : numpy.ndarray, shape (3)
-            RGB values.
+            RGB values in the range [[0,1],[0,1],[0,1]].
 
         Returns
         -------
         hsv : numpy.ndarray, shape (3)
-            HSV values.
+            HSV values in the range [[0,360],[0,1],[0,1]].
 
         """
         h,s,v = colorsys.rgb_to_hsv(rgb[0],rgb[1],rgb[2])
@@ -576,12 +582,12 @@ def _hsl2rgb(hsl: np.ndarray) -> np.ndarray:
         Parameters
         ----------
         hsl : numpy.ndarray, shape (3)
-            HSL values.
+            HSL values in the range [[0,360],[0,1],[0,1]].
 
         Returns
         -------
         rgb : numpy.ndarray, shape (3)
-            RGB values.
+            RGB values in the range [[0,1],[0,1],[0,1]].
 
         """
         return np.array(colorsys.hls_to_rgb(hsl[0]/360.,hsl[2],hsl[1]))
@@ -594,12 +600,12 @@ def _rgb2hsl(rgb: np.ndarray) -> np.ndarray:
         Parameters
         ----------
         rgb : numpy.ndarray, shape (3)
-            RGB values.
+            RGB values in the range [[0,1],[0,1],[0,1]].
 
         Returns
         -------
         hsl : numpy.ndarray, shape (3)
-            HSL values.
+            HSL values in the range [[0,360],[0,1],[0,1]].
 
         """
         h,l,s = colorsys.rgb_to_hls(rgb[0],rgb[1],rgb[2])
@@ -609,17 +615,17 @@ def _rgb2hsl(rgb: np.ndarray) -> np.ndarray:
     @staticmethod
     def _xyz2rgb(xyz: np.ndarray) -> np.ndarray:
         """
-        CIE Xyz to Red Green Blue.
+        CIEXYZ to Red Green Blue.
 
         Parameters
         ----------
         xyz : numpy.ndarray, shape (3)
-            CIE Xyz values.
+            CIEXYZ values.
 
         Returns
         -------
         rgb : numpy.ndarray, shape (3)
-            RGB values.
+            RGB values in the range [[0,1],[0,1],[0,1]].
 
         References
         ----------
@@ -639,17 +645,17 @@ def _xyz2rgb(xyz: np.ndarray) -> np.ndarray:
     @staticmethod
     def _rgb2xyz(rgb: np.ndarray) -> np.ndarray:
         """
-        Red Green Blue to CIE Xyz.
+        Red Green Blue to CIEXYZ.
 
         Parameters
         ----------
         rgb : numpy.ndarray, shape (3)
-            RGB values.
+            RGB values in the range [[0,1],[0,1],[0,1]].
 
         Returns
         -------
         xyz : numpy.ndarray, shape (3)
-            CIE Xyz values.
+            CIEXYZ values.
 
         References
         ----------
@@ -668,22 +674,23 @@ def _rgb2xyz(rgb: np.ndarray) -> np.ndarray:
     def _lab2xyz(lab: np.ndarray,
                  ref_white: np.ndarray = _REF_WHITE) -> np.ndarray:
         """
-        CIE Lab to CIE Xyz.
+        CIELAB to CIEXYZ.
 
         Parameters
         ----------
         lab : numpy.ndarray, shape (3)
-            CIE lab values.
+            CIELAB values.
         ref_white : numpy.ndarray, shape (3)
             Reference white, default value is the standard 2° observer for D65.
 
         Returns
         -------
         xyz : numpy.ndarray, shape (3)
-            CIE Xyz values.
+            CIEXYZ values.
 
         References
         ----------
+        https://en.wikipedia.org/wiki/CIELAB_color_space
         http://www.brucelindbloom.com/index.html?Eqn_Lab_to_XYZ.html
 
         """
@@ -700,26 +707,29 @@ def _lab2xyz(lab: np.ndarray,
     def _xyz2lab(xyz: np.ndarray,
                  ref_white: np.ndarray = _REF_WHITE) -> np.ndarray:
         """
-        CIE Xyz to CIE Lab.
+        CIEXYZ to CIELAB.
 
         Parameters
         ----------
         xyz : numpy.ndarray, shape (3)
-            CIE Xyz values.
+            CIEXYZ values.
         ref_white : numpy.ndarray, shape (3)
             Reference white, default value is the standard 2° observer for D65.
 
         Returns
         -------
         lab : numpy.ndarray, shape (3)
-            CIE lab values.
+            CIELAB values.
 
         References
         ----------
+        https://en.wikipedia.org/wiki/CIELAB_color_space
         http://www.brucelindbloom.com/index.html?Eqn_Lab_to_XYZ.html
 
         """
-        f = np.where(xyz/ref_white > _EPS,(xyz/ref_white)**(1./3.),(_KAPPA*xyz/ref_white+16.)/116.)
+        f = np.where(xyz/ref_white > _EPS,
+                     (xyz/ref_white)**(1./3.),
+                     (_KAPPA*xyz/ref_white+16.)/116.)
 
         return np.array([
                          116.0 *  f[1] - 16.0,
@@ -731,12 +741,12 @@ def _xyz2lab(xyz: np.ndarray,
     @staticmethod
     def _lab2msh(lab: np.ndarray) -> np.ndarray:
         """
-        CIE Lab to Msh.
+        CIELAB to Msh.
 
         Parameters
         ----------
         lab : numpy.ndarray, shape (3)
-            CIE lab values.
+            CIELAB values.
 
         Returns
         -------
@@ -759,7 +769,7 @@ def _lab2msh(lab: np.ndarray) -> np.ndarray:
     @staticmethod
     def _msh2lab(msh: np.ndarray) -> np.ndarray:
         """
-        Msh to CIE Lab.
+        Msh to CIELAB.
 
         Parameters
         ----------
@@ -769,7 +779,7 @@ def _msh2lab(msh: np.ndarray) -> np.ndarray:
         Returns
         -------
         lab : numpy.ndarray, shape (3)
-            CIE lab values.
+            CIELAB values.
 
         References
         ----------
diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py
index 8a53fd8ef..bc1371eb9 100644
--- a/python/damask/_rotation.py
+++ b/python/damask/_rotation.py
@@ -423,7 +423,7 @@ def __matmul__(self,
         >>> import damask
         >>> r = damask.Rotation.from_random(shape=(12))
         >>> o = np.ones((5,3))
-        >>> (r@o).shape                                    # (12) @ (5, 3)
+        >>> (r@o).shape                                                                             # (12) @ (5, 3)
         (12, 5, 3)
 
         Application of a (random) rotation to all twelve second-rank tensors.
@@ -432,7 +432,7 @@ def __matmul__(self,
         >>> import damask
         >>> r = damask.Rotation.from_random()
         >>> o = np.ones((12,3,3))
-        >>> (r@o).shape                                    # (1) @ (12, 3,3)
+        >>> (r@o).shape                                                                             # (1) @ (12, 3,3)
         (12, 3, 3)
 
         Application of twelve (random) rotations to the corresponding twelve second-rank tensors.
@@ -441,7 +441,7 @@ def __matmul__(self,
         >>> import damask
         >>> r = damask.Rotation.from_random(shape=(12))
         >>> o = np.ones((12,3,3))
-        >>> (r@o).shape                                    # (12) @ (3,3)
+        >>> (r@o).shape                                                                             # (12) @ (3,3)
         (12, 3, 3)
 
         Application of each of three (random) rotations to all three vectors.
@@ -450,7 +450,7 @@ def __matmul__(self,
         >>> import damask
         >>> r = damask.Rotation.from_random(shape=(3))
         >>> o = np.ones((3,3))
-        >>> (r[...,np.newaxis]@o[np.newaxis,...]).shape    # (3,1) @ (1,3, 3)
+        >>> (r[...,np.newaxis]@o[np.newaxis,...]).shape                                             # (3,1) @ (1,3, 3)
         (3, 3, 3)
 
         Application of twelve (random) rotations to all twelve second-rank tensors.
@@ -611,14 +611,14 @@ def average(self: MyType,
         https://doi.org/10.2514/1.28949
 
         """
-        def _M(quat):
+        def M(quat):
             """Intermediate representation supporting quaternion averaging."""
             return np.einsum('...i,...j',quat,quat)
 
         weights_ = np.ones(self.shape,dtype=float) if weights is None else np.array(weights,float)
 
-        eig, vec = np.linalg.eig(np.sum(_M(self.quaternion) * weights_[...,np.newaxis,np.newaxis],axis=-3)
-                                /np.sum(                      weights_[...,np.newaxis,np.newaxis],axis=-3))
+        eig, vec = np.linalg.eig(np.sum(M(self.quaternion) * weights_[...,np.newaxis,np.newaxis],axis=-3)
+                                /np.sum(                     weights_[...,np.newaxis,np.newaxis],axis=-3))
 
         return self.copy(Rotation.from_quaternion(np.real(
                                                   np.squeeze(
diff --git a/python/damask/tensor.py b/python/damask/tensor.py
index 63eb15b99..6bbb2edf1 100644
--- a/python/damask/tensor.py
+++ b/python/damask/tensor.py
@@ -59,7 +59,7 @@ def eigenvectors(T_sym: _np.ndarray,
 
     Returns
     -------
-    x : numpy.ndarray, shape (...,3,3)
+    v : numpy.ndarray, shape (...,3,3)
         Eigenvectors of T_sym sorted in ascending order of their
         associated eigenvalues.
 
diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90
index de71e8efa..30cd42011 100644
--- a/src/phase_mechanical.f90
+++ b/src/phase_mechanical.f90
@@ -1279,7 +1279,7 @@ end subroutine mechanical_restartRead
 
 
 !--------------------------------------------------------------------------------------------------
-!< @brief Get first Piola-Kirchhoff stress (for use by non-mech physics).
+!< @brief Get second Piola-Kirchhoff stress (for use by non-mech physics).
 !--------------------------------------------------------------------------------------------------
 module function mechanical_S(ph,en) result(S)
 
@@ -1335,7 +1335,7 @@ end function mechanical_F_i
 
 
 !--------------------------------------------------------------------------------------------------
-!< @brief Get second Piola-Kirchhoff stress (for use by homogenization).
+!< @brief Get first Piola-Kirchhoff stress (for use by homogenization).
 !--------------------------------------------------------------------------------------------------
 module function phase_P(co,ce) result(P)