Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

etc/config/settings.csh

@@ -41,6 +41,9 @@ alias _foamAddMan 'setenv MANPATH \!*\:${MANPATH}'
 # Set environment variables according to system type
 setenv WM_ARCH `uname -s`
 
+# Default WM_COMPILER_LIB_ARCH for 32bit
+setenv WM_COMPILER_LIB_ARCH
+
 switch ($WM_ARCH)
 case Linux:
     setenv WM_ARCH linux
@@ -269,12 +272,7 @@ case ThirdParty:
         _foamAddPath    $gccDir/bin
 
         # add compiler libraries to run-time environment
-        # 64-bit needs lib64, but 32-bit needs lib (not lib32)
+        _foamAddLib     $gccDir/lib$WM_COMPILER_LIB_ARCH
-        if ($WM_ARCH_OPTION == 64 && $?WM_COMPILER_LIB_ARCH) then
-            _foamAddLib     $gccDir/lib$WM_COMPILER_LIB_ARCH
-        else
-            _foamAddLib     $gccDir/lib
-        endif
 
         # add gmp/mpfr libraries to run-time environment
         _foamAddLib     $gmpDir/lib
@@ -371,21 +369,6 @@ unset boost_version cgal_version
 unsetenv MPI_ARCH_PATH MPI_HOME FOAM_MPI_LIBBIN
 
 switch ("$WM_MPLIB")
-case OPENMPI:
-    setenv FOAM_MPI openmpi-1.5.4
-    # optional configuration tweaks:
-    _foamSource `$WM_PROJECT_DIR/bin/foamEtcFile config/openmpi.csh`
-
-    setenv MPI_ARCH_PATH $WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER/$FOAM_MPI
-
-    # Tell OpenMPI where to find its install directory
-    setenv OPAL_PREFIX $MPI_ARCH_PATH
-
-    _foamAddPath    $MPI_ARCH_PATH/bin
-    _foamAddLib     $MPI_ARCH_PATH/lib
-    _foamAddMan     $MPI_ARCH_PATH/share/man
-    breaksw
-
 case SYSTEMOPENMPI:
     # Use the system installed openmpi, get library directory via mpicc
     setenv FOAM_MPI openmpi-system
@@ -410,13 +393,36 @@ case SYSTEMOPENMPI:
     unset libDir
     breaksw
 
+case OPENMPI:
+    setenv FOAM_MPI openmpi-1.5.4
+    # optional configuration tweaks:
+    _foamSource `$WM_PROJECT_DIR/bin/foamEtcFile config/openmpi.csh`
+
+    setenv MPI_ARCH_PATH $WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER/$FOAM_MPI
+
+    # Tell OpenMPI where to find its install directory
+    setenv OPAL_PREFIX $MPI_ARCH_PATH
+
+    _foamAddPath    $MPI_ARCH_PATH/bin
+
+    # 64-bit on OpenSuSE 12.1 uses lib64 others use lib
+    _foamAddLib     $MPI_ARCH_PATH/lib$WM_COMPILER_LIB_ARCH
+    _foamAddLib     $MPI_ARCH_PATH/lib
+
+    _foamAddMan     $MPI_ARCH_PATH/share/man
+    breaksw
+
 case MPICH:
     setenv FOAM_MPI mpich2-1.1.1p1
     setenv MPI_HOME $WM_THIRD_PARTY_DIR/$FOAM_MPI
     setenv MPI_ARCH_PATH $WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER/$FOAM_MPI
 
     _foamAddPath    $MPI_ARCH_PATH/bin
+
+    # 64-bit on OpenSuSE 12.1 uses lib64 others use lib
+    _foamAddLib     $MPI_ARCH_PATH/lib$WM_COMPILER_LIB_ARCH
     _foamAddLib     $MPI_ARCH_PATH/lib
+
     _foamAddMan     $MPI_ARCH_PATH/share/man
     breaksw
 
@@ -427,7 +433,11 @@ case MPICH-GM:
     setenv GM_LIB_PATH /opt/gm/lib64
 
     _foamAddPath    $MPI_ARCH_PATH/bin
+
+    # 64-bit on OpenSuSE 12.1 uses lib64 others use lib
+    _foamAddLib     $MPI_ARCH_PATH/lib$WM_COMPILER_LIB_ARCH
     _foamAddLib     $MPI_ARCH_PATH/lib
+
     _foamAddLib     $GM_LIB_PATH
     breaksw
 

etc/config/settings.sh

@@ -294,14 +294,7 @@ OpenFOAM | ThirdParty)
         _foamAddPath    $gccDir/bin
 
         # add compiler libraries to run-time environment
-        # 64-bit needs lib64, but 32-bit needs lib (not lib32)
+        _foamAddLib     $gccDir/lib$WM_COMPILER_LIB_ARCH
-        if [ "$WM_ARCH_OPTION" = 64 ]
-        then
-            _foamAddLib     $gccDir/lib$WM_COMPILER_LIB_ARCH
-        else
-            _foamAddLib     $gccDir/lib
-        fi
-
 
         # add gmp/mpfr libraries to run-time environment
         _foamAddLib     $gmpDir/lib
@@ -401,21 +394,6 @@ unset boost_version cgal_version
 unset MPI_ARCH_PATH MPI_HOME FOAM_MPI_LIBBIN
 
 case "$WM_MPLIB" in
-OPENMPI)
-    export FOAM_MPI=openmpi-1.5.4
-    # optional configuration tweaks:
-    _foamSource `$WM_PROJECT_DIR/bin/foamEtcFile config/openmpi.sh`
-
-    export MPI_ARCH_PATH=$WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER/$FOAM_MPI
-
-    # Tell OpenMPI where to find its install directory
-    export OPAL_PREFIX=$MPI_ARCH_PATH
-
-    _foamAddPath    $MPI_ARCH_PATH/bin
-    _foamAddLib     $MPI_ARCH_PATH/lib
-    _foamAddMan     $MPI_ARCH_PATH/share/man
-    ;;
-
 SYSTEMOPENMPI)
     # Use the system installed openmpi, get library directory via mpicc
     export FOAM_MPI=openmpi-system
@@ -441,13 +419,36 @@ SYSTEMOPENMPI)
     unset libDir
     ;;
 
+OPENMPI)
+    export FOAM_MPI=openmpi-1.5.4
+    # optional configuration tweaks:
+    _foamSource `$WM_PROJECT_DIR/bin/foamEtcFile config/openmpi.sh`
+
+    export MPI_ARCH_PATH=$WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER/$FOAM_MPI
+
+    # Tell OpenMPI where to find its install directory
+    export OPAL_PREFIX=$MPI_ARCH_PATH
+
+    _foamAddPath    $MPI_ARCH_PATH/bin
+
+    # 64-bit on OpenSuSE 12.1 uses lib64 others use lib
+    _foamAddLib     $MPI_ARCH_PATH/lib$WM_COMPILER_LIB_ARCH
+    _foamAddLib     $MPI_ARCH_PATH/lib
+
+    _foamAddMan     $MPI_ARCH_PATH/share/man
+    ;;
+
 MPICH)
     export FOAM_MPI=mpich2-1.1.1p1
     export MPI_HOME=$WM_THIRD_PARTY_DIR/$FOAM_MPI
     export MPI_ARCH_PATH=$WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER/$FOAM_MPI
 
     _foamAddPath    $MPI_ARCH_PATH/bin
+
+    # 64-bit on OpenSuSE 12.1 uses lib64 others use lib
+    _foamAddLib     $MPI_ARCH_PATH/lib$WM_COMPILER_LIB_ARCH
     _foamAddLib     $MPI_ARCH_PATH/lib
+
     _foamAddMan     $MPI_ARCH_PATH/share/man
     ;;
 
@@ -458,7 +459,11 @@ MPICH-GM)
     export GM_LIB_PATH=/opt/gm/lib64
 
     _foamAddPath    $MPI_ARCH_PATH/bin
+
+    # 64-bit on OpenSuSE 12.1 uses lib64 others use lib
+    _foamAddLib     $MPI_ARCH_PATH/lib$WM_COMPILER_LIB_ARCH
     _foamAddLib     $MPI_ARCH_PATH/lib
+
     _foamAddLib     $GM_LIB_PATH
     ;;
 

src/OpenFOAM/db/IOstreams/StringStreams/OStringStream.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -116,14 +116,14 @@ public:
 
         // Edit
 
-            //- Clear the OStringStream
+            //- Rewind the OStringStream
             void rewind()
             {
-#               if __GNUC__ < 4 && __GNUC_MINOR__ < 4
+                #if __GNUC__ < 4 && __GNUC_MINOR__ < 4
                 stdStream().rdbuf()->pubsetbuf(" ", 1);
-#               else
+                #else
                 stdStream().rdbuf()->pubseekpos(0);
-#               endif
+                #endif
             }
 
 

src/fieldSources/Make/files

@@ -15,6 +15,10 @@ basicSource/rotorDiskSource/profileModel/profileModel.C
 basicSource/rotorDiskSource/profileModel/profileModelList.C
 basicSource/rotorDiskSource/profileModel/lookup/lookupProfile.C
 basicSource/rotorDiskSource/profileModel/series/seriesProfile.C
+basicSource/rotorDiskSource/trimModel/trimModel/trimModel.C
+basicSource/rotorDiskSource/trimModel/trimModel/trimModelNew.C
+basicSource/rotorDiskSource/trimModel/fixed/fixedTrim.C
+basicSource/rotorDiskSource/trimModel/targetForce/targetForceTrim.C
 
 basicSource/actuationDiskSource/actuationDiskSource.C
 basicSource/radialActuationDiskSource/radialActuationDiskSource.C
@@ -25,4 +29,4 @@ $(interRegion)/constantHeatTransfer/constantHeatTransfer.C
 $(interRegion)/tabulatedHeatTransfer/tabulatedHeatTransfer.C
 $(interRegion)/variableHeatTransfer/variableHeatTransfer.C
 
-LIB = $(FOAM_LIBBIN)/libfieldSources
+LIB = $(FOAM_LIBBIN)/libfieldSources

src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.C

@@ -26,8 +26,8 @@ License
 #include "rotorDiskSource.H"
 #include "addToRunTimeSelectionTable.H"
 #include "mathematicalConstants.H"
+#include "trimModel.H"
 #include "unitConversion.H"
-#include "geometricOneField.H"
 #include "fvMatrices.H"
 #include "syncTools.H"
 
@@ -123,6 +123,8 @@ void Foam::rotorDiskSource::checkData()
 
 void Foam::rotorDiskSource::setFaceArea(vector& axis, const bool correct)
 {
+    area_ = 0.0;
+
     static const scalar tol = 0.8;
 
     const label nInternalFaces = mesh_.nInternalFaces();
@@ -299,12 +301,15 @@ void Foam::rotorDiskSource::createCoordinateSystem()
             reduce(axis, maxMagSqrOp<vector>());
             axis /= mag(axis);
 
-            // axis direction is somewhat arbitrary - check if user needs
+            // correct the axis direction using a point above the rotor
-            // needs to reverse
-            bool reverse(readBool(coeffs_.lookup("reverseAxis")));
-            if (reverse)
             {
-                axis *= -1.0;
+                vector pointAbove(coeffs_.lookup("pointAbove"));
+                vector dir = pointAbove - origin;
+                dir /= mag(dir);
+                if ((dir & axis) < 0)
+                {
+                    axis *= -1.0;
+                }
             }
 
             coeffs_.lookup("refDirection") >> refDir;
@@ -380,9 +385,6 @@ void Foam::rotorDiskSource::constructGeometry()
         scalar cPos = cos(beta);
         scalar sPos = sin(beta);
         invR_[i] = tensor(cPos, 0.0, -sPos, 0.0, 1.0, 0.0, sPos, 0.0, cPos);
-
-        // geometric angle of attack - not including twist [radians]
-        alphag_[i] = trim_.alphaC + trim_.A*cos(psi) + trim_.B*sin(psi);
     }
 }
 
@@ -444,13 +446,12 @@ Foam::rotorDiskSource::rotorDiskSource
     inletVelocity_(vector::zero),
     tipEffect_(1.0),
     flap_(),
-    trim_(),
+    trim_(trimModel::New(*this, coeffs_)),
     blade_(coeffs_.subDict("blade")),
     profiles_(coeffs_.subDict("profiles")),
     x_(cells_.size(), vector::zero),
     R_(cells_.size(), I),
     invR_(cells_.size(), I),
-    alphag_(cells_.size(), 0.0),
     area_(cells_.size(), 0.0),
     coordSys_(false),
     rMax_(0.0)
@@ -467,34 +468,179 @@ Foam::rotorDiskSource::~rotorDiskSource()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+void Foam::rotorDiskSource::calculate
+(
+    const vectorField& U,
+    const scalarField& alphag,
+    vectorField& force,
+    const bool divideVolume,
+    const bool output
+) const
+{
+    tmp<volScalarField> trho;
+    if (rhoName_ != "none")
+    {
+        trho = mesh_.lookupObject<volScalarField>(rhoName_);
+    }
+
+    const scalarField& V = mesh_.V();
+
+
+    // logging info
+    scalar dragEff = 0.0;
+    scalar liftEff = 0.0;
+    scalar AOAmin = GREAT;
+    scalar AOAmax = -GREAT;
+
+    forAll(cells_, i)
+    {
+        if (area_[i] > ROOTVSMALL)
+        {
+            const label cellI = cells_[i];
+
+            const scalar radius = x_[i].x();
+
+            // velocity in local cylindrical reference frame
+            vector Uc = coordSys_.localVector(U[cellI]);
+
+            // apply correction in local system due to coning
+            Uc = R_[i] & Uc;
+
+            // set radial component of velocity to zero
+            Uc.x() = 0.0;
+
+            // remove blade linear velocity from blade normal component
+            Uc.y() -= radius*omega_;
+
+            // determine blade data for this radius
+            // i1 = index of upper bound data point in blade list
+            scalar twist = 0.0;
+            scalar chord = 0.0;
+            label i1 = -1;
+            label i2 = -1;
+            scalar invDr = 0.0;
+            blade_.interpolate(radius, twist, chord, i1, i2, invDr);
+
+            // effective angle of attack
+            scalar alphaEff =
+                mathematical::pi + atan2(Uc.z(), Uc.y()) - (alphag[i] + twist);
+
+            if (alphaEff > mathematical::pi)
+            {
+                alphaEff -= mathematical::twoPi;
+            }
+            if (alphaEff < -mathematical::pi)
+            {
+                alphaEff += mathematical::twoPi;
+            }
+
+            AOAmin = min(AOAmin, alphaEff);
+            AOAmax = max(AOAmax, alphaEff);
+
+            // determine profile data for this radius and angle of attack
+            const label profile1 = blade_.profileID()[i1];
+            const label profile2 = blade_.profileID()[i2];
+
+            scalar Cd1 = 0.0;
+            scalar Cl1 = 0.0;
+            profiles_[profile1].Cdl(alphaEff, Cd1, Cl1);
+
+            scalar Cd2 = 0.0;
+            scalar Cl2 = 0.0;
+            profiles_[profile2].Cdl(alphaEff, Cd2, Cl2);
+
+            scalar Cd = invDr*(Cd2 - Cd1) + Cd1;
+            scalar Cl = invDr*(Cl2 - Cl1) + Cl1;
+
+            // apply tip effect for blade lift
+            scalar tipFactor = neg(radius/rMax_ - tipEffect_);
+
+            // calculate forces perpendicular to blade
+            scalar pDyn = 0.5*magSqr(Uc);
+            if (trho.valid())
+            {
+                pDyn *= trho()[cellI];
+            }
+
+            scalar f = pDyn*chord*nBlades_*area_[i]/mathematical::twoPi;
+            vector localForce = vector(0.0, f*Cd, tipFactor*f*Cl);
+
+            // convert force from local coning system into rotor cylindrical
+            localForce = invR_[i] & localForce;
+
+            // accumulate forces
+            dragEff += localForce.y();
+            liftEff += localForce.z();
+
+            // convert force to global cartesian co-ordinate system
+            force[cellI] = coordSys_.globalVector(localForce);
+
+            if (divideVolume)
+            {
+                force[cellI] /= V[cellI];
+            }
+        }
+    }
+
+
+    if (output)
+    {
+        reduce(AOAmin, minOp<scalar>());
+        reduce(AOAmax, maxOp<scalar>());
+        reduce(dragEff, sumOp<scalar>());
+        reduce(liftEff, sumOp<scalar>());
+
+        Info<< type() << " output:" << nl
+            << "    min/max(AOA)   = " << radToDeg(AOAmin) << ", "
+            << radToDeg(AOAmax) << nl
+            << "    Effective drag = " << dragEff << nl
+            << "    Effective lift = " << liftEff << endl;
+    }
+}
+
+
 void Foam::rotorDiskSource::addSup(fvMatrix<vector>& eqn, const label fieldI)
 {
-    // add source to rhs of eqn
+    dimensionSet dims = dimless;
-
-    const volVectorField& U = eqn.psi();
-
     if (eqn.dimensions() == dimForce)
     {
         coeffs_.lookup("rhoName") >> rhoName_;
-
+        dims.reset(dimForce/dimVolume);
-        const volScalarField& rho =
-            mesh_.lookupObject<volScalarField>(rhoName_);
-
-        eqn -= calculateForces
-            (
-                rho.internalField(),
-                inflowVelocity(U),
-                dimForce/dimVolume
-            );
     }
     else
     {
-        eqn -= calculateForces
+        dims.reset(dimForce/dimVolume/dimDensity);
-            (
+    }
-                oneField(),
+
-                inflowVelocity(U),
+    volVectorField force
-                dimForce/dimVolume/dimDensity
+    (
-            );
+        IOobject
+        (
+            "rotorForce",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedVector("zero", dims, vector::zero)
+    );
+
+    const volVectorField& U = eqn.psi();
+
+    const vectorField Uin = inflowVelocity(U);
+
+    trim_->correct(Uin, force);
+
+    calculate(Uin, trim_->thetag(), force);
+
+
+    // add source to rhs of eqn
+    eqn -= force;
+
+    if (mesh_.time().outputTime())
+    {
+        force.write();
     }
 }
 
@@ -527,12 +673,10 @@ bool Foam::rotorDiskSource::read(const dictionary& dict)
         flapCoeffs.lookup("beta1") >> flap_.beta1;
         flapCoeffs.lookup("beta2") >> flap_.beta2;
         flap_.beta0 = degToRad(flap_.beta0);
+        flap_.beta1 = degToRad(flap_.beta1);
+        flap_.beta2 = degToRad(flap_.beta2);
 
-        const dictionary& trimCoeffs(coeffs_.subDict("trimCoeffs"));
+        trim_->read(coeffs_);        
-        trimCoeffs.lookup("alphaC") >> trim_.alphaC;
-        trimCoeffs.lookup("A") >> trim_.A;
-        trimCoeffs.lookup("B") >> trim_.B;
-        trim_.alphaC = degToRad(trim_.alphaC);
 
         checkData();
 
@@ -542,7 +686,7 @@ bool Foam::rotorDiskSource::read(const dictionary& dict)
 
         if (debug)
         {
-            writeField("alphag", alphag_, true);
+            writeField("alphag", trim_->thetag()(), true);
             writeField("faceArea", area_, true);
         }
 

src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,13 +37,16 @@ Description
             fieldNames      (U);    // names of fields on which to apply source
             rhoName         rho;    // density field if compressible case
             nBlades         3;      // number of blades
-            tip effect      0.96;   // normalised radius above which lift = 0
+            tipEffect       0.96;   // normalised radius above which lift = 0
 
             inletFlowType   local;  // inlet flow type specification
 
             geometryMode    auto;   // geometry specification
 
             refDirection    (-1 0 0); // reference direction
+                                      // - used as reference for psi angle
+
+            trimModel       fixed;  // fixed || targetForce
 
             flapCoeffs
             {
@@ -51,16 +54,12 @@ Description
                 beta1           0;  // lateral flapping coeff
                 beta2           0;  // longitudinal flapping coeff
             }
-            trimCoeffs
+
-            {
-                alphac          15; // collective pitch angle [deg]
-                A               0;  // lateral cyclic coeff
-                B               0;  // longitudinal cyclic coeff
-            }
             blade
             {
                 ...
             }
+
             profiles
             {
                 ...
@@ -102,6 +101,9 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
+class trimModel;
+
 /*---------------------------------------------------------------------------*\
                       Class rotorDiskSource Declaration
 \*---------------------------------------------------------------------------*/
@@ -140,13 +142,6 @@ protected:
             scalar beta2;   // longitudinal flapping coeff
         };
 
-        struct trimData
-        {
-            scalar alphaC;  // collective pitch angle
-            scalar A;       // lateral cyclic coeff
-            scalar B;       // longitudinal cyclic coeff
-        };
-
 
     // Protected data
 
@@ -172,8 +167,8 @@ protected:
         //- Blade flap coefficients [rad/s]
         flapData flap_;
 
-        //- Blad trim coefficients
+        //- Trim model
-        trimData trim_;
+        autoPtr<trimModel> trim_;
 
         //- Blade data
         bladeModel blade_;
@@ -181,7 +176,8 @@ protected:
         //- Profile data
         profileModelList profiles_;
 
-        //- Cell centre positions in local rotor frame (Cartesian x, y, z)
+        //- Cell centre positions in local rotor frame
+        //  (Cylindrical r, theta, z)
         List<point> x_;
 
         //- Rotation tensor for flap angle
@@ -190,9 +186,6 @@ protected:
         //- Inverse rotation tensor for flap angle
         List<tensor> invR_;
 
-        //- Geometric angle of attack [deg]
-        List<scalar> alphag_;
-
         //- Area [m2]
         List<scalar> area_;
 
@@ -220,15 +213,6 @@ protected:
         //- Return the inlet flow field
         tmp<vectorField> inflowVelocity(const volVectorField& U) const;
 
-        //- Calculate forces
-        template<class RhoType>
-        tmp<volVectorField> calculateForces
-        (
-            const RhoType& rho,
-            const vectorField& U,
-            const dimensionSet& dims
-        );
-
         //- Helper function to write rotor values
         template<class Type>
         void writeField
@@ -264,6 +248,29 @@ public:
 
     // Member Functions
 
+        // Access
+            
+            //- Return the cell centre positions in local rotor frame
+            // (Cylindrical r, theta, z)
+            inline const List<point>& x() const;
+
+            //- Return the rotor co-ordinate system (r, theta, z)
+            inline const cylindricalCS& coordSys() const;
+
+
+        // Evaluation
+
+            //- Calculate forces
+            void calculate
+            (
+                const vectorField& U,
+                const scalarField& alphag,
+                vectorField& force,
+                const bool divideVolume = true,
+                const bool output = true
+            ) const;
+
+
         // Source term addition
 
             //- Source term to fvMatrix<vector>
@@ -286,6 +293,10 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "rotorDiskSourceI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #ifdef NoRepository
     #include "rotorDiskSourceTemplates.C"
 #endif

src/fieldSources/basicSource/rotorDiskSource/rotorDiskSourceI.H

@@ -0,0 +1,41 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "rotorDiskSource.H"
+
+const Foam::List<Foam::point>& Foam::rotorDiskSource::x() const
+{
+    return x_;
+}
+
+
+const Foam::cylindricalCS& Foam::rotorDiskSource::coordSys() const
+{
+    return coordSys_;
+}
+
+
+// ************************************************************************* //
+

src/fieldSources/basicSource/rotorDiskSource/rotorDiskSourceTemplates.C

@@ -24,13 +24,8 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "rotorDiskSource.H"
-#include "addToRunTimeSelectionTable.H"
-#include "mathematicalConstants.H"
-#include "unitConversion.H"
 #include "volFields.H"
 
-using namespace Foam::constant::mathematical;
-
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 template<class Type>
@@ -81,149 +76,4 @@ void Foam::rotorDiskSource::writeField
 }
 
 
-template<class RhoType>
-Foam::tmp<Foam::volVectorField> Foam::rotorDiskSource::calculateForces
-(
-    const RhoType& rho,
-    const vectorField& U,
-    const dimensionSet& dims
-)
-{
-    tmp<volVectorField> tForce
-    (
-        new volVectorField
-        (
-            IOobject
-            (
-                "rotorForce",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh_,
-            dimensionedVector("zero", dims, vector::zero)
-        )
-    );
-
-    vectorField& force = tForce().internalField();
-    const scalarField& V = mesh_.V();
-
-
-    // logging info
-    scalar dragEff = 0.0;
-    scalar liftEff = 0.0;
-    scalar AOAmin = GREAT;
-    scalar AOAmax = -GREAT;
-
-    forAll(cells_, i)
-    {
-        if (area_[i] > ROOTVSMALL)
-        {
-            const label cellI = cells_[i];
-
-            const scalar radius = x_[i].x();
-
-            // velocity in local cylindrical reference frame
-            vector Uc = coordSys_.localVector(U[cellI]);
-
-            // apply correction in local system due to coning
-            Uc = R_[i] & Uc;
-
-            // set radial component of velocity to zero
-            Uc.x() = 0.0;
-
-            // remove blade linear velocity from blade normal component
-            Uc.y() -= radius*omega_;
-
-            // velocity magnitude
-            scalar magUc = mag(Uc);
-
-            // determine blade data for this radius
-            // i1 = index of upper bound data point in blade list
-            scalar twist = 0.0;
-            scalar chord = 0.0;
-            label i1 = -1;
-            label i2 = -1;
-            scalar invDr = 0.0;
-            blade_.interpolate(radius, twist, chord, i1, i2, invDr);
-
-            // effective angle of attack
-            scalar alphaEff = pi + atan2(Uc.z(), Uc.y()) - (alphag_[i] + twist);
-            if (alphaEff > pi)
-            {
-                alphaEff -= twoPi;
-            }
-            if (alphaEff < -pi)
-            {
-                alphaEff += twoPi;
-            }
-
-            AOAmin = min(AOAmin, alphaEff);
-            AOAmax = max(AOAmax, alphaEff);
-
-            // determine profile data for this radius and angle of attack
-            const label profile1 = blade_.profileID()[i1];
-            const label profile2 = blade_.profileID()[i2];
-
-            scalar Cd1 = 0.0;
-            scalar Cl1 = 0.0;
-            profiles_[profile1].Cdl(alphaEff, Cd1, Cl1);
-
-            scalar Cd2 = 0.0;
-            scalar Cl2 = 0.0;
-            profiles_[profile2].Cdl(alphaEff, Cd2, Cl2);
-
-            scalar Cd = invDr*(Cd2 - Cd1) + Cd1;
-            scalar Cl = invDr*(Cl2 - Cl1) + Cl1;
-
-            // apply tip effect for blade lift
-            scalar tipFactor = 1.0;
-            if (radius/rMax_ > tipEffect_)
-            {
-                tipFactor = 0.0;
-            }
-
-            // calculate forces perpendicular to blade
-            scalar pDyn = 0.5*rho[cellI]*sqr(magUc);
-            scalar f = pDyn*chord*nBlades_*area_[i]/twoPi;
-            vector localForce = vector(0.0, f*Cd, tipFactor*f*Cl);
-
-            // convert force from local coning system into rotor cylindrical
-            localForce = invR_[i] & localForce;
-
-            // accumulate forces
-            dragEff += localForce.y();
-            liftEff += localForce.z();
-
-            // convert force to global cartesian co-ordinate system
-            force[cellI] = coordSys_.globalVector(localForce);
-
-            force[cellI] /= V[cellI];
-        }
-    }
-
-
-    if (mesh_.time().outputTime())
-    {
-        tForce().write();
-    }
-
-
-    reduce(AOAmin, minOp<scalar>());
-    reduce(AOAmax, maxOp<scalar>());
-    reduce(dragEff, sumOp<scalar>());
-    reduce(liftEff, sumOp<scalar>());
-
-    Info<< type() << " output:" << nl
-        << "    min/max(AOA)   = " << radToDeg(AOAmin) << ", "
-        << radToDeg(AOAmax) << nl
-        << "    Effective drag = " << dragEff << nl
-        << "    Effective lift = " << liftEff << endl;
-
-
-    return tForce;
-}
-
-
 // ************************************************************************* //

src/fieldSources/basicSource/rotorDiskSource/trimModel/fixed/fixedTrim.C

@@ -0,0 +1,96 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fixedTrim.H"
+#include "addToRunTimeSelectionTable.H"
+#include "unitConversion.H"
+#include "mathematicalConstants.H"
+
+using namespace Foam::constant;
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(fixedTrim, 0);
+
+    addToRunTimeSelectionTable(trimModel, fixedTrim, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fixedTrim::fixedTrim(const rotorDiskSource& rotor, const dictionary& dict)
+:
+    trimModel(rotor, dict, typeName),
+    thetag_(rotor.cells().size(), 0.0)
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * //
+
+Foam::fixedTrim::~fixedTrim()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::fixedTrim::read(const dictionary& dict)
+{
+    trimModel::read(dict);
+
+    scalar theta0 = degToRad(readScalar(coeffs_.lookup("theta0")));
+    scalar theta1c = degToRad(readScalar(coeffs_.lookup("theta1c")));
+    scalar theta1s = degToRad(readScalar(coeffs_.lookup("theta1s")));
+
+    const List<vector>& x = rotor_.x();
+    forAll(thetag_, i)
+    {
+        scalar psi = x[i].y();
+        if (psi < 0)
+        {
+            psi += mathematical::twoPi;
+        }
+
+        thetag_[i] = theta0 + theta1c*cos(psi) + theta1s*sin(psi);
+    }
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::fixedTrim::thetag() const
+{
+    return tmp<scalarField>(thetag_);
+}
+
+
+void Foam::fixedTrim::correct(const vectorField& U, vectorField& force)
+{
+    // do nothing
+}
+
+
+// ************************************************************************* //

src/fieldSources/basicSource/rotorDiskSource/trimModel/fixed/fixedTrim.H

@@ -0,0 +1,95 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::fixedTrim
+
+Description
+    Fixed trim coefficients
+
+SourceFiles
+    fixedTrim.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fixedTrim_H
+#define fixedTrim_H
+
+#include "trimModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class fixedTrim Declaration
+\*---------------------------------------------------------------------------*/
+
+class fixedTrim
+:
+    public trimModel
+{
+
+protected:
+
+    // Protected data
+
+        //- Geometric angle of attack [rad]
+        scalarField thetag_;
+
+
+public:
+
+    //- Run-time type information
+    TypeName("fixedTrim");
+
+    //- Constructor
+    fixedTrim(const rotorDiskSource& rotor, const dictionary& dict);
+
+    //- Destructor
+    virtual ~fixedTrim();
+
+
+    // Member functions
+
+        //- Read
+        void read(const dictionary& dict);
+
+        //- Return the geometric angle of attack [rad]
+        virtual tmp<scalarField> thetag() const;
+
+        //- Correct the model
+        virtual void correct(const vectorField& U, vectorField& force);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/fieldSources/basicSource/rotorDiskSource/trimModel/targetForce/targetForceTrim.C

@@ -0,0 +1,234 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "targetForceTrim.H"
+#include "addToRunTimeSelectionTable.H"
+#include "unitConversion.H"
+#include "mathematicalConstants.H"
+
+using namespace Foam::constant;
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(targetForceTrim, 0);
+
+    addToRunTimeSelectionTable(trimModel, targetForceTrim, dictionary);
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+Foam::vector Foam::targetForceTrim::calcForce
+(
+    const vectorField& U,
+    const scalarField& thetag,
+    vectorField& force    
+) const
+{
+    rotor_.calculate(U, thetag, force, false, false);
+
+    const labelList& cells = rotor_.cells();
+    const vectorField& C = rotor_.mesh().C();
+
+    const vector& origin = rotor_.coordSys().origin();
+    const vector& rollAxis = rotor_.coordSys().e1();
+    const vector& pitchAxis = rotor_.coordSys().e2();
+    const vector& yawAxis = rotor_.coordSys().e3();
+
+    vector f(vector::zero);
+    forAll(cells, i)
+    {
+        label cellI = cells[i];
+
+        vector moment = force[cellI]^(C[cellI] - origin);
+        f[0] += force[cellI] & yawAxis;
+        f[1] += moment & pitchAxis;
+        f[2] += moment & rollAxis;
+    }
+
+    reduce(f, sumOp<vector>());
+    
+    return f;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::targetForceTrim::targetForceTrim
+(
+    const rotorDiskSource& rotor,
+    const dictionary& dict
+)
+:
+    trimModel(rotor, dict, typeName),
+    calcFrequency_(-1),
+    target_(vector::zero),
+    theta_(vector::zero),
+    nIter_(50),
+    tol_(1e-8),
+    relax_(1.0),
+    dTheta_(degToRad(0.1))
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * //
+
+Foam::targetForceTrim::~targetForceTrim()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::targetForceTrim::read(const dictionary& dict)
+{
+    trimModel::read(dict);
+
+    const dictionary& targetDict(coeffs_.subDict("target"));
+    target_[0] = readScalar(targetDict.lookup("fThrust"));
+    target_[1] = readScalar(targetDict.lookup("mPitch"));
+    target_[2] = readScalar(targetDict.lookup("mRoll"));
+
+    const dictionary& pitchAngleDict(coeffs_.subDict("pitchAngles"));
+    theta_[0] = degToRad(readScalar(pitchAngleDict.lookup("theta0Ini")));
+    theta_[1] = degToRad(readScalar(pitchAngleDict.lookup("theta1cIni")));
+    theta_[2] = degToRad(readScalar(pitchAngleDict.lookup("theta1sIni")));
+
+    coeffs_.lookup("calcFrequency") >> calcFrequency_;
+
+    coeffs_.readIfPresent("nIter", nIter_);
+    coeffs_.readIfPresent("tol", tol_);
+    coeffs_.readIfPresent("relax", relax_);
+
+    if (coeffs_.readIfPresent("dTheta", dTheta_))
+    {
+        dTheta_ = degToRad(dTheta_);
+    }
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::targetForceTrim::thetag() const
+{
+    const List<vector>& x = rotor_.x();
+
+    tmp<scalarField> ttheta(new scalarField(x.size()));
+    scalarField& t = ttheta();
+
+    forAll(t, i)
+    {
+        scalar psi = x[i].y();
+        if (psi < 0)
+        {
+            psi += mathematical::twoPi;
+        }
+
+        t[i] = theta_[0] + theta_[1]*cos(psi) + theta_[2]*sin(psi);
+    }
+
+    return ttheta;
+}
+
+
+void Foam::targetForceTrim::correct(const vectorField& U, vectorField& force)
+{
+    if (rotor_.mesh().time().timeIndex() % calcFrequency_ == 0)
+    {
+        // iterate to find new pitch angles to achieve target force
+        scalar err = GREAT;
+        label iter = 0;
+        tensor J(tensor::zero);
+
+        while ((err > tol_) && (iter < nIter_))
+        {
+            // cache initial theta vector
+            vector theta0(theta_);
+
+            // set initial values
+            vector old = calcForce(U, thetag(), force);
+
+            // construct Jacobian by perturbing the pitch angles
+            // by +/-(dTheta_/2)
+            for (label pitchI = 0; pitchI < 3; pitchI++)
+            {
+                theta_[pitchI] -= dTheta_/2.0;
+                vector f0 = calcForce(U, thetag(), force);
+
+                theta_[pitchI] += dTheta_;
+                vector f1 = calcForce(U, thetag(), force);
+
+                vector ddTheta = (f1 - f0)/dTheta_;
+
+                J[pitchI + 0] = ddTheta[0];
+                J[pitchI + 3] = ddTheta[1];
+                J[pitchI + 6] = ddTheta[2];
+
+                theta_ = theta0;
+            }
+
+            // calculate the change in pitch angle vector
+            vector dt = inv(J) & (target_ - old);
+
+            // update pitch angles
+            vector thetaNew = theta_ + relax_*dt;
+
+            // update error
+            err = mag(thetaNew - theta_);
+
+            // update for next iteration
+            theta_ = thetaNew;
+            iter++;
+        }
+
+        if (iter == nIter_)
+        {
+            WarningIn
+            (
+                "void Foam::targetForceTrim::correct"
+                "("
+                    "const vectorField&, "
+                    "vectorField&"
+                ")"
+            )   << "Trim routine not converged in " << iter
+                << " iterations, max residual = " << err << endl;
+        }
+        else
+        {
+            Info<< type() << ": converged in " << iter
+                << " iterations" << endl;
+        }
+
+        Info<< "    new pitch angles:" << nl
+            << "        theta0  = " << radToDeg(theta_[0]) << nl
+            << "        theta1c = " << radToDeg(theta_[1]) << nl
+            << "        theta1s = " << radToDeg(theta_[2]) << nl
+            << endl;
+    }
+}
+
+
+// ************************************************************************* //

src/fieldSources/basicSource/rotorDiskSource/trimModel/targetForce/targetForceTrim.H

@@ -0,0 +1,126 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::targetForceTrim
+
+Description
+    Target force trim coefficients
+
+SourceFiles
+    targetForceTrim.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef targetForceTrim_H
+#define targetForceTrim_H
+
+#include "trimModel.H"
+#include "tensor.H"
+#include "vector.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class targetForceTrim Declaration
+\*---------------------------------------------------------------------------*/
+
+class targetForceTrim
+:
+    public trimModel
+{
+
+protected:
+
+    // Protected data
+
+        //- Number of iterations between calls to 'correct'
+        label calcFrequency_;
+
+        //- Target force [N]
+        vector target_;
+
+        //- Pitch angles (collective, roll, pitch) [rad]
+        vector theta_;
+
+        //- Maximum number of iterations in trim routine
+        label nIter_;
+
+        //- Convergence tolerance
+        scalar tol_;
+
+        //- Under-relaxation coefficient
+        scalar relax_;
+
+        //- Perturbation angle used to determine jacobian
+        scalar dTheta_;
+
+
+    // Protected member functions
+
+        //- Calculate the rotor forces
+        vector calcForce
+        (
+            const vectorField& U,
+            const scalarField& alphag,
+            vectorField& force    
+        ) const;
+
+
+public:
+
+    //- Run-time type information
+    TypeName("targetForceTrim");
+
+    //- Constructor
+    targetForceTrim(const rotorDiskSource& rotor, const dictionary& dict);
+
+    //- Destructor
+    virtual ~targetForceTrim();
+
+
+    // Member functions
+
+        //- Read
+        void read(const dictionary& dict);
+
+        //- Return the geometric angle of attack [rad]
+        virtual tmp<scalarField> thetag() const;
+
+        //- Correct the model
+        virtual void correct(const vectorField& U, vectorField& force);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/fieldSources/basicSource/rotorDiskSource/trimModel/trimModel/trimModel.C

@@ -0,0 +1,67 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "trimModel.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(trimModel, 0);
+    defineRunTimeSelectionTable(trimModel, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::trimModel::trimModel
+(
+    const rotorDiskSource& rotor,
+    const dictionary& dict,
+    const word& name
+)
+:
+    rotor_(rotor),
+    name_(name),
+    coeffs_(dictionary::null)
+{
+    read(dict);
+}
+
+// * * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * //
+
+Foam::trimModel::~trimModel()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::trimModel::read(const dictionary& dict)
+{
+    coeffs_ = dict.subDict(name_ + "Coeffs");
+}
+
+
+// ************************************************************************* //

src/fieldSources/basicSource/rotorDiskSource/trimModel/trimModel/trimModel.H

@@ -0,0 +1,135 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::trimModel
+
+Description
+    Trim model base class
+
+SourceFiles
+    trimModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef trimModel_H
+#define trimModel_H
+
+#include "rotorDiskSource.H"
+#include "dictionary.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class trimModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class trimModel
+{
+
+protected:
+
+    // Protected data
+
+        //- Reference to the rotor source model
+        const rotorDiskSource& rotor_;
+
+        //- Name of model
+        const word name_;
+
+        //- Coefficients dictionary
+        dictionary coeffs_;
+
+
+public:
+
+    //- Run-time type information
+    TypeName("trimModel");
+
+
+    // Declare runtime constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            trimModel,
+            dictionary,
+            (
+                const rotorDiskSource& rotor,
+                const dictionary& dict
+            ),
+            (rotor, dict)
+        );
+
+
+    // Constructors
+
+        //- Construct from components
+        trimModel
+        (
+            const rotorDiskSource& rotor,
+            const dictionary& dict,
+            const word& name
+        );
+
+
+    // Selectors
+
+        //- Return a reference to the selected trim model
+        static autoPtr<trimModel> New
+        (
+            const rotorDiskSource& rotor,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~trimModel();
+
+
+    // Member functions
+
+        //- Read
+        virtual void read(const dictionary& dict);
+
+        //- Return the geometric angle of attack [rad]
+        virtual tmp<scalarField> thetag() const = 0;
+
+        //- Correct the model
+        virtual void correct(const vectorField& U, vectorField& force) = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/fieldSources/basicSource/rotorDiskSource/trimModel/trimModel/trimModelNew.C

@@ -0,0 +1,59 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "trimModel.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::trimModel> Foam::trimModel::New
+(
+    const rotorDiskSource& rotor,
+    const dictionary& dict
+)
+{
+    const word modelType(dict.lookup(typeName));
+
+    Info<< "    Selecting " << typeName << " " << modelType << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(modelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "trimModel::New(const rotorDiskSource&, const dictionary&)"
+        )   << "Unknown " << typeName << " type "
+            << modelType << nl << nl
+            << "Valid " << typeName << " types are:" << nl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<trimModel>(cstrIter()(rotor, dict));
+}
+
+
+// ************************************************************************* //

wmake/rules/General/mplibMPICH

@@ -1,3 +1,3 @@
 PFLAGS     = -DMPICH_SKIP_MPICXX
 PINC       = -I$(MPI_ARCH_PATH)/include
-PLIBS      = -L$(MPI_ARCH_PATH)/lib -lmpich -lrt
+PLIBS      = -L$(MPI_ARCH_PATH)/lib$(WM_COMPILER_LIB_ARCH) -L$(MPI_ARCH_PATH)/lib -lmpich -lrt

wmake/rules/General/mplibMPICH-GM

@@ -1,3 +1,3 @@
 PFLAGS     =
 PINC       = -I$(MPI_ARCH_PATH)/include
-PLIBS      = -L$(MPI_ARCH_PATH)/lib -lmpich -L$(GM_LIB_PATH) -lgm
+PLIBS      = -L$(MPI_ARCH_PATH)/lib$(WM_COMPILER_LIB_ARCH) -L$(MPI_ARCH_PATH)/lib -lmpich -L$(GM_LIB_PATH) -lgm

wmake/rules/General/mplibOPENMPI

@@ -1,3 +1,3 @@
 PFLAGS     = -DOMPI_SKIP_MPICXX
 PINC       = -I$(MPI_ARCH_PATH)/include
-PLIBS      = -L$(MPI_ARCH_PATH)/lib -lmpi
+PLIBS      = -L$(MPI_ARCH_PATH)/lib$(WM_COMPILER_LIB_ARCH) -L$(MPI_ARCH_PATH)/lib -lmpi