code: def udot(op1, op2):
      
      dot = np.dot(op1.d, op2.d)
      units = op1.units * op2.units
      if dot.shape == ():
          return unyt_quantity(dot, units)
      return unyt_array(dot, units)
docs: Matrix or vector dot product that preserves units
  
      This is a wrapper around np.dot that preserves units.
  
      Examples
      --------
      >>> from unyt import km, s
      >>> a = np.eye(2)*km
      >>> b = (np.ones((2, 2)) * 2)*s
      >>> print(udot(a, b))
      [[2. 2.]
       [2. 2.]] km*s