Propagation.H

Go to the documentation of this file.

 
#ifndef _PROPAGATION_H_
#define _PROPAGATION_H_
 
#include "../Core/CRPacket.H"
#include "../Core/FieldQty.H"
#include "../Gas/Gas.H"
#include "../IO/ParmParser.H"
#include "../Utils/Types.H"
#include "../Utils/Vec3.H"
 
namespace criptic {
 
  typedef enum {
    noFieldQty = 0,   
    needFieldQty,     
    needFieldQtyGrad  
  } FieldQtyNeedType;
  
  namespace propagation {
    
    typedef struct {
      Real kPar;         
      Real kPerp;        
      Real kPP;          
#ifdef TRACK_PITCH_ANGLE
      Real kMu;          
#else
      Real vStr;         
#endif
      Real dkPP_dp;      
      Real dvStr_dp;     
      RealVec kParGrad;  
      RealVec kPerpGrad; 
#ifdef TRACK_PITCH_ANGLE
      Real dkMudMu;      
#else
      RealVec vStrGrad;  
#endif      
    } PropagationData;
 
    class Propagation {
 
    public:
 
      Propagation(const ParmParser& pp) {
        if (!pp.query("cr.enforceBohmLimit", enforceBohmLimit)) {
          enforceBohmLimit = 0;
        }
        if (!pp.query("cr.kParMax", kParMax)) kParMax = 0;
        if (!pp.query("cr.kPerpMax", kPerpMax)) kPerpMax = 0;
        if (!pp.query("cr.kPPMax", kPPMax)) kPPMax = 0;
        if (!pp.query("cr.vDriftMax", vDriftMax)) vDriftMax = 1.0;
      }
      
      virtual ~Propagation() { };
    
      virtual FieldQtyNeedType fieldQtyNeed() const = 0;
      
      virtual PropagationData
      operator()(const RealVec& x,
         const Real t,
         const gas::GasData& gd,
         const CRPacket& packet,
         const FieldQty& qty,
         const FieldQtyGrad& qtyGrad) const = 0;
 
      void applyLimits(
            const gas::GasData& gd,
          const CRPacket& packet,
            PropagationData& pd) const
      {
            if (enforceBohmLimit) {
                  const Real Bmag = gd.B.mag();
                  const Real kBohm = packet.rG(Bmag) * 
                        packet.v() * constants::c / 3;
                  if (pd.kPar < kBohm) {
                        pd.kPar = kBohm;
                        pd.kParGrad = -kBohm * gd.BGrad.contract2(gd.B) / (Bmag*Bmag);
                  }
            }
            if (vDriftMax > 0.0) {
                  const Real facPar = pd.kParGrad.mag() / (vDriftMax * constants::c);
                  if (facPar > 1.0) {
                        pd.kParGrad /= facPar;
                        pd.kPerpGrad /= facPar;
                  }
                  const Real facPerp = pd.kPerpGrad.mag() / (vDriftMax * constants::c);
                  if (facPerp > 1.0) {
                        pd.kParGrad /= facPerp;
                        pd.kPerpGrad /= facPerp;
                  }
                  const Real facStr = fabs(pd.vStr + packet.p * pd.dvStr_dp /3) /
                        (vDriftMax * constants::c);
                  if (facStr > 1.0) {
                        pd.vStr /= facStr;
                        pd.dvStr_dp /= facStr;
                  }
            }
            if (kParMax > 0 && pd.kPar > kParMax) {
                  pd.kPar = kParMax;
                  pd.kParGrad = zeroVec;
            }
            if (kPerpMax > 0 && pd.kPerp > kPerpMax) {
                  pd.kPerp = kPerpMax;
                  pd.kPerpGrad = zeroVec;
            }
            if (kPPMax > 0 && pd.kPP > kPPMax) {
                  pd.kPP = kPPMax;
                  pd.dkPP_dp = 0;
            }
      }
      
 
    private:
 
      int enforceBohmLimit;    
      Real kParMax;            
      Real kPerpMax;           
      Real kPPMax;             
      Real vDriftMax;          
    };
  }
}
 
#endif
// _PROPAGATION_H_