code: def _t_run_reactive(self, x):
          
          if x is None:
              x = np.zeros(shape=[self.Np, ], dtype=float)
          self[self.settings['quantity']] = x
          relax = self.settings['relaxation_quantity']
          res = 1e+06
          for itr in range(int(self.settings['max_iter'])):
              if res >= self.settings['r_tolerance']:
                  logger.info('Tolerance not met: ' + str(res))
                  self[self.settings['quantity']] = x
                  self._A = (self._A_t).copy()
                  self._b = (self._b_t).copy()
                  self._apply_sources()
                  x_new = self._solve()
                  # Relaxation
                  x_new = relax*x_new + (1-relax)*self[self.settings['quantity']]
                  self[self.settings['quantity']] = x_new
                  res = np.sum(np.absolute(x**2 - x_new**2))
                  x = x_new
              if (res < self.settings['r_tolerance'] or
                      self.settings['sources'] == []):
                  logger.info('Solution converged: ' + str(res))
                  break
          return x_new
docs: r
          Repeatedly updates transient 'A', 'b', and the solution guess within
          each time step according to the applied source term then calls '_solve'
          to solve the resulting system of linear equations. Stops when the
          residual falls below 'r_tolerance'.
  
          Parameters
          ----------
          x : ND-array
              Initial guess of unknown variable
  
          Returns
          -------
          x_new : ND-array
              Solution array.
  
          Notes
          -----
          Description of 'relaxation_quantity' and 'max_iter' settings can be
          found in the parent class 'ReactiveTransport' documentation.