ENH: Changed mixed coupled bc to use combination of zerogradient and fixedValue

Is more convergent&stable than switching.

src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentTemperatureCoupledBaffleMixed/turbulentTemperatureCoupledBaffleMixedFvPatchScalarField.C

@@ -42,71 +42,6 @@ namespace compressible
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-bool turbulentTemperatureCoupledBaffleMixedFvPatchScalarField::interfaceOwner
-(
-    const polyMesh& nbrRegion,
-    const polyPatch& nbrPatch
-) const
-{
-    const fvMesh& myRegion = patch().boundaryMesh().mesh();
-
-    if (nbrRegion.name() == myRegion.name())
-    {
-        return patch().index() < nbrPatch.index();
-    }
-    else
-    {
-        const regionProperties& props =
-            myRegion.objectRegistry::parent().lookupObject<regionProperties>
-            (
-                "regionProperties"
-            );
-
-        label myIndex = findIndex(props.fluidRegionNames(), myRegion.name());
-        if (myIndex == -1)
-        {
-            label i = findIndex(props.solidRegionNames(), myRegion.name());
-
-            if (i == -1)
-            {
-                FatalErrorIn
-                (
-                    "turbulentTemperatureCoupledBaffleMixedFvPatchScalarField"
-                    "::interfaceOwner(const polyMesh&"
-                    ", const polyPatch&)const"
-                )   << "Cannot find region " << myRegion.name()
-                    << " neither in fluids " << props.fluidRegionNames()
-                    << " nor in solids " << props.solidRegionNames()
-                    << exit(FatalError);
-            }
-            myIndex = props.fluidRegionNames().size() + i;
-        }
-        label nbrIndex = findIndex
-        (
-            props.fluidRegionNames(),
-            nbrRegion.name()
-        );
-        if (nbrIndex == -1)
-        {
-            label i = findIndex(props.solidRegionNames(), nbrRegion.name());
-
-            if (i == -1)
-            {
-                FatalErrorIn
-                (
-                    "coupleManager::interfaceOwner"
-                    "(const polyMesh&, const polyPatch&) const"
-                )   << "Cannot find region " << nbrRegion.name()
-                    << " neither in fluids " << props.fluidRegionNames()
-                    << " nor in solids " << props.solidRegionNames()
-                    << exit(FatalError);
-            }
-            nbrIndex = props.fluidRegionNames().size() + i;
-        }
-        return myIndex < nbrIndex;
-    }
-}
-
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -124,7 +59,6 @@ turbulentTemperatureCoupledBaffleMixedFvPatchScalarField
     this->refValue() = 0.0;
     this->refGrad() = 0.0;
     this->valueFraction() = 1.0;
-    this->fixesValue_ = true;
 }
 
 
@@ -139,8 +73,7 @@ turbulentTemperatureCoupledBaffleMixedFvPatchScalarField
 :
     mixedFvPatchScalarField(ptf, p, iF, mapper),
     neighbourFieldName_(ptf.neighbourFieldName_),
-    KName_(ptf.KName_),
-    fixesValue_(ptf.fixesValue_)
+    KName_(ptf.KName_)
 {}
 
 
@@ -183,7 +116,6 @@ turbulentTemperatureCoupledBaffleMixedFvPatchScalarField
         refValue() = scalarField("refValue", dict, p.size());
         refGrad() = scalarField("refGradient", dict, p.size());
         valueFraction() = scalarField("valueFraction", dict, p.size());
-        fixesValue_ = readBool(dict.lookup("fixesValue"));
     }
     else
     {
@@ -191,7 +123,6 @@ turbulentTemperatureCoupledBaffleMixedFvPatchScalarField
         refValue() = *this;
         refGrad() = 0.0;
         valueFraction() = 1.0;
-        fixesValue_ = true;
     }
 }
 
@@ -205,8 +136,7 @@ turbulentTemperatureCoupledBaffleMixedFvPatchScalarField
 :
     mixedFvPatchScalarField(wtcsf, iF),
     neighbourFieldName_(wtcsf.neighbourFieldName_),
-    KName_(wtcsf.KName_),
-    fixesValue_(wtcsf.fixesValue_)
+    KName_(wtcsf.KName_)
 {}
 
 
@@ -281,95 +211,80 @@ void turbulentTemperatureCoupledBaffleMixedFvPatchScalarField::updateCoeffs()
 
     // Force recalculation of mapping and schedule
     const mapDistribute& distMap = mpp.map();
-    (void)distMap.schedule();
 
     tmp<scalarField> intFld = patchInternalField();
 
-    if (interfaceOwner(nbrMesh, nbrPatch.patch()))
-    {
-        // Note: other side information could be cached - it only needs
-        // to be updated the first time round the iteration (i.e. when
-        // switching regions) but unfortunately we don't have this information.
-
 
-        // Calculate the temperature by harmonic averaging
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Calculate the temperature by harmonic averaging
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-        const turbulentTemperatureCoupledBaffleMixedFvPatchScalarField&
-        nbrField =
-        refCast
-        <
-            const turbulentTemperatureCoupledBaffleMixedFvPatchScalarField
-        >
-        (
-            nbrPatch.lookupPatchField<volScalarField, scalar>
-            (
-                neighbourFieldName_
-            )
-        );
-
-        // Swap to obtain full local values of neighbour internal field
-        scalarField nbrIntFld = nbrField.patchInternalField();
-        mapDistribute::distribute
-        (
-            Pstream::defaultCommsType,
-            distMap.schedule(),
-            distMap.constructSize(),
-            distMap.subMap(),           // what to send
-            distMap.constructMap(),     // what to receive
-            nbrIntFld
-        );
-
-        // Swap to obtain full local values of neighbour K*delta
-        scalarField nbrKDelta = nbrField.K()*nbrPatch.deltaCoeffs();
-        mapDistribute::distribute
-        (
-            Pstream::defaultCommsType,
-            distMap.schedule(),
-            distMap.constructSize(),
-            distMap.subMap(),           // what to send
-            distMap.constructMap(),     // what to receive
-            nbrKDelta
-        );
-
-        tmp<scalarField> myKDelta = K()*patch().deltaCoeffs();
-
-        // Calculate common wall temperature. Reuse *this to store common value.
-        scalarField Twall
-        (
-            (myKDelta()*intFld() + nbrKDelta*nbrIntFld)
-          / (myKDelta() + nbrKDelta)
-        );
-        // Assign to me
-        fvPatchScalarField::operator=(Twall);
-        // Distribute back and assign to neighbour
-        mapDistribute::distribute
-        (
-            Pstream::defaultCommsType,
-            distMap.schedule(),
-            nbrField.size(),
-            distMap.constructMap(),     // reverse : what to send
-            distMap.subMap(),
-            Twall
-        );
-        const_cast<turbulentTemperatureCoupledBaffleMixedFvPatchScalarField&>
+    const turbulentTemperatureCoupledBaffleMixedFvPatchScalarField& nbrField =
+    refCast
+    <
+        const turbulentTemperatureCoupledBaffleMixedFvPatchScalarField
+    >
+    (
+        nbrPatch.lookupPatchField<volScalarField, scalar>
         (
-            nbrField
-        ).fvPatchScalarField::operator=(Twall);
-    }
+            neighbourFieldName_
+        )
+    );
+
+    // Swap to obtain full local values of neighbour internal field
+    scalarField nbrIntFld = nbrField.patchInternalField();
+    mapDistribute::distribute
+    (
+        Pstream::defaultCommsType,
+        distMap.schedule(),
+        distMap.constructSize(),
+        distMap.subMap(),           // what to send
+        distMap.constructMap(),     // what to receive
+        nbrIntFld
+    );
 
+    // Swap to obtain full local values of neighbour K*delta
+    scalarField nbrKDelta = nbrField.K()*nbrPatch.deltaCoeffs();
+    mapDistribute::distribute
+    (
+        Pstream::defaultCommsType,
+        distMap.schedule(),
+        distMap.constructSize(),
+        distMap.subMap(),           // what to send
+        distMap.constructMap(),     // what to receive
+        nbrKDelta
+    );
 
-    // Switch between fixed value (of harmonic avg) or gradient
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    tmp<scalarField> myKDelta = K()*patch().deltaCoeffs();
 
-    label nFixed = 0;
 
-    // Like snGrad but bypass switching on refValue/refGrad.
-    tmp<scalarField> normalGradient = (*this-intFld())*patch().deltaCoeffs();
+    // Both sides agree on
+    // - temperature : (myKDelta*fld + nbrKDelta*nbrFld)/(myKDelta+nbrKDelta)
+    // - gradient    : (temperature-fld)*delta
+    // We've got a degree of freedom in how to implement this in a mixed bc.
+    // (what gradient, what fixedValue and mixing coefficient)
+    // Two reasonable choices:
+    // 1. specify above temperature on one side (preferentially the high side)
+    //    and above gradient on the other. So this will switch between pure
+    //    fixedvalue and pure fixedgradient
+    // 2. specify gradient and temperature such that the equations are the
+    //    same on both sides. This leads to the choice of
+    //    - refGradient = zero gradient
+    //    - refValue = neighbour value
+    //    - mixFraction = nbrKDelta / (nbrKDelta + myKDelta())
+
+
+    this->refValue() = nbrIntFld;
+
+    this->refGrad() = 0.0;
+
+    this->valueFraction() = nbrKDelta / (nbrKDelta + myKDelta());
+
+    mixedFvPatchScalarField::updateCoeffs();
+
 
     if (debug)
     {
-        scalar Q = gSum(K()*patch().magSf()*normalGradient());
+        scalar Q = gSum(K()*patch().magSf()*snGrad());
 
         Info<< patch().boundaryMesh().mesh().name() << ':'
             << patch().name() << ':'
@@ -384,48 +299,6 @@ void turbulentTemperatureCoupledBaffleMixedFvPatchScalarField::updateCoeffs()
             << " avg:" << gAverage(*this)
             << endl;
     }
-
-    forAll(*this, i)
-    {
-        // if outgoing flux use fixed value.
-        if (normalGradient()[i] < 0.0)
-        {
-            this->refValue()[i] = operator[](i);
-            this->refGrad()[i] = 0.0;   // not used by me
-            this->valueFraction()[i] = 1.0;
-            nFixed++;
-        }
-        else
-        {
-            // Fixed gradient. Make sure to have valid refValue (even though
-            // I am not using it - other boundary conditions might)
-            this->refValue()[i] = operator[](i);
-            this->refGrad()[i] = normalGradient()[i];
-            this->valueFraction()[i] = 0.0;
-        }
-    }
-
-    reduce(nFixed, sumOp<label>());
-
-    fixesValue_ = (nFixed > 0);
-
-    if (debug)
-    {
-        label nTotSize = returnReduce(this->size(), sumOp<label>());
-
-        Info<< patch().boundaryMesh().mesh().name() << ':'
-            << patch().name() << ':'
-            << this->dimensionedInternalField().name() << " -> "
-            << nbrMesh.name() << ':'
-            << nbrPatch.name() << ':'
-            << this->dimensionedInternalField().name() << " :"
-            << " patch:" << patch().name()
-            << " out of:" << nTotSize
-            << " fixedBC:" << nFixed
-            << " gradient:" << nTotSize-nFixed << endl;
-    }
-
-    mixedFvPatchScalarField::updateCoeffs();
 }
 
 
@@ -438,7 +311,6 @@ void turbulentTemperatureCoupledBaffleMixedFvPatchScalarField::write
     os.writeKeyword("neighbourFieldName")<< neighbourFieldName_
         << token::END_STATEMENT << nl;
     os.writeKeyword("K") << KName_ << token::END_STATEMENT << nl;
-    os.writeKeyword("fixesValue") << fixesValue_ << token::END_STATEMENT << nl;
 }
 
 

src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentTemperatureCoupledBaffleMixed/turbulentTemperatureCoupledBaffleMixedFvPatchScalarField.H

@@ -92,14 +92,6 @@ class turbulentTemperatureCoupledBaffleMixedFvPatchScalarField
         //- Name of thermal conductivity field
         const word KName_;
 
-        bool fixesValue_;
-
-
-    // Private Member Functions
-
-        //- Am I or neighbour owner of interface
-        bool interfaceOwner(const polyMesh&, const polyPatch&) const;
-
 
 public:
 
@@ -176,14 +168,6 @@ public:
         //- Get corresponding K field
         tmp<scalarField> K() const;
 
-        //- Return true if this patch field fixes a value.
-        //  Needed to check if a level has to be specified while solving
-        //  Poissons equations.
-        virtual bool fixesValue() const
-        {
-            return fixesValue_;
-        }
-
         //- Update the coefficients associated with the patch field
         virtual void updateCoeffs();