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

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Global
  25     meanMomentumEnergyAndNMols.H
  26
  27 Description
  28     Calculates and prints the mean momentum and energy in the system
  29     and the number of molecules.
  30
  31 \*---------------------------------------------------------------------------*/
  32
  33
  34 vector singleStepTotalLinearMomentum(vector::zero);
  35
  36 vector singleStepTotalAngularMomentum(vector::zero);
  37
  38 scalar singleStepMaxVelocityMag = 0.0;
  39
  40 scalar singleStepTotalMass = 0.0;
  41
  42 scalar singleStepTotalLinearKE = 0.0;
  43
  44 scalar singleStepTotalAngularKE = 0.0;
  45
  46 scalar singleStepTotalPE = 0.0;
  47
  48 scalar singleStepTotalrDotf = 0.0;
  49
  50 //vector singleStepCentreOfMass(vector::zero);
  51
  52 label singleStepNMols = molecules.size();
  53
  54 label singleStepDOFs = 0;
  55
  56 {
  57     forAllConstIter(IDLList<molecule>, molecules, mol)
  58     {
  59         const label molId = mol().id();
  60
  61         scalar molMass(molecules.constProps(molId).mass());
  62
  63         singleStepTotalMass += molMass;
  64
  65         //singleStepCentreOfMass += mol().position()*molMass;
  66     }
  67
  68     // if (singleStepNMols)
  69     // {
  70     //     singleStepCentreOfMass /= singleStepTotalMass;
  71     // }
  72
  73     forAllConstIter(IDLList<molecule>, molecules, mol)
  74     {
  75         const label molId = mol().id();
  76
  77         const molecule::constantProperties cP(molecules.constProps(molId));
  78
  79         scalar molMass(cP.mass());
  80
  81         const diagTensor& molMoI(cP.momentOfInertia());
  82
  83         const vector& molV(mol().v());
  84
  85         const vector& molOmega(inv(molMoI) & mol().pi());
  86
  87         vector molPiGlobal = mol().Q() & mol().pi();
  88
  89         singleStepTotalLinearMomentum += molV * molMass;
  90
  91         singleStepTotalAngularMomentum += molPiGlobal;
  92         //+((mol().position() - singleStepCentreOfMass) ^ (molV * molMass));
  93
  94         if (mag(molV) > singleStepMaxVelocityMag)
  95         {
  96             singleStepMaxVelocityMag = mag(molV);
  97         }
  98
  99         singleStepTotalLinearKE += 0.5*molMass*magSqr(molV);
 100
 101         singleStepTotalAngularKE += 0.5*(molOmega & molMoI & molOmega);
 102
 103         singleStepTotalPE += mol().potentialEnergy();
 104
 105         singleStepTotalrDotf += tr(mol().rf());
 106
 107         singleStepDOFs += cP.degreesOfFreedom();
 108     }
 109 }
 110
 111 if (Pstream::parRun())
 112 {
 113     reduce(singleStepTotalLinearMomentum, sumOp<vector>());
 114
 115     reduce(singleStepTotalAngularMomentum, sumOp<vector>());
 116
 117     reduce(singleStepMaxVelocityMag, maxOp<scalar>());
 118
 119     reduce(singleStepTotalMass, sumOp<scalar>());
 120
 121     reduce(singleStepTotalLinearKE, sumOp<scalar>());
 122
 123     reduce(singleStepTotalAngularKE, sumOp<scalar>());
 124
 125     reduce(singleStepTotalPE, sumOp<scalar>());
 126
 127     reduce(singleStepTotalrDotf, sumOp<scalar>());
 128
 129     reduce(singleStepNMols, sumOp<label>());
 130
 131     reduce(singleStepDOFs, sumOp<label>());
 132 }
 133
 134 if (singleStepNMols)
 135 {
 136     Info<< "Number of molecules in system = "
 137         << singleStepNMols << nl
 138         << "Overall number density = "
 139         << singleStepNMols/meshVolume << nl
 140         << "Overall mass density = "
 141         << singleStepTotalMass/meshVolume << nl
 142         << "Average linear momentum per molecule = "
 143         << singleStepTotalLinearMomentum/singleStepNMols << ' '
 144         << mag(singleStepTotalLinearMomentum)/singleStepNMols << nl
 145         << "Average angular momentum per molecule = "
 146         << singleStepTotalAngularMomentum << ' '
 147         << mag(singleStepTotalAngularMomentum)/singleStepNMols << nl
 148         << "Maximum |velocity| = "
 149         << singleStepMaxVelocityMag << nl
 150         << "Average linear KE per molecule = "
 151         << singleStepTotalLinearKE/singleStepNMols << nl
 152         << "Average angular KE per molecule = "
 153         << singleStepTotalAngularKE/singleStepNMols << nl
 154         << "Average PE per molecule = "
 155         << singleStepTotalPE/singleStepNMols << nl
 156         << "Average TE per molecule = "
 157         <<
 158         (
 159             singleStepTotalLinearKE
 160           + singleStepTotalAngularKE
 161           + singleStepTotalPE
 162         )
 163         /singleStepNMols
 164         << endl;
 165
 166         // Info<< singleStepNMols << " "
 167         //     << singleStepTotalMomentum/singleStepTotalMass << " "
 168         //     << singleStepMaxVelocityMag << " "
 169         //     << singleStepTotalKE/singleStepNMols << " "
 170         //     << singleStepTotalPE/singleStepNMols << " "
 171         //     << (singleStepTotalKE + singleStepTotalPE)
 172         //        /singleStepNMols << endl;
 173 }
 174 else
 175 {
 176     Info<< "No molecules in system" << endl;
 177 }
 178
 179
 180 // ************************************************************************* //