src/lagrangian/molecularDynamics/molecule/mdTools/averageMDFields.H

   1 if (runTime.outputTime())
   2 {
   3     /*-----------------------------------------------------------------------*\
   4         Number density
   5     \*-----------------------------------------------------------------------*/
   6
   7     scalarField totalRhoN_sum(mesh.nCells(), 0.0);
   8
   9     forAll(allSpeciesRhoN, rN)
  10     {
  11         allSpeciesRhoN[rN].internalField() =
  12             allSpeciesN_RU[rN]
  13             /mesh.cellVolumes()
  14             /nAveragingSteps;
  15
  16         totalRhoN_sum += allSpeciesRhoN[rN].internalField();
  17     }
  18
  19     totalRhoN.internalField() = totalRhoN_sum;
  20
  21
  22     /*-----------------------------------------------------------------------*\
  23         Mass density
  24     \*-----------------------------------------------------------------------*/
  25
  26     scalarField totalRhoM_sum(mesh.nCells(), 0.0);
  27
  28     forAll(allSpeciesRhoM, rM)
  29     {
  30         allSpeciesRhoM[rM].internalField() =
  31             allSpeciesM_RU[rM]
  32             /mesh.cellVolumes()
  33             /nAveragingSteps;
  34
  35         totalRhoM_sum += allSpeciesRhoM[rM].internalField();
  36     }
  37
  38     totalRhoM.internalField() = totalRhoM_sum;
  39
  40     /*-----------------------------------------------------------------------*\
  41         Bulk velocity
  42     \*-----------------------------------------------------------------------*/
  43
  44     vectorField totalMomentum_sum(mesh.nCells(), vector::zero);
  45
  46     scalarField totalMass_sum(mesh.nCells(), 0.0);
  47
  48     forAll(allSpeciesVelocity, v)
  49     {
  50         // A check for 1/0 molecules is required.
  51
  52         vectorField& singleSpeciesVelocity
  53         (
  54             allSpeciesVelocity[v].internalField()
  55         );
  56
  57         forAll(singleSpeciesVelocity, sSV)
  58         {
  59             if (allSpeciesN_RU[v][sSV])
  60             {
  61                 singleSpeciesVelocity[sSV] =
  62                     allSpeciesVelocitySum_RU[v][sSV]
  63                     /allSpeciesN_RU[v][sSV];
  64
  65                 totalMomentum_sum[sSV] +=
  66                     allSpeciesM_RU[v][sSV]
  67                     /allSpeciesN_RU[v][sSV]
  68                     *allSpeciesVelocitySum_RU[v][sSV];
  69
  70                 totalMass_sum[sSV] += allSpeciesM_RU[v][sSV];
  71             }
  72             else
  73             {
  74                 singleSpeciesVelocity[sSV] = vector::zero;
  75             }
  76         }
  77     }
  78
  79     forAll(totalVelocity.internalField(), tV)
  80     {
  81         if (totalMass_sum[tV] > VSMALL)
  82         {
  83             totalVelocity.internalField()[tV] =
  84                 totalMomentum_sum[tV]
  85                 /totalMass_sum[tV];
  86         }
  87         else
  88         {
  89             totalVelocity.internalField()[tV] =
  90             vector::zero;
  91         }
  92     }
  93
  94     /*-----------------------------------------------------------------------*\
  95         Kinetic temperature
  96     \*-----------------------------------------------------------------------*/
  97
  98     scalarField totalTemperatureVTerms_sum(mesh.nCells(), 0.0);
  99
 100     scalarField totalN_sum(mesh.nCells(), 0.0);
 101
 102     forAll(allSpeciesTemperature, t)
 103     {
 104         // A check for 1/0 molecules is required.
 105
 106         scalarField& singleSpeciesTemp
 107         (
 108             allSpeciesTemperature[t].internalField()
 109         );
 110
 111         forAll(singleSpeciesTemp, sST)
 112         {
 113             if (allSpeciesN_RU[t][sST])
 114             {
 115                 singleSpeciesTemp[sST] =
 116                     allSpeciesM_RU[t][sST]
 117                     /allSpeciesN_RU[t][sST]
 118                     /(3.0 * moleculeCloud::kb * allSpeciesN_RU[t][sST])
 119                    *(
 120                         allSpeciesVelocityMagSquaredSum_RU[t][sST]
 121                         -
 122                         (
 123                             allSpeciesVelocitySum_RU[t][sST]
 124                             &
 125                             allSpeciesVelocitySum_RU[t][sST]
 126                         )
 127                         /allSpeciesN_RU[t][sST]
 128                     );
 129
 130                 totalTemperatureVTerms_sum[sST] +=
 131                     allSpeciesM_RU[t][sST]
 132                    /allSpeciesN_RU[t][sST]
 133                    *(
 134                         allSpeciesVelocityMagSquaredSum_RU[t][sST]
 135                       -
 136                         (
 137                             allSpeciesVelocitySum_RU[t][sST]
 138                             &
 139                             allSpeciesVelocitySum_RU[t][sST]
 140                         )
 141                         /allSpeciesN_RU[t][sST]
 142                     );
 143
 144                 totalN_sum[sST] += allSpeciesN_RU[t][sST];
 145             }
 146             else
 147             {
 148                 singleSpeciesTemp[sST] = 0.0;
 149             }
 150         }
 151     }
 152
 153     forAll(totalTemperature.internalField(), tT)
 154     {
 155         if (totalN_sum[tT] > 0)
 156         {
 157             totalTemperature.internalField()[tT] =
 158                 totalTemperatureVTerms_sum[tT]
 159                 /(3.0 * moleculeCloud::kb * totalN_sum[tT]);
 160         }
 161         else
 162         {
 163             totalTemperature.internalField()[tT] = 0.0;
 164         }
 165     }
 166
 167     /*-----------------------------------------------------------------------*\
 168         Mean kinetic energy
 169     \*-----------------------------------------------------------------------*/
 170
 171     scalarField totalKE_sum(mesh.nCells(), 0.0);
 172
 173     forAll(allSpeciesMeanKE, mKE)
 174     {
 175         // A check for 1/0 molecules is required.
 176
 177         scalarField& singleSpeciesMeanKE
 178         (
 179             allSpeciesMeanKE[mKE].internalField()
 180         );
 181
 182         forAll(singleSpeciesMeanKE, sSMKE)
 183         {
 184             if (allSpeciesN_RU[mKE][sSMKE])
 185             {
 186                 singleSpeciesMeanKE[sSMKE] =
 187                     allSpeciesM_RU[mKE][sSMKE]
 188                    /allSpeciesN_RU[mKE][sSMKE]
 189                    /(2.0*allSpeciesN_RU[mKE][sSMKE])
 190                    *(
 191                         allSpeciesVelocityMagSquaredSum_RU[mKE][sSMKE]
 192                     );
 193
 194                 totalKE_sum[sSMKE] +=
 195                     allSpeciesM_RU[mKE][sSMKE]
 196                     /allSpeciesN_RU[mKE][sSMKE]
 197                     /2.0
 198                    *(
 199                         allSpeciesVelocityMagSquaredSum_RU[mKE][sSMKE]
 200                     );
 201             }
 202             else
 203             {
 204                 singleSpeciesMeanKE[sSMKE] = 0.0;
 205             }
 206         }
 207     }
 208
 209     forAll(totalMeanKE.internalField(), tMKE)
 210     {
 211         if (totalN_sum[tMKE] > 0)
 212         {
 213             totalMeanKE.internalField()[tMKE] =
 214                 totalKE_sum[tMKE]
 215                 /totalN_sum[tMKE];
 216         }
 217         else
 218         {
 219             totalMeanKE.internalField()[tMKE] = 0.0;
 220         }
 221     }
 222 }
 223