src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchField.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 Author
  26     Hrvoje Jasak, Wikki Ltd.  All rights reserved
  27
  28 \*---------------------------------------------------------------------------*/
  29
  30 #include "regionCoupleFvPatchField.H"
  31
  32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  33
  34 namespace Foam
  35 {
  36
  37 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  38
  39 template<class Type>
  40 regionCoupleFvPatchField<Type>::regionCoupleFvPatchField
  41 (
  42     const fvPatch& p,
  43     const DimensionedField<Type, volMesh>& iF
  44 )
  45 :
  46     coupledFvPatchField<Type>(p, iF),
  47     regionCouplePatch_(refCast<const regionCoupleFvPatch>(p)),
  48     remoteFieldName_(iF.name()),
  49     matrixUpdateBuffer_()
  50 {}
  51
  52
  53 template<class Type>
  54 regionCoupleFvPatchField<Type>::regionCoupleFvPatchField
  55 (
  56     const fvPatch& p,
  57     const DimensionedField<Type, volMesh>& iF,
  58     const dictionary& dict
  59 )
  60 :
  61     coupledFvPatchField<Type>(p, iF, dict),
  62     regionCouplePatch_(refCast<const regionCoupleFvPatch>(p)),
  63     remoteFieldName_(dict.lookup("remoteField")),
  64     matrixUpdateBuffer_()
  65 {
  66     if (!isType<regionCoupleFvPatch>(p))
  67     {
  68         FatalIOErrorIn
  69         (
  70             "regionCoupleFvPatchField<Type>::regionCoupleFvPatchField\n"
  71             "(\n"
  72             "    const fvPatch& p,\n"
  73             "    const DimensionedField<Type, volMesh>& iF,\n"
  74             "    const dictionary& dict\n"
  75             ")\n",
  76             dict
  77         )   << "patch " << this->patch().index() << " not regionCouple type. "
  78             << "Patch type = " << p.type()
  79             << exit(FatalIOError);
  80     }
  81 }
  82
  83
  84 template<class Type>
  85 regionCoupleFvPatchField<Type>::regionCoupleFvPatchField
  86 (
  87     const regionCoupleFvPatchField<Type>& ptf,
  88     const fvPatch& p,
  89     const DimensionedField<Type, volMesh>& iF,
  90     const fvPatchFieldMapper& mapper
  91 )
  92 :
  93     coupledFvPatchField<Type>(ptf, p, iF, mapper),
  94     regionCouplePatch_(refCast<const regionCoupleFvPatch>(p)),
  95     remoteFieldName_(ptf.remoteFieldName_),
  96     matrixUpdateBuffer_()
  97 {
  98     if (!isType<regionCoupleFvPatch>(this->patch()))
  99     {
 100         FatalErrorIn
 101         (
 102             "regionCoupleFvPatchField<Type>::regionCoupleFvPatchField\n"
 103             "(\n"
 104             "    const regionCoupleFvPatchField<Type>& ptf,\n"
 105             "    const fvPatch& p,\n"
 106             "    const DimensionedField<Type, volMesh>& iF,\n"
 107             "    const fvPatchFieldMapper& mapper\n"
 108             ")\n"
 109         )   << "Field type does not correspond to patch type for patch "
 110             << this->patch().index() << "." << endl
 111             << "Field type: " << typeName << endl
 112             << "Patch type: " << this->patch().type()
 113             << exit(FatalError);
 114     }
 115 }
 116
 117
 118 template<class Type>
 119 regionCoupleFvPatchField<Type>::regionCoupleFvPatchField
 120 (
 121     const regionCoupleFvPatchField<Type>& ptf,
 122     const DimensionedField<Type, volMesh>& iF
 123 )
 124 :
 125     regionCoupleLduInterfaceField(),
 126     coupledFvPatchField<Type>(ptf, iF),
 127     regionCouplePatch_(refCast<const regionCoupleFvPatch>(ptf.patch())),
 128     remoteFieldName_(ptf.remoteFieldName_),
 129     matrixUpdateBuffer_()
 130 {}
 131
 132
 133 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 134
 135 // Return neighbour field
 136 template<class Type>
 137 const regionCoupleFvPatchField<Type>&
 138 regionCoupleFvPatchField<Type>::shadowPatchField() const
 139 {
 140     // Lookup neighbour field
 141     typedef GeometricField<Type, fvPatchField, volMesh> GeoField;
 142
 143     const GeoField& coupleField =
 144         regionCouplePatch_.shadowRegion().
 145         objectRegistry::lookupObject<GeoField>(remoteFieldName_);
 146
 147     return refCast<const regionCoupleFvPatchField<Type> >
 148     (
 149         coupleField.boundaryField()[regionCouplePatch_.shadowIndex()]
 150     );
 151 }
 152
 153
 154 // Return neighbour field
 155 template<class Type>
 156 tmp<Field<Type> > regionCoupleFvPatchField<Type>::patchNeighbourField() const
 157 {
 158     return regionCouplePatch_.interpolate
 159     (
 160         shadowPatchField().patchInternalField()
 161     );
 162 }
 163
 164
 165 template<class Type>
 166 void regionCoupleFvPatchField<Type>::evaluate
 167 (
 168     const Pstream::commsTypes
 169 )
 170 {
 171     // Implement weights-based stabilised harmonic interpolation using
 172     // magnitude of type
 173     // Algorithm:
 174     // 1) calculate magnitude of internal field and neighbour field
 175     // 2) calculate harmonic mean magnitude
 176     // 3) express harmonic mean magnitude as: mean = w*mOwn + (1 - w)*mNei
 177     // 4) Based on above, calculate w = (mean - mNei)/(mOwn - mNei)
 178     // 5) Use weights to interpolate values
 179
 180     Field<Type> fOwn = this->patchInternalField();
 181     Field<Type> fNei = this->patchNeighbourField();
 182
 183     scalarField magFOwn = mag(fOwn);
 184     scalarField magFNei = mag(fNei);
 185
 186     // Calculate internal weights using field magnitude
 187     scalarField weights(fOwn.size());
 188
 189     forAll (weights, faceI)
 190     {
 191         scalar mOwn = magFOwn[faceI];
 192         scalar mNei = magFNei[faceI];
 193
 194         scalar den = mOwn - mNei;
 195
 196         if (mag(den) > SMALL)
 197         {
 198             scalar mean = 2.0*mOwn*mNei/(mOwn + mNei);
 199             weights[faceI] = (mean - mNei)/den;
 200         }
 201         else
 202         {
 203             weights[faceI] = 0.5;
 204         }
 205     }
 206
 207     // Do interpolation
 208     Field<Type>::operator=(weights*fOwn + (1.0 - weights)*fNei);
 209 }
 210
 211
 212 // Initialise neighbour processor internal cell data
 213 template<class Type>
 214 void regionCoupleFvPatchField<Type>::initInterfaceMatrixUpdate
 215 (
 216     const scalarField& psiInternal,
 217     scalarField&,
 218     const lduMatrix&,
 219     const scalarField&,
 220     const direction,
 221     const Pstream::commsTypes
 222 ) const
 223 {
 224     matrixUpdateBuffer_ = this->patch().patchInternalField(psiInternal);
 225 }
 226
 227
 228 // Return matrix product for coupled boundary
 229 template<class Type>
 230 void regionCoupleFvPatchField<Type>::updateInterfaceMatrix
 231 (
 232     const scalarField& psiInternal,
 233     scalarField& result,
 234     const lduMatrix&,
 235     const scalarField& coeffs,
 236     const direction cmpt,
 237     const Pstream::commsTypes
 238 ) const
 239 {
 240     scalarField pnf =
 241         regionCouplePatch_.interpolate
 242         (
 243             this->shadowPatchField().matrixUpdateBuffer()
 244         );
 245
 246     // Multiply the field by coefficients and add into the result
 247     const unallocLabelList& fc = regionCouplePatch_.faceCells();
 248
 249     forAll(fc, elemI)
 250     {
 251         result[fc[elemI]] -= coeffs[elemI]*pnf[elemI];
 252     }
 253 }
 254
 255
 256 // Write
 257 template<class Type>
 258 void regionCoupleFvPatchField<Type>::write(Ostream& os) const
 259 {
 260     fvPatchField<Type>::write(os);
 261     os.writeKeyword("remoteField")
 262         << remoteFieldName_ << token::END_STATEMENT << nl;
 263     this->writeEntry("value", os);
 264 }
 265
 266
 267 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 268
 269 } // End namespace Foam
 270
 271 // ************************************************************************* //