code: def _apply_sources(self):
          
          if self.settings['t_scheme'] == 'cranknicolson':
              f1 = 0.5
          else:
              f1 = 1
          phase = self.project.phases()[self.settings['phase']]
          relax = self.settings['relaxation_source']
          for item in self.settings['sources']:
              Ps = self.pores(item)
              # Add S1 to diagonal of A
              # TODO: We need this to NOT overwrite the A and b, but create
              # copy, otherwise we have to regenerate A and b on each loop
              datadiag = self._A.diagonal().copy()
              # Source term relaxation
              S1_old = phase[item+'.'+'S1'][Ps].copy()
              S2_old = phase[item+'.'+'S2'][Ps].copy()
              self._update_physics()
              S1 = phase[item+'.'+'S1'][Ps]
              S2 = phase[item+'.'+'S2'][Ps]
              S1 = relax*S1 + (1-relax)*S1_old
              S2 = relax*S2 + (1-relax)*S2_old
              phase[item+'.'+'S1'][Ps] = S1
              phase[item+'.'+'S2'][Ps] = S2
              datadiag[Ps] = datadiag[Ps] - f1*S1
              # Add S1 to A
              self._A.setdiag(datadiag)
              # Add S2 to b
              self._b[Ps] = self._b[Ps] + f1*S2
docs: r
          Update 'A' and 'b' applying source terms to specified pores
  
          Notes
          -----
          Applying source terms to 'A' and 'b' is performed after (optionally)
          under-relaxing the source term to improve numerical stability. Physics
          are also updated before applying source terms to ensure that source
          terms values are associated with the current value of 'quantity'.
          In the case of a transient simulation, the updates in 'A' and 'b'
          also depend on the time scheme.