FieldQty.H

Go to the documentation of this file.

 
#ifndef _FIELDQTY_H_
#define _FIELDQTY_H_
 
#include <array>
#include <string>
#include <vector>
#include "CRPacket.H"
#include "../Definitions.H"
#include "../Utils/Constants.H"
#include "../Utils/RealTensor2.H"
#include "../Utils/Types.H"
#include "../Utils/Units.H"
#include "../Utils/Vec3.H"
 
namespace criptic {
 
  typedef enum {
    ndenIdx = 0,  
    presIdx,      
    edenIdx,      
    nFieldQty     
  } FieldQtyType;
  
  // FieldQty class
 
  class FieldQty {
 
  public:
 
    std::array<Real,nFieldQty> q;  
    // Names and units of field quantities -- used in ASCII output
    static constexpr int qtyNameLen = 8; 
    static constexpr std::array<char[qtyNameLen],nFieldQty> qtyNames =
      {{ "nden", "pres", "eden" }};  
    static constexpr std::array<Real,nFieldQty> qtyUnits = {
      1.0, constants::mp_c2 / constants::kB, constants::mp_c2 / units::eV
    };  
    FieldQty() {
      for (int i = 0; i < nFieldQty; i++) q[i] = 0;
    }
    
    void setZero() {
      for (int i=0; i<nFieldQty; i++) q[i] = 0.0;
    }
 
    std::array<Real,nFieldQty> getWgt(const CRPacket& pkt) const {
      std::array<Real,nFieldQty> w;
      w[ndenIdx] = pkt.w;
      w[presIdx] = pkt.w * pkt.v() * pkt.p;
      w[edenIdx] = pkt.w * pkt.T();
      return w;
    }
    
    void addPacket(const CRPacket& pkt,
           const RealVec& r,
           const RealTensor2& h2Inv) {
      Real k = std::exp( -0.5 * h2Inv.contract(r, r) );
      q[ndenIdx] += k * pkt.w;
      q[presIdx] += k * pkt.w * pkt.v() * pkt.p;
      q[edenIdx] += k * pkt.w * pkt.T();
    }
 
    Real& operator [] (int i) {
      return q[i];
    }
    const Real& operator [] (int i) const {
      return q[i];
    }
    
    // Binary operations with two FieldQty objects
    
    FieldQty operator+(const FieldQty& a) const {
      FieldQty fq;
      for (int i=0; i<nFieldQty; i++) fq.q[i] = q[i] + a.q[i];
      return fq;
    }
    FieldQty operator-(const FieldQty& a) const {
      FieldQty fq;
      for (int i=0; i<nFieldQty; i++) fq.q[i] = q[i] - a.q[i];
      return fq;
    }
    std::array<Real,nFieldQty> operator/(const FieldQty& a) const {
      std::array<Real,nFieldQty> arr;
      for (int i=0; i<nFieldQty; i++) arr[i] = q[i] / a.q[i];
      return arr;
    }
 
    // Binary operations with scalars
    FieldQty operator*(const Real a) const {
      FieldQty fq;
      for (int i=0; i<nFieldQty; i++) fq.q[i] = q[i] * a;
      return fq;
    }
    FieldQty operator/(const Real a) const {
      FieldQty fq;
      for (int i=0; i<nFieldQty; i++) fq.q[i] = q[i] / a;
      return fq;
    }
    
    // Unary operations with scalars
    FieldQty& operator*=(const Real& a) {
      for (int i=0; i<nFieldQty; i++) q[i] *= a;
      return *this;
    }
    FieldQty& operator/=(const Real& a) {
      for (int i=0; i<nFieldQty; i++) q[i] /= a;
      return *this;
    }    
    
    // Unary operations with field quantities
    FieldQty& operator+=(const FieldQty& a) {
      for (int i=0; i<nFieldQty; i++) q[i] += a.q[i];
      return *this;
    }
    FieldQty& operator-=(const FieldQty& a) {
      for (int i=0; i<nFieldQty; i++) q[i] -= a.q[i];
      return *this;
    }
 
  }; // end class FieldQty
  
  inline criptic::FieldQty operator*(const criptic::Real a,
                     const criptic::FieldQty& fq) {
    FieldQty f;
    for (int i=0; i<nFieldQty; i++) f.q[i] = a * fq.q[i];
    return f;
  }
 
  inline criptic::FieldQty min(const criptic::FieldQty& q1,
                   const criptic::FieldQty& q2) {
    FieldQty f;
    for (int i=0; i<nFieldQty; i++)
      f.q[i] = std::min(q1.q[i], q2.q[i]);
    return f;
  }
 
  inline criptic::FieldQty max(const criptic::FieldQty& q1,
                   const criptic::FieldQty& q2) {
    FieldQty f;
    for (int i=0; i<nFieldQty; i++)
      f.q[i] = std::max(q1.q[i], q2.q[i]);
    return f;
  }
  
  // FieldQtyGrad class
 
  class FieldQtyGrad {
 
  public:
    
    std::array<RealVec,nFieldQty> qGrad; 
    FieldQtyGrad() {
      for (int i = 0; i < nFieldQty; i++) qGrad[i] = 0;
    }
 
    void setZero() {
      for (int i=0; i<nFieldQty; i++) qGrad[i] = 0.0;
    }
 
    void addPacket(const CRPacket& pkt,
           const RealVec& r,
           const RealTensor2& h2Inv) {
      RealVec hr = h2Inv * r;
      RealVec dk = -hr * std::exp( -hr.dot(r) / 2);
      qGrad[ndenIdx] += dk * pkt.w;
      qGrad[presIdx] += dk * pkt.w * pkt.v() * pkt.p;
      qGrad[edenIdx] += dk * pkt.w * pkt.T();
    }
    
    RealVec& operator [] (int i) {
      return qGrad[i];
    }
    const RealVec& operator [] (int i) const {
      return qGrad[i];
    }
 
    // Binary operations involving two FieldQtyGrad objects
    FieldQtyGrad operator+(const FieldQtyGrad& a) const {
      FieldQtyGrad fqg;
      for (int i=0; i<nFieldQty; i++)
    fqg.qGrad[i] = qGrad[i] + a.qGrad[i];
      return fqg;
    }
    FieldQtyGrad operator-(const FieldQtyGrad& a) const {
      FieldQtyGrad fqg;
      for (int i=0; i<nFieldQty; i++)
    fqg.qGrad[i] = qGrad[i] - a.qGrad[i];
      return fqg;
    }
    std::array<RealVec,nFieldQty> operator/(const FieldQtyGrad& a) const {
      std::array<RealVec,nFieldQty> arr;
      for (int i=0; i<nFieldQty; i++) arr[i] = qGrad[i] / a.qGrad[i];
      return arr;
    }
 
    // Binary operations with scalars
    FieldQtyGrad operator*(const Real a) const {
      FieldQtyGrad fqg;
      for (int i=0; i<nFieldQty; i++)
    fqg.qGrad[i] = qGrad[i] * a;
      return fqg;
    }
    FieldQtyGrad operator/(const Real a) const {
      FieldQtyGrad fqg;
      for (int i=0; i<nFieldQty; i++)
    fqg.qGrad[i] = qGrad[i] / a;
      return fqg;
    }
 
    // Binary operations with vectors
    FieldQtyGrad operator*(const RealVec& a) const {
      FieldQtyGrad fqg;
      for (int i=0; i<nFieldQty; i++)
    fqg.qGrad[i] = qGrad[i] * a;
      return fqg;
    }
    
    // Unary operations with scalars
    FieldQtyGrad& operator*=(const Real& a) {
      for (int i=0; i<nFieldQty; i++) qGrad[i] *= a;
      return *this;
    }
    FieldQtyGrad& operator/=(const Real& a) {
      for (int i=0; i<nFieldQty; i++) qGrad[i] /= a;
      return *this;
    }
    
    // Unary operations with field quantity gradients
    FieldQtyGrad& operator+=(const FieldQtyGrad& a) {
      for (int i=0; i<nFieldQty; i++) qGrad[i] += a.qGrad[i];
      return *this;
    }
    FieldQtyGrad& operator-=(const FieldQtyGrad& a) {
      for (int i=0; i<nFieldQty; i++) qGrad[i] -= a.qGrad[i];
      return *this;
    }
    
  }; // end FieldQtyGrad
 
  inline criptic::FieldQtyGrad operator*(const criptic::Real a,
                     const criptic::FieldQtyGrad& fqg) {
    FieldQtyGrad f;
    for (int i=0; i<nFieldQty; i++) f.qGrad[i] = a * fqg.qGrad[i];
    return f;
  }
 
  inline criptic::FieldQtyGrad operator*(const criptic::RealVec& a,
                     const criptic::FieldQtyGrad& fqg) {
    FieldQtyGrad f;
    for (int i=0; i<nFieldQty; i++) f.qGrad[i] = a * fqg.qGrad[i];
    return f;
  }
  
  inline criptic::FieldQtyGrad min(const criptic::FieldQtyGrad& q1,
                   const criptic::FieldQtyGrad& q2) {
    FieldQtyGrad f;
    for (int i=0; i<nFieldQty; i++)
      f[i] = min(q1[i], q2[i]);
    return f;
  }
 
  inline criptic::FieldQtyGrad max(const criptic::FieldQtyGrad& q1,
                   const criptic::FieldQtyGrad& q2) {
    FieldQtyGrad f;
    for (int i=0; i<nFieldQty; i++)
      f[i] = max(q1[i], q2[i]);
    return f;
  }
  
} // end namespace criptic
 
 
 
// Overload the >> and << operators to make it easy to do IO on field
// quantities
std::istream& operator>>(std::istream& is,
             criptic::FieldQty& fQ);
 
std::istream& operator>>(std::istream& is,
             criptic::FieldQtyGrad& fQG);
 
std::ostream& operator<<(std::ostream& os,
             const criptic::FieldQty& fQ);
 
std::ostream& operator<<(std::ostream& os,
             const criptic::FieldQtyGrad& fQG);
 
#endif
// _FIELDQTY_H_