AdvancePacket.H

Go to the documentation of this file.

 
#ifndef _ADVANCEPACKET_H_
#define _ADVANCEPACKET_H_
 
#include "CRPacket.H"
#include "FieldQty.H"
#include "Geometry.H"
#include "../Gas/Gas.H"
#include "../Losses/Losses.H"
#include "../Propagation/Propagation.H"
#include "../Utils/RealTensor2.H"
#include "../Utils/RKF.H"
#include "../Utils/RngThread.H"
#include "../Utils/Types.H"
#include "../Utils/Vec3.H"
#include <sstream>
 
namespace criptic {
 
  criptic::Real advancePacket(const criptic::Real tStart,
                  const criptic::Real tStop,
                  const criptic::Geometry& geom,
                  const criptic::gas::Gas& gasBG,
                  const criptic::propagation::Propagation& prop,
                  const criptic::Losses& loss,
                  const criptic::RngThread& rng,
                  const criptic::FieldQty& qty,
                  const criptic::FieldQtyGrad& qtyGrad,
                  const criptic::Real cStep,
                  const criptic::Real errTol,
                  const criptic::Real wFracMin,
                  const criptic::Real pMin,
                  const criptic::Real TMin,
                  const criptic::Real packetDtMin,
                  const criptic::Real& h,
                  criptic::RealVec& x,
                  criptic::CRPacket& packet,
                  criptic::Real& dt,
                  std::vector<RealVec>& secondaryPositions,
                  std::vector<CRPacket>& secondaryPackets) {
 
    // Set time to starting time for this packet, and initialize
    // minimum time step recorder to numerical infinity
    Real t = std::max(packet.tInj, tStart);
    dt = maxReal;
    
    // Loop until update is complete
    bool packetDone = false;
    while (!packetDone) {
 
      // Step 1: grab the gas state at the start of the step, and
      // compute the TNB basis vectors from it
      gas::GasData gd = gasBG.gasData(x, t);
      TNBBasis b = TNBVectors(gd.B, gd.BGrad);
 
      // Step 2: grab the propagation parameters at the starting
      // position and time
      propagation::PropagationData pd =
    prop(x, t, gd, packet, qty, qtyGrad);
      prop.applyLimits(gd, packet, pd);
 
      // Step 3: compute momentum loss rate, catastrophic loss rate,
      // and properties of secondary creation channels
      Real dpdtCts, lossRate, wSec;
      Losses::MechMatrix<Real> pSec;
      loss.computeLoss(gd, packet, dpdtCts, lossRate, wSec, pSec);
 
      // Step 4: construct the drift vector in the TNB basis,
      // co-moving with the gas. The mathematical relationships here
      // are:
      //
      // For the case with streaming, not following pitch angle:
      // drift[0] = dKpar / dxT + vStr + (p/3) dvStr/dp
      // drift[1] = dKperp / dxN
      // drift[2] = dKperp / dxB
      // pDrift = dKpp / dp + (2/p) kPP - \dot{p}_cts -
      //      (p/3) * (del . v + grad(vStr) . eT + vStr del . eT),
      //
      // For the case following pitch angle evolution:
      // drift[0] = dKpar / dxT + mu * v
      // drift[1] = dKperp / dxN
      // drift[2] = dKperp / dxB
      // pDrift = dKpp / dp + (2/p) kPP - \dot{p}_cts -
      //      (p/3) * (del . v)
      // muDrift = -2 * mu * kMu + (1 - mu^2) dkMu / dMu
      //
      // where d/dxT indicates derivative in the T direction, and
      // similarly for d/dxN and d/dxB. We evalute these derivatives
      // using the relationship that, for an arbitrary scalar q, we
      // have dq/dxT = grad(q) . eT, and similarly for the N and B
      // directions. We also evaluate the final term del . eT, the
      // divergence of the eT vector, using the identity
      // div(eT) = -eT_i eT_j (dB_i / dx_j) / |B|
      RealVec xDrift;
#ifndef TRACK_PITCH_ANGLE
      xDrift[0] = pd.vStr + packet.p * pd.dvStr_dp / 3.0
    + pd.kParGrad.dot(b.eT);
#else
      xDrift[0] = packet.mu * packet.v() * constants::c + pd.kParGrad.dot(b.eT);
#endif
      xDrift[1] = pd.kPerpGrad.dot(b.eN);
      xDrift[2] = pd.kPerpGrad.dot(b.eB);
#ifdef TRACK_PITCH_ANGLE
      Real pDrift = pd.dkPP_dp + 2 * pd.kPP/packet.p - dpdtCts
    - packet.p / 3 * gd.vGrad.trace();
      Real muDrift = -2 * packet.mu * pd.kMu +
    (1.0 - packet.mu * packet.mu) * pd.dkMudMu;
#else
      Real pDrift = pd.dkPP_dp + 2 * pd.kPP/packet.p - dpdtCts
    - packet.p / 3 *
    (gd.vGrad.trace() + pd.vStrGrad.dot(b.eT) -
     pd.vStr * gd.BGrad.contract(b.eT, b.eT) / gd.B.mag());
#endif
      
      // Step 5: compute the time step
      Real dx = gd.dx;                    // Gas length scale
      if (h != 0) dx = std::min(dx, h);   // Packet length scale
      Real xDriftMag = xDrift.mag();
      Real dtAdv = (xDriftMag != 0.0 && dx != 0) ?
    cStep * dx / xDriftMag : maxReal;
      Real kMag = std::max(pd.kPar, pd.kPerp);
      Real dtDiff = (kMag != 0.0 && dx != 0.0) ?
    cStep * dx * dx / kMag : maxReal;
      Real dtPAdv = pDrift != 0.0 ?
    cStep * packet.p / std::abs(pDrift) : maxReal;
      Real dtPP =
    pd.kPP != 0.0 ? cStep * packet.p * packet.p / pd.kPP : maxReal;
      Real dtLoss = lossRate != 0.0 ? cStep / lossRate : maxReal;
#ifdef TRACK_PITCH_ANGLE
      Real dtMuAdv = fabs(muDrift) != 0.0 ?
    0.02 / fabs(muDrift) : maxReal;
      Real dtMuDiff = pd.kMu != 0.0 ?
    0.02 * cStep / pd.kMu : maxReal;
#endif
      // Safety check: make sure we have at least one valid time
      // scale; issue a warning if not, but do not halt, since in some
      // problems this will occur only very early in the calculation
      // when the number of packets is too small to compute smoothing
      // lengths or gradients, but the problem will resolve as soon as
      // the number of packets becomes reasonable
      if (dtAdv == maxReal && dtDiff == maxReal &&
      dtPAdv == maxReal && dtPP == maxReal &&
      dtPP == maxReal && dtLoss == maxReal
#ifdef TRACK_PITCH_ANGLE
      && dtMuAdv == maxReal && dtMuDiff == maxReal
#endif
      ) {
    std::cout << "AdvancePacket: warning: unable to compute time "
          << "step size; note that in problems where the gas "
          << "does not contain a characteristic length scale, "
          << "the length scale dx must be set by hand in order "
          << "to compute valid spatial drift or diffusion "
          << "timescales; calculation will continue, but time "
          << "step will not be self-consistently computed for "
          << "this packet advance";
      }
      
      Real dtStep =  1.0 / (1.0/dtAdv + 1.0/dtDiff + 1.0/dtPAdv
                + 1.0/dtPP + 1.0/dtLoss
#ifdef TRACK_PITCH_ANGLE
                + 1.0/dtMuAdv + 1.0/dtMuDiff
#endif
                );
      
      // Check for unacceptably small packet time step
      if (!(dtStep >= packetDtMin)) {
    std::stringstream ss;
    ss << "AdvancePacket: packet dt = " << dtStep
       << " < min dt = " << packetDtMin
       << "; individual time steps are: "
       << "dtAdv = " << dtAdv
       << ", dtDiff = " << dtDiff
       << ", dtPAdv = " << dtPAdv
       << ", dtPP = " << dtPP
       << ", dtLoss = " << dtLoss
#ifdef TRACK_PITCH_ANGLE
       << ", dtMuAdv = " << dtMuAdv
       << ", dtMuDiff = " << dtMuDiff
#endif
       << "; packet x = " << x
       << ", p = " << packet.p
       << ", type = " << packet.type
       << "; propagation parameters "
#ifdef TRACK_PITCH_ANGLE
       << "kMu = " << pd.kMu
       << ", dkMu_dMu = " << pd.dkMudMu << ", "
#else
       << "vStr = " << pd.vStr
       << ", vStrGrad = " << pd.vStrGrad
       << ", dvStr_dp = " << pd.dvStr_dp
#endif
       << ", kPar = " << pd.kPar
       << ", kPerp = " << pd.kPerp
       << ", kPP = " << pd.kPP
       << ", kParGrad = " << pd.kParGrad
       << ", kPerpGrad = " << pd.kPerpGrad
       << ", dkPP_dp = " << pd.dkPP_dp;
    MPIUtil::haltRun(ss.str(), errTimeStepTooSmall);
      }
      if (dt > dtStep) {
    dt = dtStep;   // Record minimum internal time step
      }
      if (t + dtStep >= tStop) {
    dtStep = tStop - t;
    packetDone = true;
      }
 
      // Step 6: update packet grammage, momentum, and weight; for
      // pitch angle scattering, also update the pitch angle
      packet.gr += gd.den * packet.v() * constants::c * dtStep;
      Real wFac = std::exp( -lossRate * dtStep);
      packet.w *= wFac;
      packet.p += pDrift * dtStep +
    std::sqrt(2 * pd.kPP * dtStep) * rng.normal();
#ifdef TRACK_PITCH_ANGLE
      packet.mu += dtStep * muDrift +
    std::sqrt(2 * pd.kMu * (1.0 - packet.mu * packet.mu) * dtStep)
    * rng.normal();
      while (packet.mu > 2) packet.mu -= 2;
      while (packet.mu < -2) packet.mu += 2;
      if (packet.mu > 1) packet.mu = 2 - packet.mu;
      if (packet.mu < -1) packet.mu = -2 - packet.mu;
#endif
 
      // Step 7: create secondaries
      for (int i = 0; i < lossMech::nMech; i++) {
    for (int j = 0; j < partTypes::nPartType; j++) {
      if (pSec[i][j] == 0.0) continue;
      if (rng.uniform() < pSec[i][j] * (1 - wFac)) {
        secondaryPositions.push_back(x);
        secondaryPackets.
          push_back(loss.drawSecondary(gd, packet, t, wSec / wFac,
                       (lossMech::mechType) i,
                       (partTypes::pType) j));
      }
    }
      }
 
      // Step 8: compute effective velocities for position update
      Real vT = xDrift[0]
    + std::sqrt( 2 * pd.kPar / dtStep) * rng.normal();
      Real vN = xDrift[1] +
    std::sqrt( 2 * pd.kPerp / dtStep ) * rng.normal();
      Real vB = xDrift[2] +
    + std::sqrt( 2 * pd.kPerp / dtStep ) * rng.normal();
      Real vMag = gd.v.mag() + std::sqrt( vT*vT + vN*vN + vB*vB );
 
      // Step 9: update packet position; this may require several
      // sub-cycles to ensure field lines are followed within the
      // required tolerance
      {
 
    // Initial number of substeps; for gas descriptors that have
    // ghost zones, ensure that this is small enough that we
    // cannot move outside the ghost region in a single time step
    int nDiv = 1;
    if (gasBG.dxGhost() > 0.0) {
      while (dtStep * vMag > nDiv * gasBG.dxGhost()) nDiv *= 2;
    }
 
    // Save starting position and time
    RealVec xSave = x;
 
    // Try to advance through substeps
    bool success = false;
    while (!success) {
 
      // Step substep size for this attempt
      Real dtTmp = dtStep / nDiv;
      success = true;
 
      // Save position at beginning of substep
      RealVec xInit = x;
 
      // Advance through substeps; at each substep, recompute the
      // comoving TNB frame, then move the packet, including applying
      // periodic or reflecting boundary conditions as needed; we then
      // check if the accumulated error is acceptably small at the
      // end
      Real err2 = 0.0;
      for (int n=0; n<nDiv; n++) {
 
        // Stage 1
        RealVec vCom;
        gasBG.frame(x, t + dtTmp * n, vCom, b);
        RealVec k1 = vCom + vT*b.eT + vN*b.eN + vB*b.eB;
        x = xInit + dtTmp * rkf::ah[0] * k1;
        geom.applyBC(x
#ifdef TRACK_PITCH_ANGLE
            , packet
#endif
        );
 
        // Stage 2
        gasBG.frame(x, t + dtTmp * (n + rkf::ah[0]), vCom, b);
        RealVec k2 = vCom + vT*b.eT + vN*b.eN + vB*b.eB;
        x = xInit + dtTmp * (rkf::b3[0] * k1 + rkf::b3[1] * k2);
        geom.applyBC(x
#ifdef TRACK_PITCH_ANGLE
            , packet
#endif      
        );
 
        // Stage 3
        gasBG.frame(x, t + dtTmp * (n + rkf::ah[1]), vCom, b);
        RealVec k3 = vCom + vT*b.eT + vN*b.eN + vB*b.eB;
        x = xInit + dtTmp * (rkf::b4[0] * k1 + rkf::b4[1] * k2 +
                 rkf::b4[2] * k3);
        geom.applyBC(x
#ifdef TRACK_PITCH_ANGLE
            , packet
#endif
        );
 
        // Stage 4
        gasBG.frame(x, t + dtTmp * (n + rkf::ah[2]), vCom, b);
        RealVec k4 = vCom + vT*b.eT + vN*b.eN + vB*b.eB;
        x = xInit + dtTmp * (rkf::b5[0] * k1 + rkf::b5[1] * k2 +
                 rkf::b5[2] * k3 + rkf::b5[3] * k4);
        geom.applyBC(x
#ifdef TRACK_PITCH_ANGLE
            , packet
#endif
        );
 
        // Stage 5
        gasBG.frame(x, t + dtTmp * (n + rkf::ah[3]), vCom, b);
        RealVec k5 = vCom + vT*b.eT + vN*b.eN + vB*b.eB;
        x = xInit + dtTmp * (rkf::b6[0] * k1 + rkf::b6[1] * k2 +
                 rkf::b6[2] * k3 + rkf::b6[3] * k4 +
                 rkf::b6[4] * k5);
        geom.applyBC(x
#ifdef TRACK_PITCH_ANGLE
            , packet
#endif
        );
 
        // Stage 6
        gasBG.frame(x, t + dtTmp * (n + rkf::ah[4]), vCom, b);
        RealVec k6 = vCom + vT*b.eT + vN*b.eN + vB*b.eB;
        x = xInit + dtTmp * (rkf::c1 * k1 + 
                 rkf::c3 * k3 +
                 rkf::c4 * k4 +
                 rkf::c5 * k5 +
                 rkf::c6 * k6);
        geom.applyBC(x
#ifdef TRACK_PITCH_ANGLE
            , packet
#endif
        );
 
        // Compute error on this step, add to total
        err2 += (dtTmp/dx) * (dtTmp/dx) *
          (rkf::ec[1] * k1 +
           rkf::ec[3] * k3 +
           rkf::ec[4] * k4 +
           rkf::ec[5] * k5 +
           rkf::ec[6] * k6).mag2();
 
        // Check if we have exceeded error budget; if so, restore
        // starting position and time, increase number of
        // divisions, and try again
        if (err2 > errTol * errTol) {
          x = xSave;
          nDiv *= 2;
          success = false;
          break;
        }
 
        // Check to ensure that packet is still in domain; if not,
        // cease update and flag packet for deletion by flipping its
        // weight negative
        if (!geom.contains(x)) {
          packet.w = -packet.w;
          return t;
        }
      }
    }
      }
 
      // Advance time
      t += dtStep;
 
      // Step 10: check if packet has fallen below minimum weight or
      // momentum; if so, flag for deletion by flipping weight
      // negative, and return
      if (packet.w / packet.wInj < wFracMin ||
      packet.p < pMin ||
      packet.T() < TMin) {
    packet.w = -packet.w;
    return t;
      }
    }
 
    // Return time
    return t;
  }
 
}
 
#endif
// _ADANCEPACKET_H_