smile_Kd_eps.hpp

Go to the documentation of this file.

/*
 * smile_Kd.hpp
 */
 
#ifndef _SMILE_KD_HPP_
#define _SMILE_KD_HPP_
 
#include "lib_disc/spatial_disc/user_data/linker/linker.h"
 
#include "smile_chem_eq.hpp"
#include "smile_cloud.hpp"
 
namespace ug {
namespace smile {
 
template <int dim>
class SmartKd
  : public StdDataLinker< SmartKd<dim>, number, dim>
{
  typedef StdDataLinker< SmartKd<dim>, number, dim> base_type;
  
  static inline void compute_I_pH_m3
  (
    number s, number H, number Ca, number Al, number DIC, number SO4, // [i]
    number& I, number& pH // [o]
  )
  {
    Compute_I_pH (NaCl_sat*s/1000, H/1000, Ca/1000, Al/1000, DIC/1000, SO4/1000, I, pH);
  }
 
  public:
  SmartKd(Cloud *cloud,
      SmartPtr<CplUserData<number, dim> > spH,
      SmartPtr<CplUserData<number, dim> > spS,
      SmartPtr<CplUserData<number, dim> > spCa,
      SmartPtr<CplUserData<number, dim> > spDIC,
      SmartPtr<CplUserData<number, dim> > spSO4,
      SmartPtr<CplUserData<number, dim> > spAl,
      number entry)
 
    : m_cloud(cloud), m_spH(spH), m_spS(spS), m_spCa(spCa), m_spDIC(spDIC),
      m_spSO4(spSO4), m_spAl(spAl), m_entry(entry)
  {
    this->set_num_input(6);
    m_spDH   =   m_spH.template cast_dynamic<DependentUserData<number, dim> >();
    m_spDS   =   m_spS.template cast_dynamic<DependentUserData<number, dim> >();
    m_spDCa  =  m_spCa.template cast_dynamic<DependentUserData<number, dim> >();
    m_spDDIC = m_spDIC.template cast_dynamic<DependentUserData<number, dim> >();
    m_spDSO4 = m_spSO4.template cast_dynamic<DependentUserData<number, dim> >();
    m_spDAl  =  m_spAl.template cast_dynamic<DependentUserData<number, dim> >();
    base_type::set_input(0, spH,   spH  );
    base_type::set_input(1, spS,   spS  );
    base_type::set_input(2, spCa,  spCa );
    base_type::set_input(3, spDIC, spDIC);
    base_type::set_input(4, spSO4, spSO4);
    base_type::set_input(5, spAl,  spAl );
  }
 
  inline void evaluate (number& Kd,
              const MathVector<dim>& globIP,
              number time, int si) const
  {
    number H;
    (*m_spH)  (H,   globIP, time, si);
    number s;
    (*m_spS)  (s,   globIP, time, si);
    number Ca;
    (*m_spCa) (Ca,  globIP, time, si);
    number DIC;
    (*m_spDIC)(DIC, globIP, time, si);
    number SO4;
    (*m_spSO4)(SO4, globIP, time, si);
    number Al;
    (*m_spAl)(Al,  globIP, time, si);
 
    number I, pH;
    compute_I_pH_m3(s, H, Ca, Al, DIC, SO4, I, pH);
 
    number x[6];
    x[0] = pH; x[1] = DIC; x[2] = I; x[3] = Ca; x[4] = Al; x[5] = SO4;
 
    Kd = m_cloud->interpolate(x, m_entry);
  }
 
  template <int refDim>
  inline void evaluate(number vKd[],
             const MathVector<dim> vGlobIP[],
             number time, int si,
             GridObject* elem,
             const MathVector<dim> vCornerCoords[],
             const MathVector<refDim> vLocIP[],
             const size_t nip,
             LocalVector* u,
             const MathMatrix<refDim, dim>* vJT = NULL) const
  {
    std::vector<number> vH(nip);
    (*m_spH)  (&vH[0],   vGlobIP, time, si, elem, vCornerCoords, vLocIP, nip, u, vJT);
    std::vector<number> vS(nip);
    (*m_spS)  (&vS[0],   vGlobIP, time, si, elem, vCornerCoords, vLocIP, nip, u, vJT);
    std::vector<number> vCa(nip);
    (*m_spCa) (&vCa[0],  vGlobIP, time, si, elem, vCornerCoords, vLocIP, nip, u, vJT);
    std::vector<number> vDIC(nip);
    (*m_spDIC)(&vDIC[0], vGlobIP, time, si, elem, vCornerCoords, vLocIP, nip, u, vJT);
    std::vector<number> vSO4(nip);
    (*m_spSO4)(&vSO4[0], vGlobIP, time, si, elem, vCornerCoords, vLocIP, nip, u, vJT);
    std::vector<number> vAl(nip);
    (*m_spAl) (&vAl[0],  vGlobIP, time, si, elem, vCornerCoords, vLocIP, nip, u, vJT);
 
    for(size_t ip = 0; ip < nip; ++ip)
    {
      number I, pH;
      compute_I_pH_m3(vS[ip], vH[ip], vCa[ip], vAl[ip], vDIC[ip], vSO4[ip], I, pH);
 
      number x[6];
      x[0] = pH; x[1] = vDIC[ip]; x[2] = I; x[3] = vCa[ip]; x[4] = vAl[ip]; x[5] = vSO4[ip];
 
      vKd[ip] = m_cloud->interpolate(x, m_entry);
    }
  }
 
  template <int refDim>
  void eval_and_deriv(number vKd[],
            const MathVector<dim> vGlobIP[],
            number time, int si,
            GridObject* elem,
            const MathVector<dim> vCornerCoords[],
            const MathVector<refDim> vLocIP[],
            const size_t nip,
            LocalVector* u,
            bool bDeriv,
            int s,
            std::vector<std::vector<number> > vvvDeriv[],
            const MathMatrix<refDim, dim>* vJT = NULL) const
  {
    const number* vH   = m_spH  ->values(s);
    const number* vS   = m_spS  ->values(s);
    const number* vCa  = m_spCa ->values(s);
    const number* vDIC = m_spDIC->values(s);
    const number* vSO4 = m_spSO4->values(s);
    const number* vAl  = m_spAl ->values(s);
 
    for(size_t ip = 0; ip < nip; ++ip)
    {
      number I, pH;
      compute_I_pH_m3(vS[ip], vH[ip], vCa[ip], vAl[ip], vDIC[ip], vSO4[ip], I, pH);
 
      number x[6];
      x[0] = pH; x[1] = vDIC[ip]; x[2] = I; x[3] = vCa[ip]; x[4] = vAl[ip]; x[5] = vSO4[ip];
 
      vKd[ip] = m_cloud->interpolate(x, m_entry);
    }
    
    if (bDeriv)
    {
      // clear all derivative values
      this->set_zero(vvvDeriv, nip);
 
      // Derivatives w.r.t. H
      if (m_spDH.valid() && !m_spDH->zero_derivative())
        for(size_t ip = 0; ip < nip; ++ip)
          for(size_t fct = 0; fct < m_spDH->num_fct(); ++fct)
          {
            const number* vDH = m_spDH->deriv(s, ip, fct);
            const size_t commonFct = this->input_common_fct(0, fct);
            for(size_t sh = 0; sh < this->num_sh(commonFct); ++sh)
            {
              number eps = 1e-3*vH[ip];
 
              number I, pH;
              compute_I_pH_m3(vS[ip], vH[ip] + eps, vCa[ip], vAl[ip], vDIC[ip], vSO4[ip], I, pH);
 
              number x[6];
              x[0] = pH; x[1] = vDIC[ip]; x[2] = I; x[3] = vCa[ip]; x[4] = vAl[ip]; x[5] = vSO4[ip];
              number Kd1 = m_cloud->interpolate(x, m_entry);
              
              compute_I_pH_m3(vS[ip], vH[ip] - eps, vCa[ip], vAl[ip], vDIC[ip], vSO4[ip], I, pH);
 
              x[0] = pH; x[1] = vDIC[ip]; x[2] = I; x[3] = vCa[ip]; x[4] = vAl[ip]; x[5] = vSO4[ip];
              number Kd2 = m_cloud->interpolate(x, m_entry);
 
              vvvDeriv[ip][commonFct][sh] += vDH[sh]*(Kd1 - Kd2)/(2*eps);
            }
          }
 
      // Derivatives w.r.t. S
      if (m_spDS.valid() && !m_spDS->zero_derivative())
        for(size_t ip = 0; ip < nip; ++ip)
          for(size_t fct = 0; fct < m_spDS->num_fct(); ++fct)
          {
            const number* vDS = m_spDS->deriv(s, ip, fct);
            const size_t commonFct = this->input_common_fct(1, fct);
            for(size_t sh = 0; sh < this->num_sh(commonFct); ++sh)
            {
              number eps1, eps2;
              if (vS[ip] > 0.0)
                eps1 = eps2 = 1e-3*vS[ip];
              else {
                eps1 = 1e-3;
                eps2 = 0.0;
              }
 
              number I, pH;
              compute_I_pH_m3(vS[ip] + eps1, vH[ip], vCa[ip], vAl[ip], vDIC[ip], vSO4[ip], I, pH);
              
              number x[6];
              x[0] = pH; x[1] = vDIC[ip]; x[2] = I; x[3] = vCa[ip]; x[4] = vAl[ip]; x[5] = vSO4[ip];
              number Kd1 = m_cloud->interpolate(x, m_entry);
              
              compute_I_pH_m3(vS[ip] - eps2, vH[ip], vCa[ip], vAl[ip], vDIC[ip], vSO4[ip], I, pH);
              
              x[0] = pH; x[1] = vDIC[ip]; x[2] = I; x[3] = vCa[ip]; x[4] = vAl[ip]; x[5] = vSO4[ip];
              number Kd2 = m_cloud->interpolate(x, m_entry);
              
              vvvDeriv[ip][commonFct][sh] += vDS[sh]*(Kd1 - Kd2)/(eps1 + eps2); 
            }
          }
 
      // Derivatives w.r.t. Ca
      if (m_spDCa.valid() && !m_spDCa->zero_derivative())
        for(size_t ip = 0; ip < nip; ++ip)
          for(size_t fct = 0; fct < m_spDCa->num_fct(); ++fct)
          {
            const number* vDCa = m_spDCa->deriv(s, ip, fct);
            const size_t commonFct = this->input_common_fct(2, fct);
            for(size_t sh = 0; sh < this->num_sh(commonFct); ++sh)
            {
              number eps = 1e-3*vCa[ip];
 
              number I, pH;
              compute_I_pH_m3(vS[ip], vH[ip], vCa[ip] + eps, vAl[ip], vDIC[ip], vSO4[ip], I, pH);
 
              number x[6];
              x[0] = pH; x[1] = vDIC[ip]; x[2] = I; x[3] = vCa[ip] + eps; x[4] = vAl[ip]; x[5] = vSO4[ip];
              number Kd1 = m_cloud->interpolate(x, m_entry);
              
              compute_I_pH_m3(vS[ip], vH[ip], vCa[ip] - eps, vAl[ip], vDIC[ip], vSO4[ip], I, pH);
 
              x[0] = pH; x[1] = vDIC[ip]; x[2] = I; x[3] = vCa[ip] - eps; x[4] = vAl[ip]; x[5] = vSO4[ip];
              number Kd2 = m_cloud->interpolate(x, m_entry);
 
              vvvDeriv[ip][commonFct][sh] += vDCa[sh]*(Kd1 - Kd2)/(2*eps);
            }
          }
      
      // Derivatives w.r.t. DIC
      if (m_spDDIC.valid() && !m_spDDIC->zero_derivative())
        for(size_t ip = 0; ip < nip; ++ip)
          for(size_t fct = 0; fct < m_spDDIC->num_fct(); ++fct)
          {
            const number* vDDIC = m_spDDIC->deriv(s, ip, fct);
            const size_t commonFct = this->input_common_fct(3, fct);
            for(size_t sh = 0; sh < this->num_sh(commonFct); ++sh)
            {
              number eps = 1e-3*vDIC[ip];
 
              number I, pH;
              compute_I_pH_m3(vS[ip], vH[ip], vCa[ip], vAl[ip], vDIC[ip] + eps, vSO4[ip], I, pH);
 
              number x[6];
              x[0] = pH; x[1] = vDIC[ip] + eps; x[2] = I; x[3] = vCa[ip]; x[4] = vAl[ip]; x[5] = vSO4[ip];
              number Kd1 = m_cloud->interpolate(x, m_entry);
              
              compute_I_pH_m3(vS[ip], vH[ip], vCa[ip], vAl[ip], vDIC[ip] - eps, vSO4[ip], I, pH);
 
              x[0] = pH; x[1] = vDIC[ip] - eps; x[2] = I; x[3] = vCa[ip]; x[4] = vAl[ip]; x[5] = vSO4[ip];
              number Kd2 = m_cloud->interpolate(x, m_entry);
 
              vvvDeriv[ip][commonFct][sh] += vDDIC[sh]*(Kd1 - Kd2)/(2*eps);
            }
          }
 
      // Derivatives w.r.t. SO4
      if (m_spDSO4.valid() && !m_spDSO4->zero_derivative())
        for(size_t ip = 0; ip < nip; ++ip)
          for(size_t fct = 0; fct < m_spDSO4->num_fct(); ++fct)
          {
            const number* vDSO4 = m_spDSO4->deriv(s, ip, fct);
            const size_t commonFct = this->input_common_fct(4, fct);
            for(size_t sh = 0; sh < this->num_sh(commonFct); ++sh)
            {
              number eps1, eps2;
              if (vSO4[ip] > 0.0)
                eps1 = eps2 = 1e-3*vSO4[ip];
              else {
                eps1 = 1e-3;
                eps2 = 0.0;
              }
 
              number I, pH;
              compute_I_pH_m3(vS[ip], vH[ip], vCa[ip], vAl[ip], vDIC[ip], vSO4[ip] + eps1, I, pH);
 
              number x[6];
              x[0] = pH; x[1] = vDIC[ip]; x[2] = I; x[3] = vCa[ip]; x[4] = vAl[ip]; x[5] = vSO4[ip] + eps1;
              number Kd1 = m_cloud->interpolate(x, m_entry);
              
              compute_I_pH_m3(vS[ip], vH[ip], vCa[ip], vAl[ip], vDIC[ip], vSO4[ip] - eps2, I, pH);
 
              x[0] = pH; x[1] = vDIC[ip]; x[2] = I; x[3] = vCa[ip]; x[4] = vAl[ip]; x[5] = vSO4[ip] - eps2;
              number Kd2 = m_cloud->interpolate(x, m_entry);
 
              vvvDeriv[ip][commonFct][sh] += vDSO4[sh]*(Kd1 - Kd2)/(eps1 + eps2);
            }
          }
 
      // Derivatives w.r.t. Al
      if (m_spDAl.valid() && !m_spDAl->zero_derivative())
        for(size_t ip = 0; ip < nip; ++ip)
          for(size_t fct = 0; fct < m_spDAl->num_fct(); ++fct)
          {
            const number* vDDIC = m_spDAl->deriv(s, ip, fct);
            const size_t commonFct = this->input_common_fct(5, fct);
            for(size_t sh = 0; sh < this->num_sh(commonFct); ++sh)
            {
              number eps1, eps2;
                    if (vAl[ip] > 0.0)
                eps1 = eps2 = 1e-3*vAl[ip];
              else {
                eps1 = 1e-3;
                eps2 = 0.0;
              }
 
              number I, pH;
              compute_I_pH_m3(vS[ip], vH[ip], vCa[ip], vAl[ip] + eps1, vDIC[ip], vSO4[ip], I, pH);
 
              number x[6];
              x[0] = pH; x[1] = vDIC[ip]; x[2] = I; x[3] = vCa[ip]; x[4] = vAl[ip] + eps1; x[5] = vSO4[ip];
              number Kd1 = m_cloud->interpolate(x, m_entry);
              
              compute_I_pH_m3(vS[ip], vH[ip], vCa[ip], vAl[ip] - eps2, vDIC[ip], vSO4[ip], I, pH);
 
              x[0] = pH; x[1] = vDIC[ip]; x[2] = I; x[3] = vCa[ip]; x[4] = vAl[ip] - eps2; x[5] = vSO4[ip];
              number Kd2 = m_cloud->interpolate(x, m_entry);
 
              vvvDeriv[ip][commonFct][sh] += vDDIC[sh]*(Kd1 - Kd2)/(eps1 + eps2);
            }
          }
    } 
  }
 
  protected:
 
  Cloud *m_cloud;
  SmartPtr<CplUserData<number, dim> > m_spH;
  SmartPtr<DependentUserData<number, dim> > m_spDH;
  SmartPtr<CplUserData<number, dim> > m_spS;
  SmartPtr<DependentUserData<number, dim> > m_spDS;
  SmartPtr<CplUserData<number, dim> > m_spCa;
  SmartPtr<DependentUserData<number, dim> > m_spDCa;
  SmartPtr<CplUserData<number, dim> > m_spDIC;
  SmartPtr<DependentUserData<number, dim> > m_spDDIC;
  SmartPtr<CplUserData<number, dim> > m_spSO4;
  SmartPtr<DependentUserData<number, dim> > m_spDSO4;
  SmartPtr<CplUserData<number, dim> > m_spAl;
  SmartPtr<DependentUserData<number, dim> > m_spDAl;
  const number m_entry;
};
  
} // end namespace smile
} // end namespace ug
 
#endif