[chrisb2244]

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

volScalarField alpha2
(
	IOobject
	(
		"alpha2",
		runTime.timeName(),
		mesh
	),
	1.0 - alpha1
);
	
twoPhaseProperties.alpha1() = alpha1;
twoPhaseProperties.alpha2() = alpha2;

instead of

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

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

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:

        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)

{
    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;
}

[chrisb2244]

Looks like my problem lay in the way I was creating the alpha1 field within the solver.

By copying all of the code lines from the regenerateAlpha() function (above) into the solver, so that it reads

wordList wantedTypes;
wantedTypes.append(cyclicPolyPatch::typeName);
wantedTypes.append(cyclicPolyPatch::typeName);
wantedTypes.append(zeroGradientFvPatchField<scalar>::typeName);
wantedTypes.append(zeroGradientFvPatchField<scalar>::typeName);
wantedTypes.append(emptyPolyPatch::typeName);

volScalarField alpha1
(
	IOobject
	(
		"alpha1",
		runTime.timeName(),
		mesh,
		IOobject::NO_READ,
		IOobject::AUTO_WRITE,
		true // registry
	),
	mesh,
	dimensionedScalar
	(
		"alpha1",
		dimensionSet(0, 0, 0, 0, 0, 0, 0),
		0.0
	),
	wantedTypes
);

alpha1 = regAClass.regenerateAlpha();
alpha1.write();

the solver now works as expected