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) 2008-2009 OpenCFD Ltd.
   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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 Global
  26     meanMomentumEnergyAndNMols.H
  27
  28 Description
  29     Calculates and prints the mean momentum and energy in the system
  30     and the number of molecules.
  31
  32 \*---------------------------------------------------------------------------*/
  33
  34
  35 vector singleStepTotalLinearMomentum(vector::zero);
  36
  37 vector singleStepTotalAngularMomentum(vector::zero);
  38
  39 scalar singleStepMaxVelocityMag = 0.0;
  40
  41 scalar singleStepTotalMass = 0.0;
  42
  43 scalar singleStepTotalLinearKE = 0.0;
  44
  45 scalar singleStepTotalAngularKE = 0.0;
  46
  47 scalar singleStepTotalPE = 0.0;
  48
  49 scalar singleStepTotalrDotf = 0.0;
  50
  51 //vector singleStepCentreOfMass(vector::zero);
  52
  53 label singleStepNMols = molecules.size();
  54
  55 label singleStepDOFs = 0;
  56
  57 {
  58     IDLList<molecule>::iterator mol(molecules.begin());
  59
  60     for
  61     (
  62         mol = molecules.begin();
  63         mol != molecules.end();
  64         ++mol
  65     )
  66     {
  67         label molId = mol().id();
  68
  69         scalar molMass(molecules.constProps(molId).mass());
  70
  71         singleStepTotalMass += molMass;
  72
  73         //singleStepCentreOfMass += mol().position()*molMass;
  74     }
  75
  76     // if(singleStepNMols)
  77     // {
  78     //     singleStepCentreOfMass /= singleStepTotalMass;
  79     // }
  80
  81     for
  82     (
  83         mol = molecules.begin();
  84         mol != molecules.end();
  85         ++mol
  86     )
  87     {
  88         label molId = mol().id();
  89
  90         const molecule::constantProperties cP(molecules.constProps(molId));
  91
  92         scalar molMass(cP.mass());
  93
  94         const diagTensor& molMoI(cP.momentOfInertia());
  95
  96         const vector& molV(mol().v());
  97
  98         const vector& molOmega(inv(molMoI) & mol().pi());
  99
 100         vector molPiGlobal = mol().Q() & mol().pi();
 101
 102         singleStepTotalLinearMomentum += molV * molMass;
 103
 104         singleStepTotalAngularMomentum += molPiGlobal;
 105         //+((mol().position() - singleStepCentreOfMass) ^ (molV * molMass));
 106
 107         if(mag(molV) > singleStepMaxVelocityMag)
 108         {
 109             singleStepMaxVelocityMag = mag(molV);
 110         }
 111
 112         singleStepTotalLinearKE += 0.5*molMass*magSqr(molV);
 113
 114         singleStepTotalAngularKE += 0.5*(molOmega & molMoI & molOmega);
 115
 116         singleStepTotalPE += mol().potentialEnergy();
 117
 118         singleStepTotalrDotf += tr(mol().rf());
 119
 120         singleStepDOFs += cP.degreesOfFreedom();
 121     }
 122 }
 123
 124 if (Pstream::parRun())
 125 {
 126     reduce(singleStepTotalLinearMomentum, sumOp<vector>());
 127
 128     reduce(singleStepTotalAngularMomentum, sumOp<vector>());
 129
 130     reduce(singleStepMaxVelocityMag, maxOp<scalar>());
 131
 132     reduce(singleStepTotalMass, sumOp<scalar>());
 133
 134     reduce(singleStepTotalLinearKE, sumOp<scalar>());
 135
 136     reduce(singleStepTotalAngularKE, sumOp<scalar>());
 137
 138     reduce(singleStepTotalPE, sumOp<scalar>());
 139
 140     reduce(singleStepTotalrDotf, sumOp<scalar>());
 141
 142     reduce(singleStepNMols, sumOp<label>());
 143
 144     reduce(singleStepDOFs, sumOp<label>());
 145 }
 146
 147 if (singleStepNMols)
 148 {
 149     Info<< "Number of molecules in system = "
 150         << singleStepNMols << nl
 151         << "Overall number density = "
 152         << singleStepNMols/meshVolume << nl
 153         << "Overall mass density = "
 154         << singleStepTotalMass/meshVolume << nl
 155         << "Average linear momentum per molecule = "
 156         << singleStepTotalLinearMomentum/singleStepNMols << ' '
 157         << mag(singleStepTotalLinearMomentum)/singleStepNMols << nl
 158         << "Average angular momentum per molecule = "
 159         << singleStepTotalAngularMomentum << ' '
 160         << mag(singleStepTotalAngularMomentum)/singleStepNMols << nl
 161         << "Maximum |velocity| = "
 162         << singleStepMaxVelocityMag << nl
 163         << "Average linear KE per molecule = "
 164         << singleStepTotalLinearKE/singleStepNMols << nl
 165         << "Average angular KE per molecule = "
 166         << singleStepTotalAngularKE/singleStepNMols << nl
 167         << "Average PE per molecule = "
 168         << singleStepTotalPE/singleStepNMols << nl
 169         << "Average TE per molecule = "
 170         <<
 171         (
 172             singleStepTotalLinearKE
 173           + singleStepTotalAngularKE
 174           + singleStepTotalPE
 175         )
 176         /singleStepNMols
 177         << endl;
 178
 179         // Info << singleStepNMols << " "
 180         //     << singleStepTotalMomentum/singleStepTotalMass << " "
 181         //     << singleStepMaxVelocityMag << " "
 182         //     << singleStepTotalKE/singleStepNMols << " "
 183         //     << singleStepTotalPE/singleStepNMols << " "
 184         //     << (singleStepTotalKE + singleStepTotalPE)
 185         //        /singleStepNMols << endl;
 186 }
 187 else
 188 {
 189     Info<< "No molecules in system" << endl;
 190 }
 191
 192
 193 // ************************************************************************* //