| chrisb2244 |
March 18, 2014 23:29 |
Altered interDyMFoam hangs during MULES::explicitSolve
I'm running a slightly (although presumably critically, and badly) modified version of interDyMFoam in OF-2.2.2.
My createFields.H is modified to read
Code:
volScalarField alpha2
(
IOobject
(
"alpha2",
runTime.timeName(),
mesh
),
1.0 - alpha1
);
twoPhaseProperties.alpha1() = alpha1;
twoPhaseProperties.alpha2() = alpha2;
instead of
Code:
volScalarField& alpha1(twoPhaseProperties.alpha1());
volScalarField& alpha2(twoPhaseProperties.alpha2());
and volScalarField alpha1 is created using a library I have written.
Examining the alpha1 file written out by alpha1.write(), I can find no changes (apart from the values of alpha1) compared to a normally generated file.
When I run the solver, it stops (without error message) during the MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0); step in alphaEqn.H
I have modified it to read
Code:
Pout<< "Beginning explicitSolve(..)" << endl;
MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
Pout<< "explicitSolve(..) finished" << endl;
and the output to logs on running the solver reads
Code:
Execution time for mesh.update() = 0.06 s
time step continuity errors : sum local = 0, global = 0, cumulative = 0
GAMGPCG: Solving for pcorr, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 0, global = 0, cumulative = 0
hexRef4::setInstance(const fileName& inst) : Resetting file instance to "0.001"
Beginning explicitSolve(..)
and then ends.
Any suggestions as to what I should investigate here would be appreciated.
The code used to generate alpha1 is as follows:
Code:
wordList wantedTypes;
wantedTypes.append(cyclicPolyPatch::typeName);
wantedTypes.append(cyclicPolyPatch::typeName);
wantedTypes.append(zeroGradientFvPatchField<scalar>::typeName);
wantedTypes.append(zeroGradientFvPatchField<scalar>::typeName);
wantedTypes.append(emptyPolyPatch::typeName);
if (debug) Pout<< "wantedTypes = " << wantedTypes << endl;
volScalarField alpha1
(
IOobject
(
"alpha1",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::AUTO_WRITE,
true // registry
),
mesh_,
dimensionedScalar
(
"alpha1",
dimensionSet(0, 0, 0, 0, 0, 0, 0),
0.0
),
wantedTypes
);
if (debug) Pout<< "Field created " << endl;
double yHeight = yMid_; // This value is passed in at construction from the calling program
forAll(mesh_.cellCentres(), i)
{
double cellHeight = sqrt(mesh_.cellVolumes()[i] / cellDepth_);
yHeight = yMid_;
Foam::vector position(mesh_.cellCentres()[i]);
for (unsigned int k=0; k < cosineVector.size(); k++)
{
yHeight += cosineVector[k](position[0]);
}
double deltaY = yHeight - position[1];
if (deltaY < -cellHeight)
{
alpha1.internalField()[i] = 1.0;
}
else if (deltaY > cellHeight)
{
alpha1.internalField()[i] = 0.0;
}
else
{
double value_Alpha = (0.5+(-0.5 * (deltaY/cellHeight)));
alpha1.internalField()[i] = value_Alpha;
}
}
return alpha1;
cosineVector is a std::vector containing a series of functors, one for each chosen wavenumber, taking the value of x as an argument to the () operator.
The full alphaEqn.H file is below (but is identical to the one used by interDyMFoam)
Code:
{
word alphaScheme("div(phi,alpha)");
word alpharScheme("div(phirb,alpha)");
surfaceScalarField phic(mag(phi/mesh.magSf()));
phic = min(interface.cAlpha()*phic, max(phic));
surfaceScalarField phir(phic*interface.nHatf());
for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
{
surfaceScalarField phiAlpha
(
fvc::flux
(
phi,
alpha1,
alphaScheme
)
+ fvc::flux
(
-fvc::flux(-phir, alpha2, alpharScheme),
alpha1,
alpharScheme
)
);
Pout<< "Beginning explicitSolve(..)" << endl;
MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
Pout<< "explicitSolve(..) finished" << endl;
alpha2 = 1.0 - alpha1;
rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
}
Info<< "Phase-1 volume fraction = "
<< alpha1.weightedAverage(mesh.Vsc()).value()
<< " Min(alpha1) = " << min(alpha1).value()
<< " Max(alpha1) = " << max(alpha1).value()
<< endl;
}
|
