diff --git a/applications/solvers/combustion/XiFoam/XiDyMFoam/Make/files b/applications/solvers/combustion/XiFoam/XiDyMFoam/Make/files
new file mode 100644
index 0000000000..55ea77cbdf
--- /dev/null
+++ b/applications/solvers/combustion/XiFoam/XiDyMFoam/Make/files
@@ -0,0 +1,3 @@
+XiDyMFoam.C
+
+EXE = $(FOAM_APPBIN)/XiDyMFoam
diff --git a/applications/solvers/combustion/XiFoam/XiDyMFoam/Make/options b/applications/solvers/combustion/XiFoam/XiDyMFoam/Make/options
new file mode 100644
index 0000000000..a726de69b3
--- /dev/null
+++ b/applications/solvers/combustion/XiFoam/XiDyMFoam/Make/options
@@ -0,0 +1,33 @@
+EXE_INC = \
+ -I.. \
+ -I$(LIB_SRC)/finiteVolume/lnInclude \
+ -I$(LIB_SRC)/fvOptions/lnInclude \
+ -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+ -I$(LIB_SRC)/dynamicMesh/lnInclude \
+ -I$(LIB_SRC)/meshTools/lnInclude \
+ -I$(LIB_SRC)/sampling/lnInclude \
+ -I$(LIB_SRC)/engine/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+ -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+ -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+ -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+ -I$(LIB_SRC)/thermophysicalModels/laminarFlameSpeed/lnInclude
+
+EXE_LIBS = \
+ -lfiniteVolume \
+ -lfvOptions \
+ -ldynamicFvMesh \
+ -ltopoChangerFvMesh \
+ -ldynamicMesh \
+ -lmeshTools \
+ -lsampling \
+ -lengine \
+ -lturbulenceModels \
+ -lcompressibleTurbulenceModels \
+ -lcompressibleTransportModels \
+ -lfluidThermophysicalModels \
+ -lreactionThermophysicalModels \
+ -lspecie \
+ -llaminarFlameSpeedModels
diff --git a/applications/solvers/combustion/XiFoam/XiDyMFoam/XiDyMFoam.C b/applications/solvers/combustion/XiFoam/XiDyMFoam/XiDyMFoam.C
new file mode 100644
index 0000000000..d4cb0666e2
--- /dev/null
+++ b/applications/solvers/combustion/XiFoam/XiDyMFoam/XiDyMFoam.C
@@ -0,0 +1,180 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 .
+
+Application
+ XiFoam
+
+Description
+ Solver for compressible premixed/partially-premixed combustion with
+ turbulence modelling.
+
+ Combusting RANS code using the b-Xi two-equation model.
+ Xi may be obtained by either the solution of the Xi transport
+ equation or from an algebraic exression. Both approaches are
+ based on Gulder's flame speed correlation which has been shown
+ to be appropriate by comparison with the results from the
+ spectral model.
+
+ Strain effects are encorporated directly into the Xi equation
+ but not in the algebraic approximation. Further work need to be
+ done on this issue, particularly regarding the enhanced removal rate
+ caused by flame compression. Analysis using results of the spectral
+ model will be required.
+
+ For cases involving very lean Propane flames or other flames which are
+ very strain-sensitive, a transport equation for the laminar flame
+ speed is present. This equation is derived using heuristic arguments
+ involving the strain time scale and the strain-rate at extinction.
+ the transport velocity is the same as that for the Xi equation.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "dynamicFvMesh.H"
+#include "psiuReactionThermo.H"
+#include "turbulentFluidThermoModel.H"
+#include "laminarFlameSpeed.H"
+#include "ignition.H"
+#include "Switch.H"
+#include "pimpleControl.H"
+#include "CorrectPhi.H"
+#include "fvIOoptionList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+ #include "setRootCase.H"
+
+ #include "createTime.H"
+ #include "createDynamicFvMesh.H"
+ #include "initContinuityErrs.H"
+
+ pimpleControl pimple(mesh);
+
+ #include "readCombustionProperties.H"
+ #include "readGravitationalAcceleration.H"
+ #include "createFields.H"
+ #include "createMRF.H"
+ #include "createFvOptions.H"
+ #include "createRhoUf.H"
+ #include "createControls.H"
+ #include "initContinuityErrs.H"
+ #include "compressibleCourantNo.H"
+ #include "setInitialDeltaT.H"
+
+ // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ Info<< "\nStarting time loop\n" << endl;
+
+ while (runTime.run())
+ {
+ #include "createTimeControls.H"
+
+ {
+ // Store divrhoU from the previous mesh so that it can be mapped
+ // and used in correctPhi to ensure the corrected phi has the
+ // same divergence
+ volScalarField divrhoU
+ (
+ "divrhoU",
+ fvc::div(fvc::absolute(phi, rho, U))
+ );
+
+ #include "compressibleCourantNo.H"
+ #include "setDeltaT.H"
+
+ runTime++;
+
+ Info<< "Time = " << runTime.timeName() << nl << endl;
+
+ // Store momentum to set rhoUf for introduced faces.
+ volVectorField rhoU("rhoU", rho*U);
+
+ // Do any mesh changes
+ mesh.update();
+
+ if (mesh.changing() && correctPhi)
+ {
+ // Calculate absolute flux from the mapped surface velocity
+ phi = mesh.Sf() & rhoUf;
+
+ #include "correctPhi.H"
+
+ // Make the fluxes relative to the mesh-motion
+ fvc::makeRelative(phi, rho, U);
+ }
+ }
+
+ if (mesh.changing() && checkMeshCourantNo)
+ {
+ #include "meshCourantNo.H"
+ }
+
+ #include "rhoEqn.H"
+ Info<< "rhoEqn max/min : " << max(rho).value()
+ << " " << min(rho).value() << endl;
+
+ // --- Pressure-velocity PIMPLE corrector loop
+ while (pimple.loop())
+ {
+ #include "UEqn.H"
+
+ #include "ftEqn.H"
+ #include "bEqn.H"
+ #include "EauEqn.H"
+ #include "EaEqn.H"
+
+ if (!ign.ignited())
+ {
+ thermo.heu() == thermo.he();
+ }
+
+ // --- Pressure corrector loop
+ while (pimple.correct())
+ {
+ #include "pEqn.H"
+ }
+
+ if (pimple.turbCorr())
+ {
+ turbulence->correct();
+ }
+ }
+
+ rho = thermo.rho();
+
+ runTime.write();
+
+ Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+ << " ClockTime = " << runTime.elapsedClockTime() << " s"
+ << nl << endl;
+ }
+
+ Info<< "End\n" << endl;
+
+ return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/combustion/XiFoam/XiDyMFoam/correctPhi.H b/applications/solvers/combustion/XiFoam/XiDyMFoam/correctPhi.H
new file mode 100644
index 0000000000..37072312ff
--- /dev/null
+++ b/applications/solvers/combustion/XiFoam/XiDyMFoam/correctPhi.H
@@ -0,0 +1,11 @@
+CorrectPhi
+(
+ U,
+ phi,
+ p,
+ rho,
+ psi,
+ dimensionedScalar("rAUf", dimTime, 1),
+ divrhoU,
+ pimple
+);
diff --git a/applications/solvers/combustion/XiFoam/XiDyMFoam/createControls.H b/applications/solvers/combustion/XiFoam/XiDyMFoam/createControls.H
new file mode 100644
index 0000000000..aed0e76956
--- /dev/null
+++ b/applications/solvers/combustion/XiFoam/XiDyMFoam/createControls.H
@@ -0,0 +1,11 @@
+#include "createTimeControls.H"
+
+bool correctPhi
+(
+ pimple.dict().lookupOrDefault("correctPhi", true)
+);
+
+bool checkMeshCourantNo
+(
+ pimple.dict().lookupOrDefault("checkMeshCourantNo", false)
+);
diff --git a/applications/solvers/combustion/XiFoam/XiDyMFoam/pEqn.H b/applications/solvers/combustion/XiFoam/XiDyMFoam/pEqn.H
new file mode 100644
index 0000000000..438396ae90
--- /dev/null
+++ b/applications/solvers/combustion/XiFoam/XiDyMFoam/pEqn.H
@@ -0,0 +1,99 @@
+rho = thermo.rho();
+
+volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn.H();
+
+if (pimple.transonic())
+{
+ surfaceScalarField phid
+ (
+ "phid",
+ fvc::interpolate(psi)
+ *(
+ (fvc::interpolate(HbyA) & mesh.Sf())
+ + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
+ )
+ );
+
+ fvc::makeRelative(phid, psi, U);
+ MRF.makeRelative(fvc::interpolate(psi), phid);
+
+ while (pimple.correctNonOrthogonal())
+ {
+ fvScalarMatrix pEqn
+ (
+ fvm::ddt(psi, p)
+ + fvm::div(phid, p)
+ - fvm::laplacian(rhorAUf, p)
+ ==
+ fvOptions(psi, p, rho.name())
+ );
+
+ pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+
+ if (pimple.finalNonOrthogonalIter())
+ {
+ phi == pEqn.flux();
+ }
+ }
+}
+else
+{
+ surfaceScalarField phiHbyA
+ (
+ "phiHbyA",
+ (
+ (fvc::interpolate(rho*HbyA) & mesh.Sf())
+ + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
+ )
+ );
+
+ fvc::makeRelative(phiHbyA, rho, U);
+ MRF.makeRelative(phiHbyA);
+
+ while (pimple.correctNonOrthogonal())
+ {
+ fvScalarMatrix pEqn
+ (
+ fvm::ddt(psi, p)
+ + fvc::div(phiHbyA)
+ - fvm::laplacian(rhorAUf, p)
+ ==
+ fvOptions(psi, p, rho.name())
+ );
+
+ pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+
+ if (pimple.finalNonOrthogonalIter())
+ {
+ phi = phiHbyA + pEqn.flux();
+ }
+ }
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+U = HbyA - rAU*fvc::grad(p);
+U.correctBoundaryConditions();
+fvOptions.correct(U);
+K = 0.5*magSqr(U);
+
+{
+ rhoUf = fvc::interpolate(rho*U);
+ surfaceVectorField n(mesh.Sf()/mesh.magSf());
+ rhoUf += n*(fvc::absolute(phi, rho, U)/mesh.magSf() - (n & rhoUf));
+}
+
+if (thermo.dpdt())
+{
+ dpdt = fvc::ddt(p);
+
+ if (mesh.moving())
+ {
+ dpdt -= fvc::div(fvc::meshPhi(rho, U), p);
+ }
+}
diff --git a/applications/solvers/combustion/XiFoam/XiDyMFoam/readControls.H b/applications/solvers/combustion/XiFoam/XiDyMFoam/readControls.H
new file mode 100644
index 0000000000..cf50980460
--- /dev/null
+++ b/applications/solvers/combustion/XiFoam/XiDyMFoam/readControls.H
@@ -0,0 +1,7 @@
+ #include "readTimeControls.H"
+
+ bool correctPhi =
+ pimple.dict().lookupOrDefault("correctPhi", true);
+
+ bool checkMeshCourantNo =
+ pimple.dict().lookupOrDefault("checkMeshCourantNo", false);
diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicAMI/cyclicAMIFvPatch.C b/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicAMI/cyclicAMIFvPatch.C
index 5d8935139d..0f6064471f 100644
--- a/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicAMI/cyclicAMIFvPatch.C
+++ b/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicAMI/cyclicAMIFvPatch.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -34,6 +34,13 @@ namespace Foam
{
defineTypeNameAndDebug(cyclicAMIFvPatch, 0);
addToRunTimeSelectionTable(fvPatch, cyclicAMIFvPatch, polyPatch);
+ addNamedToRunTimeSelectionTable
+ (
+ fvPatch,
+ cyclicAMIFvPatch,
+ polyPatch,
+ cyclicPeriodicAMI
+ );
}
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
index 0e7a762423..7bad5f58c8 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
@@ -327,9 +327,30 @@ void Foam::AMIInterpolation::agglomerate
// So now we have agglomeration of the target side in
// allRestrict:
- // 0..size-1 : local agglomeration (= targetRestrictAddressing)
+ // 0..size-1 : local agglomeration (= targetRestrictAddressing
+ // (but potentially permutated))
// size.. : agglomeration data from other processors
+
+ // The trickiness in this algorithm is finding out the compaction
+ // of the remote data (i.e. allocation of the coarse 'slots'). We could
+ // either send across the slot compaction maps or just make sure
+ // that we allocate the slots in exactly the same order on both sending
+ // and receiving side (e.g. if the submap is set up to send 4 items,
+ // the constructMap is also set up to receive 4 items.
+
+
+ // Short note about the various types of indices:
+ // - face indices : indices into the geometry.
+ // - coarse face indices : how the faces get agglomerated
+ // - transferred data : how mapDistribute sends/receives data,
+ // - slots : indices into data after distribution (e.g. stencil,
+ // srcAddress/tgtAddress). Note: for fully local addressing
+ // the slots are equal to face indices.
+ // A mapDistribute has:
+ // - a subMap : these are face indices
+ // - a constructMap : these are from 'transferred-date' to slots
+
labelListList tgtSubMap(Pstream::nProcs());
// Local subMap is just identity
@@ -341,11 +362,17 @@ void Foam::AMIInterpolation::agglomerate
{
if (procI != Pstream::myProcNo())
{
- // Combine entries that point to the same coarse element. All
- // the elements refer to local data so index into
- // targetRestrictAddressing or allRestrict (since the same
- // for local data).
+ // Combine entries that point to the same coarse element.
+ // The important bit is to loop over the data (and hand out
+ // compact indices ) in 'transferred data' order. This
+ // guarantees that we're doing exactly the
+ // same on sending and receiving side - e.g. the fourth element
+ // in the subMap is the fourth element received in the
+ // constructMap
+
const labelList& elems = map.subMap()[procI];
+ const labelList& elemsMap =
+ map.constructMap()[Pstream::myProcNo()];
labelList& newSubMap = tgtSubMap[procI];
newSubMap.setSize(elems.size());
@@ -354,7 +381,7 @@ void Foam::AMIInterpolation::agglomerate
forAll(elems, i)
{
- label fineElem = elems[i];
+ label fineElem = elemsMap[elems[i]];
label coarseElem = allRestrict[fineElem];
if (oldToNew[coarseElem] == -1)
{
@@ -372,16 +399,31 @@ void Foam::AMIInterpolation::agglomerate
// the sending map
labelListList tgtConstructMap(Pstream::nProcs());
- labelList tgtCompactMap;
// Local constructMap is just identity
{
tgtConstructMap[Pstream::myProcNo()] =
identity(targetCoarseSize);
-
- tgtCompactMap = targetRestrictAddressing;
}
- tgtCompactMap.setSize(map.constructSize());
+
+ labelList tgtCompactMap(map.constructSize());
+
+ {
+ // Note that in special cases (e.g. 'appending' two AMIs) the
+ // local size after distributing can be longer than the number
+ // of faces. I.e. it duplicates elements.
+ // Since we don't know this size instead we loop over all
+ // reachable elements (using the local constructMap)
+
+ const labelList& elemsMap = map.constructMap()[Pstream::myProcNo()];
+ forAll(elemsMap, i)
+ {
+ label fineElem = elemsMap[i];
+ label coarseElem = allRestrict[fineElem];
+ tgtCompactMap[fineElem] = coarseElem;
+ }
+ }
+
label compactI = targetCoarseSize;
// Compact data from other processors
@@ -440,7 +482,6 @@ void Foam::AMIInterpolation::agglomerate
}
}
-
srcAddress.setSize(sourceCoarseSize);
srcWeights.setSize(sourceCoarseSize);
@@ -493,7 +534,7 @@ void Foam::AMIInterpolation::agglomerate
forAll(fineSrcAddress, faceI)
{
// All the elements contributing to faceI. Are slots in
- // mapDistribute'd data.
+ // target data.
const labelList& elems = fineSrcAddress[faceI];
const scalarList& weights = fineSrcWeights[faceI];
const scalar fineArea = fineSrcMagSf[faceI];
@@ -1041,28 +1082,7 @@ void Foam::AMIInterpolation::update
);
// weights normalisation
- normaliseWeights
- (
- srcMagSf_,
- "source",
- srcAddress_,
- srcWeights_,
- srcWeightsSum_,
- AMIPtr->conformal(),
- true,
- lowWeightCorrection_
- );
- normaliseWeights
- (
- tgtMagSf_,
- "target",
- tgtAddress_,
- tgtWeights_,
- tgtWeightsSum_,
- AMIPtr->conformal(),
- true,
- lowWeightCorrection_
- );
+ normaliseWeights(AMIPtr->conformal(), true);
// cache maps and reset addresses
List > cMap;
@@ -1100,28 +1120,7 @@ void Foam::AMIInterpolation::update
tgtWeights_
);
- normaliseWeights
- (
- srcMagSf_,
- "source",
- srcAddress_,
- srcWeights_,
- srcWeightsSum_,
- AMIPtr->conformal(),
- true,
- lowWeightCorrection_
- );
- normaliseWeights
- (
- tgtMagSf_,
- "target",
- tgtAddress_,
- tgtWeights_,
- tgtWeightsSum_,
- AMIPtr->conformal(),
- true,
- lowWeightCorrection_
- );
+ normaliseWeights(AMIPtr->conformal(), true);
}
if (debug)
@@ -1137,6 +1136,217 @@ void Foam::AMIInterpolation::update
}
+template
+void Foam::AMIInterpolation::append
+(
+ const SourcePatch& srcPatch,
+ const TargetPatch& tgtPatch
+)
+{
+ // Create a new interpolation
+ autoPtr > newPtr
+ (
+ new AMIInterpolation
+ (
+ srcPatch,
+ tgtPatch,
+ triMode_,
+ requireMatch_,
+ methodName_,
+ lowWeightCorrection_,
+ reverseTarget_
+ )
+ );
+
+ // If parallel then combine the mapDistribution and re-index
+ if (singlePatchProc_ == -1)
+ {
+ labelListList& srcSubMap = srcMapPtr_->subMap();
+ labelListList& srcConstructMap = srcMapPtr_->constructMap();
+
+ labelListList& tgtSubMap = tgtMapPtr_->subMap();
+ labelListList& tgtConstructMap = tgtMapPtr_->constructMap();
+
+ labelListList& newSrcSubMap = newPtr->srcMapPtr_->subMap();
+ labelListList& newSrcConstructMap = newPtr->srcMapPtr_->constructMap();
+
+ labelListList& newTgtSubMap = newPtr->tgtMapPtr_->subMap();
+ labelListList& newTgtConstructMap = newPtr->tgtMapPtr_->constructMap();
+
+ // Re-calculate the source indices
+ {
+ labelList mapMap(0), newMapMap(0);
+ forAll(srcSubMap, procI)
+ {
+ mapMap.append
+ (
+ identity(srcConstructMap[procI].size())
+ + mapMap.size() + newMapMap.size()
+ );
+ newMapMap.append
+ (
+ identity(newSrcConstructMap[procI].size())
+ + mapMap.size() + newMapMap.size()
+ );
+ }
+
+ forAll(srcSubMap, procI)
+ {
+ forAll(srcConstructMap[procI], srcI)
+ {
+ srcConstructMap[procI][srcI] =
+ mapMap[srcConstructMap[procI][srcI]];
+ }
+ }
+
+ forAll(srcSubMap, procI)
+ {
+ forAll(newSrcConstructMap[procI], srcI)
+ {
+ newSrcConstructMap[procI][srcI] =
+ newMapMap[newSrcConstructMap[procI][srcI]];
+ }
+ }
+
+ forAll(tgtAddress_, tgtI)
+ {
+ forAll(tgtAddress_[tgtI], tgtJ)
+ {
+ tgtAddress_[tgtI][tgtJ] =
+ mapMap[tgtAddress_[tgtI][tgtJ]];
+ }
+ }
+
+ forAll(newPtr->tgtAddress_, tgtI)
+ {
+ forAll(newPtr->tgtAddress_[tgtI], tgtJ)
+ {
+ newPtr->tgtAddress_[tgtI][tgtJ] =
+ newMapMap[newPtr->tgtAddress_[tgtI][tgtJ]];
+ }
+ }
+ }
+
+ // Re-calculate the target indices
+ {
+ labelList mapMap(0), newMapMap(0);
+ forAll(srcSubMap, procI)
+ {
+ mapMap.append
+ (
+ identity(tgtConstructMap[procI].size())
+ + mapMap.size() + newMapMap.size()
+ );
+ newMapMap.append
+ (
+ identity(newTgtConstructMap[procI].size())
+ + mapMap.size() + newMapMap.size()
+ );
+ }
+
+ forAll(srcSubMap, procI)
+ {
+ forAll(tgtConstructMap[procI], tgtI)
+ {
+ tgtConstructMap[procI][tgtI] =
+ mapMap[tgtConstructMap[procI][tgtI]];
+ }
+ }
+
+ forAll(srcSubMap, procI)
+ {
+ forAll(newTgtConstructMap[procI], tgtI)
+ {
+ newTgtConstructMap[procI][tgtI] =
+ newMapMap[newTgtConstructMap[procI][tgtI]];
+ }
+ }
+
+ forAll(srcAddress_, srcI)
+ {
+ forAll(srcAddress_[srcI], srcJ)
+ {
+ srcAddress_[srcI][srcJ] =
+ mapMap[srcAddress_[srcI][srcJ]];
+ }
+ }
+
+ forAll(newPtr->srcAddress_, srcI)
+ {
+ forAll(newPtr->srcAddress_[srcI], srcJ)
+ {
+ newPtr->srcAddress_[srcI][srcJ] =
+ newMapMap[newPtr->srcAddress_[srcI][srcJ]];
+ }
+ }
+ }
+
+ // Sum the construction sizes
+ srcMapPtr_->constructSize() += newPtr->srcMapPtr_->constructSize();
+ tgtMapPtr_->constructSize() += newPtr->tgtMapPtr_->constructSize();
+
+ // Combine the maps
+ forAll(srcSubMap, procI)
+ {
+ srcSubMap[procI].append(newSrcSubMap[procI]);
+ srcConstructMap[procI].append(newSrcConstructMap[procI]);
+
+ tgtSubMap[procI].append(newTgtSubMap[procI]);
+ tgtConstructMap[procI].append(newTgtConstructMap[procI]);
+ }
+ }
+
+ // Combine new and current source data
+ forAll(srcMagSf_, srcFaceI)
+ {
+ srcAddress_[srcFaceI].append(newPtr->srcAddress()[srcFaceI]);
+ srcWeights_[srcFaceI].append(newPtr->srcWeights()[srcFaceI]);
+ srcWeightsSum_[srcFaceI] += newPtr->srcWeightsSum()[srcFaceI];
+ }
+
+ // Combine new and current target data
+ forAll(tgtMagSf_, tgtFaceI)
+ {
+ tgtAddress_[tgtFaceI].append(newPtr->tgtAddress()[tgtFaceI]);
+ tgtWeights_[tgtFaceI].append(newPtr->tgtWeights()[tgtFaceI]);
+ tgtWeightsSum_[tgtFaceI] += newPtr->tgtWeightsSum()[tgtFaceI];
+ }
+}
+
+
+template
+void Foam::AMIInterpolation::normaliseWeights
+(
+ const bool conformal,
+ const bool output
+)
+{
+ normaliseWeights
+ (
+ srcMagSf_,
+ "source",
+ srcAddress_,
+ srcWeights_,
+ srcWeightsSum_,
+ conformal,
+ output,
+ lowWeightCorrection_
+ );
+
+ normaliseWeights
+ (
+ tgtMagSf_,
+ "target",
+ tgtAddress_,
+ tgtWeights_,
+ tgtWeightsSum_,
+ conformal,
+ output,
+ lowWeightCorrection_
+ );
+}
+
+
template
template
void Foam::AMIInterpolation::interpolateToTarget
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H
index 27d077f126..67205b69a7 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -418,6 +418,16 @@ public:
const TargetPatch& tgtPatch
);
+ //- Append additional addressing and weights
+ void append
+ (
+ const SourcePatch& srcPatch,
+ const TargetPatch& tgtPatch
+ );
+
+ //- Normalise the weights
+ void normaliseWeights(const bool conformal, const bool output);
+
// Evaluation
diff --git a/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPointPatch/cyclicAMIPointPatch.C b/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPointPatch/cyclicAMIPointPatch.C
index 2b4e6b8ad3..98fe75a0b7 100644
--- a/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPointPatch/cyclicAMIPointPatch.C
+++ b/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPointPatch/cyclicAMIPointPatch.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -38,6 +38,13 @@ namespace Foam
cyclicAMIPointPatch,
polyPatch
);
+ addNamedToRunTimeSelectionTable
+ (
+ facePointPatch,
+ cyclicAMIPointPatch,
+ polyPatch,
+ cyclicPeriodicAMI
+ );
}
diff --git a/src/meshTools/AMIInterpolation/patches/cyclicPeriodicAMI/cyclicPeriodicAMIPolyPatch/cyclicPeriodicAMIPolyPatch.C b/src/meshTools/AMIInterpolation/patches/cyclicPeriodicAMI/cyclicPeriodicAMIPolyPatch/cyclicPeriodicAMIPolyPatch.C
new file mode 100644
index 0000000000..dd9b431c43
--- /dev/null
+++ b/src/meshTools/AMIInterpolation/patches/cyclicPeriodicAMI/cyclicPeriodicAMIPolyPatch/cyclicPeriodicAMIPolyPatch.C
@@ -0,0 +1,723 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2015 OpenCFD Ltd.
+ \\/ 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 .
+
+\*---------------------------------------------------------------------------*/
+
+#include "cyclicPeriodicAMIPolyPatch.H"
+#include "addToRunTimeSelectionTable.H"
+
+// For debugging
+#include "OBJstream.H"
+#include "PatchTools.H"
+#include "Time.H"
+//Note: cannot use vtkSurfaceWriter here - circular linkage
+//#include "vtkSurfaceWriter.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ defineTypeNameAndDebug(cyclicPeriodicAMIPolyPatch, 0);
+
+ addToRunTimeSelectionTable(polyPatch, cyclicPeriodicAMIPolyPatch, word);
+ addToRunTimeSelectionTable
+ (
+ polyPatch,
+ cyclicPeriodicAMIPolyPatch,
+ dictionary
+ );
+}
+
+
+// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
+
+void Foam::cyclicPeriodicAMIPolyPatch::syncTransforms() const
+{
+ if (owner())
+ {
+ // At this point we guarantee that the transformations have been
+ // updated. There is one particular case, where if the periodic patch
+ // locally has zero faces then its transformation will not be set. We
+ // need to synchronise the transforms over the zero-sized patches as
+ // well.
+ //
+ // We can't put the logic into cyclicPolyPatch as
+ // processorCyclicPolyPatch uses cyclicPolyPatch functionality.
+ // Synchronisation will fail because processor-type patches do not exist
+ // on all processors.
+ //
+ // The use in cyclicPeriodicAMI is special; we use the patch even
+ // though we have no faces in it. Ideally, in future, the transformation
+ // logic will be abstracted out, and we won't need a periodic patch
+ // here. Until then, we can just fix the behaviour of the zero-sized
+ // coupled patches here
+
+ // Get the periodic patch
+ const coupledPolyPatch& periodicPatch
+ (
+ refCast
+ (
+ boundaryMesh()[periodicPatchID()]
+ )
+ );
+
+ // If there are any zero-sized periodic patches
+ if (returnReduce((size() && !periodicPatch.size()), orOp()))
+ {
+ if (periodicPatch.separation().size() > 1)
+ {
+ FatalErrorIn
+ (
+ "cyclicPeriodicAMIPolyPatch::resetAMI"
+ "(const AMIPatchToPatchInterpolation::interpolationMethod&"
+ ") const"
+ ) << "Periodic patch " << periodicPatchName_
+ << " has non-uniform separation vector "
+ << periodicPatch.separation()
+ << "This is not allowed inside " << type()
+ << " patch " << name()
+ << exit(FatalError);
+ }
+
+ if (periodicPatch.forwardT().size() > 1)
+ {
+ FatalErrorIn
+ (
+ "cyclicPeriodicAMIPolyPatch::resetAMI"
+ "(const AMIPatchToPatchInterpolation::interpolationMethod&"
+ ") const"
+ ) << "Periodic patch " << periodicPatchName_
+ << " has non-uniform transformation tensor "
+ << periodicPatch.forwardT()
+ << "This is not allowed inside " << type()
+ << " patch " << name()
+ << exit(FatalError);
+ }
+
+ // Note that zero-sized patches will have zero-sized fields for the
+ // separation vector, forward and reverse transforms. These need
+ // replacing with the transformations from other processors.
+
+ // Parallel in this context refers to a parallel transformation,
+ // rather than a rotational transformation.
+
+ // Note that a cyclic with zero faces is considered parallel so
+ // explicitly check for that.
+ bool isParallel =
+ (
+ periodicPatch.size()
+ && periodicPatch.parallel()
+ );
+ reduce(isParallel, orOp());
+
+ if (isParallel)
+ {
+ // Sync a list of separation vectors
+ List sep(Pstream::nProcs());
+ sep[Pstream::myProcNo()] = periodicPatch.separation();
+ Pstream::gatherList(sep);
+ Pstream::scatterList(sep);
+
+ List coll(Pstream::nProcs());
+ coll[Pstream::myProcNo()] = periodicPatch.collocated();
+ Pstream::gatherList(coll);
+ Pstream::scatterList(coll);
+
+ // If locally we have zero faces pick the first one that has a
+ // separation vector
+ if (!periodicPatch.size())
+ {
+ forAll(sep, procI)
+ {
+ if (sep[procI].size())
+ {
+ const_cast
+ (
+ periodicPatch.separation()
+ ) = sep[procI];
+ const_cast
+ (
+ periodicPatch.collocated()
+ ) = coll[procI];
+
+ break;
+ }
+ }
+ }
+ }
+ else
+ {
+ // Sync a list of forward and reverse transforms
+ List forwardT(Pstream::nProcs());
+ forwardT[Pstream::myProcNo()] = periodicPatch.forwardT();
+ Pstream::gatherList(forwardT);
+ Pstream::scatterList(forwardT);
+
+ List reverseT(Pstream::nProcs());
+ reverseT[Pstream::myProcNo()] = periodicPatch.reverseT();
+ Pstream::gatherList(reverseT);
+ Pstream::scatterList(reverseT);
+
+ // If locally we have zero faces pick the first one that has a
+ // transformation vector
+ if (!periodicPatch.size())
+ {
+ forAll(forwardT, procI)
+ {
+ if (forwardT[procI].size())
+ {
+ const_cast
+ (
+ periodicPatch.forwardT()
+ ) = forwardT[procI];
+ const_cast
+ (
+ periodicPatch.reverseT()
+ ) = reverseT[procI];
+
+ break;
+ }
+ }
+ }
+ }
+ }
+ }
+}
+
+
+void Foam::cyclicPeriodicAMIPolyPatch::writeOBJ
+(
+ const primitivePatch& p,
+ OBJstream& str
+) const
+{
+ // Collect faces and points
+ pointField allPoints;
+ faceList allFaces;
+ labelList pointMergeMap;
+ PatchTools::gatherAndMerge
+ (
+ -1.0, // do not merge points
+ p,
+ allPoints,
+ allFaces,
+ pointMergeMap
+ );
+
+ if (Pstream::master())
+ {
+ // Write base geometry
+ str.write(allFaces, allPoints);
+ }
+}
+
+
+void Foam::cyclicPeriodicAMIPolyPatch::resetAMI
+(
+ const AMIPatchToPatchInterpolation::interpolationMethod& AMIMethod
+) const
+{
+ if (owner())
+ {
+ // Get the periodic patch
+ const coupledPolyPatch& periodicPatch
+ (
+ refCast
+ (
+ boundaryMesh()[periodicPatchID()]
+ )
+ );
+
+ // Synchronise the transforms
+ syncTransforms();
+
+ // Create copies of both patches' points, transformed to the owner
+ pointField thisPoints0(localPoints());
+ pointField nbrPoints0(neighbPatch().localPoints());
+ transformPosition(nbrPoints0);
+
+ // Reset the stored number of periodic transformations to a lower
+ // absolute value if possible
+ if (nSectors_ > 0)
+ {
+ if (nTransforms_ > nSectors_/2)
+ {
+ nTransforms_ -= nSectors_;
+ }
+ else if (nTransforms_ < - nSectors_/2)
+ {
+ nTransforms_ += nSectors_;
+ }
+ }
+
+ // Apply the stored number of periodic transforms
+ for (label i = 0; i < nTransforms_; ++ i)
+ {
+ periodicPatch.transformPosition(thisPoints0);
+ }
+ for (label i = 0; i > nTransforms_; -- i)
+ {
+ periodicPatch.transformPosition(nbrPoints0);
+ }
+
+ autoPtr ownStr;
+ autoPtr neiStr;
+ if (debug)
+ {
+ const Time& runTime = boundaryMesh().mesh().time();
+
+ fileName dir(runTime.rootPath()/runTime.globalCaseName());
+ fileName postfix("_" + runTime.timeName()+"_expanded.obj");
+
+ ownStr.reset(new OBJstream(dir/name() + postfix));
+ neiStr.reset(new OBJstream(dir/neighbPatch().name() + postfix));
+
+ Info<< "cyclicPeriodicAMIPolyPatch::resetAMI : patch:" << name()
+ << " writing accumulated AMI to " << ownStr().name()
+ << " and " << neiStr().name() << endl;
+ }
+
+ // Create another copy
+ pointField thisPoints(thisPoints0);
+ pointField nbrPoints(nbrPoints0);
+
+ // Create patches for all the points
+
+ // Source patch at initial location
+ const primitivePatch thisPatch0
+ (
+ SubList(localFaces(), size()),
+ thisPoints0
+ );
+ // Source patch that gets moved
+ primitivePatch thisPatch
+ (
+ SubList(localFaces(), size()),
+ thisPoints
+ );
+ // Target patch at initial location
+ const primitivePatch nbrPatch0
+ (
+ SubList(neighbPatch().localFaces(), neighbPatch().size()),
+ nbrPoints0
+ );
+ // Target patch that gets moved
+ primitivePatch nbrPatch
+ (
+ SubList(neighbPatch().localFaces(), neighbPatch().size()),
+ nbrPoints
+ );
+
+ // Construct a new AMI interpolation between the initial patch locations
+ AMIPtr_.reset
+ (
+ new AMIPatchToPatchInterpolation
+ (
+ thisPatch0,
+ nbrPatch0,
+ surfPtr(),
+ faceAreaIntersect::tmMesh,
+ false,
+ AMIPatchToPatchInterpolation::imPartialFaceAreaWeight,
+ AMILowWeightCorrection_,
+ AMIReverse_
+ )
+ );
+
+ // Number of geometry replications
+ label iter(0);
+ label nTransformsOld(nTransforms_);
+
+ if (ownStr.valid())
+ {
+ writeOBJ(thisPatch0, ownStr());
+ }
+ if (neiStr.valid())
+ {
+ writeOBJ(nbrPatch0, neiStr());
+ }
+
+
+ // Weight sum averages
+ scalar srcSum(gAverage(AMIPtr_->srcWeightsSum()));
+ scalar tgtSum(gAverage(AMIPtr_->tgtWeightsSum()));
+
+ // Direction (or rather side of AMI : this or nbr patch) of
+ // geometry replication
+ bool direction = nTransforms_ >= 0;
+
+ // Increase in the source weight sum for the last iteration in the
+ // opposite direction. If the current increase is less than this, the
+ // direction (= side of AMI to transform) is reversed.
+ // We switch the side to replicate instead of reversing the transform
+ // since at the coupledPolyPatch level there is no
+ // 'reverseTransformPosition' functionality.
+ scalar srcSumDiff = 0;
+
+ if (debug)
+ {
+ Info<< "cyclicPeriodicAMIPolyPatch::resetAMI : patch:" << name()
+ << " srcSum:" << srcSum
+ << " tgtSum:" << tgtSum
+ << " direction:" << direction
+ << endl;
+ }
+
+ // Loop, replicating the geometry
+ while
+ (
+ (iter < maxIter_)
+ && (
+ (1 - srcSum > matchTolerance())
+ || (1 - tgtSum > matchTolerance())
+ )
+ )
+ {
+ if (direction)
+ {
+ periodicPatch.transformPosition(thisPoints);
+
+ if (debug)
+ {
+ Info<< "cyclicPeriodicAMIPolyPatch::resetAMI : patch:"
+ << name()
+ << " moving this side from:"
+ << gAverage(thisPatch.points())
+ << " to:" << gAverage(thisPoints) << endl;
+ }
+
+ thisPatch.movePoints(thisPoints);
+
+ if (debug)
+ {
+ Info<< "cyclicPeriodicAMIPolyPatch::resetAMI : patch:"
+ << name()
+ << " appending weights with untransformed slave side"
+ << endl;
+ }
+
+ AMIPtr_->append(thisPatch, nbrPatch0);
+
+ if (ownStr.valid())
+ {
+ writeOBJ(thisPatch, ownStr());
+ }
+ }
+ else
+ {
+ periodicPatch.transformPosition(nbrPoints);
+
+ if (debug)
+ {
+ Info<< "cyclicPeriodicAMIPolyPatch::resetAMI : patch:"
+ << name()
+ << " moving neighbour side from:"
+ << gAverage(nbrPatch.points())
+ << " to:" << gAverage(nbrPoints) << endl;
+ }
+
+ nbrPatch.movePoints(nbrPoints);
+
+ AMIPtr_->append(thisPatch0, nbrPatch);
+
+ if (neiStr.valid())
+ {
+ writeOBJ(nbrPatch, neiStr());
+ }
+ }
+
+ const scalar srcSumNew = gAverage(AMIPtr_->srcWeightsSum());
+ const scalar srcSumDiffNew = srcSumNew - srcSum;
+
+ if (srcSumDiffNew < srcSumDiff || srcSumDiffNew < SMALL)
+ {
+ direction = !direction;
+
+ srcSumDiff = srcSumDiffNew;
+ }
+
+ srcSum = srcSumNew;
+ tgtSum = gAverage(AMIPtr_->tgtWeightsSum());
+
+ nTransforms_ += direction ? +1 : -1;
+
+ ++ iter;
+
+ if (debug)
+ {
+ Info<< "cyclicPeriodicAMIPolyPatch::resetAMI : patch:" << name()
+ << " iteration:" << iter
+ << " srcSum:" << srcSum
+ << " tgtSum:" << tgtSum
+ << " direction:" << direction
+ << endl;
+ }
+ }
+
+
+ // Close debug streams
+ if (ownStr.valid())
+ {
+ ownStr.clear();
+ }
+ if (neiStr.valid())
+ {
+ neiStr.clear();
+ }
+
+
+
+ // Average the number of transformstions
+ nTransforms_ = (nTransforms_ + nTransformsOld)/2;
+
+ // Check that the match is complete
+ if (iter == maxIter_)
+ {
+ // The matching algorithm has exited without getting the
+ // srcSum and tgtSum above 1. This can happen because
+ // - of an incorrect setup
+ // - or because of non-exact meshes and truncation errors
+ // (transformation, accumulation of cutting errors)
+ // so for now this situation is flagged as a SeriousError instead of
+ // a FatalError since the default matchTolerance is quite strict
+ // (0.001) and can get triggered far into the simulation.
+ SeriousErrorIn
+ (
+ "void Foam::cyclicPeriodicAMIPolyPatch::resetPeriodicAMI"
+ "("
+ "const AMIPatchToPatchInterpolation::interpolationMethod&"
+ ") const"
+ )
+ << "Patches " << name() << " and " << neighbPatch().name()
+ << " do not couple to within a tolerance of "
+ << matchTolerance()
+ << " when transformed according to the periodic patch "
+ << periodicPatch.name() << "." << nl
+ << "The current sum of weights are for owner " << name()
+ << " : " << srcSum << " and for neighbour "
+ << neighbPatch().name() << " : " << tgtSum << nl
+ << "This is only acceptable during post-processing"
+ << "; not during running. Improve your mesh or increase"
+ << " the 'matchTolerance' setting in the patch specification."
+ << endl;
+ }
+
+ // Check that both patches have replicated an integer number of times
+ if
+ (
+ mag(srcSum - floor(srcSum + 0.5)) > srcSum*matchTolerance()
+ || mag(tgtSum - floor(tgtSum + 0.5)) > tgtSum*matchTolerance()
+ )
+ {
+ // This condition is currently enforced until there is more
+ // experience with the matching algorithm and numerics.
+ // This check means that e.g. different numbers of stator and
+ // rotor partitions are not allowed.
+ // Again see the section above about tolerances.
+ SeriousErrorIn
+ (
+ "void Foam::cyclicPeriodicAMIPolyPatch::resetPeriodicAMI"
+ "("
+ "const AMIPatchToPatchInterpolation::interpolationMethod&"
+ ") const"
+ )
+ << "Patches " << name() << " and " << neighbPatch().name()
+ << " do not overlap an integer number of times when transformed"
+ << " according to the periodic patch "
+ << periodicPatch.name() << "." << nl
+ << "The current matchTolerance : " << matchTolerance()
+ << ", sum of owner weights : " << srcSum
+ << ", sum of neighbour weights : " << tgtSum
+ << "." << nl
+ << "This is only acceptable during post-processing"
+ << "; not during running. Improve your mesh or increase"
+ << " the 'matchTolerance' setting in the patch specification."
+ << endl;
+ }
+
+ // Normalise the weights
+ AMIPtr_->normaliseWeights(true, false);
+ }
+}
+
+
+// * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * * //
+
+Foam::cyclicPeriodicAMIPolyPatch::cyclicPeriodicAMIPolyPatch
+(
+ const word& name,
+ const label size,
+ const label start,
+ const label index,
+ const polyBoundaryMesh& bm,
+ const word& patchType,
+ const transformType transform
+)
+:
+ cyclicAMIPolyPatch(name, size, start, index, bm, patchType, transform),
+ periodicPatchName_(word::null),
+ periodicPatchID_(-1),
+ nTransforms_(0),
+ nSectors_(0),
+ maxIter_(36)
+{}
+
+
+Foam::cyclicPeriodicAMIPolyPatch::cyclicPeriodicAMIPolyPatch
+(
+ const word& name,
+ const dictionary& dict,
+ const label index,
+ const polyBoundaryMesh& bm,
+ const word& patchType
+)
+:
+ cyclicAMIPolyPatch(name, dict, index, bm, patchType),
+ periodicPatchName_(dict.lookup("periodicPatch")),
+ periodicPatchID_(-1),
+ nTransforms_(dict.lookupOrDefault