ENH: pimpleFoam, rhoPimpleFoam, interDyMFoam: Rationalized mesh-motion support

Added support for mesh-motion update within PIMPLE loop in pimpleFoam and rhoPimpleFoam.
This commit is contained in:
Henry Weller 2017-11-30 13:07:42 +00:00 committed by Andrew Heather
parent 889791060c
commit e50af751a4
37 changed files with 240 additions and 325 deletions

View File

@ -1,5 +1,3 @@
InfoInFunction << "consistent" << endl;
if (!pimple.SIMPLErho())
{
rho = thermo.rho();

View File

@ -1,10 +0,0 @@
#include "readTimeControls.H"
bool correctPhi = pimple.dict().lookupOrDefault<Switch>
(
"correctPhi",
!isA<staticFvMesh>(mesh)
);
bool checkMeshCourantNo =
pimple.dict().lookupOrDefault<Switch>("checkMeshCourantNo", false);

View File

@ -57,12 +57,11 @@ int main(int argc, char *argv[])
#include "setRootCase.H"
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "createControl.H"
#include "createDyMControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "createFieldRefs.H"
#include "createRhoUfIfPresent.H"
#include "createTimeControls.H"
turbulence->validate();
@ -78,81 +77,81 @@ int main(int argc, char *argv[])
while (runTime.run())
{
#include "readControls.H"
#include "readDyMControls.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
autoPtr<volScalarField> divrhoU;
if (correctPhi)
{
// 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
autoPtr<volScalarField> divrhoU;
if (correctPhi)
{
divrhoU.reset
divrhoU.reset
(
new volScalarField
(
new volScalarField
(
"divrhoU",
fvc::div(fvc::absolute(phi, rho, U))
)
);
}
if (LTS)
{
#include "setRDeltaT.H"
}
else
{
#include "compressibleCourantNo.H"
#include "setDeltaT.H"
}
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
// Store momentum to set rhoUf for introduced faces.
autoPtr<volVectorField> rhoU;
if (rhoUf.valid())
{
rhoU.reset(new volVectorField("rhoU", rho*U));
}
// Do any mesh changes
mesh.update();
#include "updateRhoUf.H"
if (mesh.changing())
{
MRF.update();
if (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 (checkMeshCourantNo)
{
#include "meshCourantNo.H"
}
}
"divrhoU",
fvc::div(fvc::absolute(phi, rho, U))
)
);
}
if (pimple.nCorrPIMPLE() <= 1)
if (LTS)
{
#include "rhoEqn.H"
#include "setRDeltaT.H"
}
else
{
#include "compressibleCourantNo.H"
#include "setDeltaT.H"
}
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
if (pimple.firstIter() || moveMeshOuterCorrectors)
{
// Store momentum to set rhoUf for introduced faces.
autoPtr<volVectorField> rhoU;
if (rhoUf.valid())
{
rhoU.reset(new volVectorField("rhoU", rho*U));
}
// Do any mesh changes
mesh.update();
if (mesh.changing())
{
MRF.update();
if (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 (checkMeshCourantNo)
{
#include "meshCourantNo.H"
}
}
}
if (pimple.nCorrPIMPLE() <= 1)
{
#include "rhoEqn.H"
}
#include "UEqn.H"
#include "EEqn.H"

View File

@ -50,8 +50,7 @@ int main(int argc, char *argv[])
#include "setRootCase.H"
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "createControl.H"
#include "createTimeControls.H"
#include "createDyMControls.H"
#include "createFields.H"
#include "createFieldRefs.H"
#include "createRhoUf.H"
@ -67,7 +66,7 @@ int main(int argc, char *argv[])
while (runTime.run())
{
#include "readControls.H"
#include "readDyMControls.H"
{
// Store divrhoU from the previous mesh so that it can be mapped

View File

@ -1,12 +0,0 @@
#include "createControl.H"
#include "createTimeControls.H"
bool correctPhi
(
pimple.dict().lookupOrDefault("correctPhi", false)
);
bool checkMeshCourantNo
(
pimple.dict().lookupOrDefault("checkMeshCourantNo", false)
);

View File

@ -85,7 +85,7 @@ int main(int argc, char *argv[])
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "initContinuityErrs.H"
#include "createControls.H"
#include "createDyMControls.H"
#include "createFields.H"
#include "createUfIfPresent.H"
#include "CourantNo.H"
@ -99,7 +99,7 @@ int main(int argc, char *argv[])
while (runTime.run())
{
#include "readControls.H"
#include "readDyMControls.H"
#include "CourantNo.H"
#include "setDeltaT.H"
@ -107,34 +107,36 @@ int main(int argc, char *argv[])
Info<< "Time = " << runTime.timeName() << nl << endl;
mesh.update();
#include "updateUf.H"
if (mesh.changing())
{
MRF.update();
if (correctPhi)
{
// Calculate absolute flux from the mapped surface velocity
phi = mesh.Sf() & Uf();
#include "correctPhi.H"
// Make the flux relative to the mesh motion
fvc::makeRelative(phi, U);
}
if (checkMeshCourantNo)
{
#include "meshCourantNo.H"
}
}
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
if (pimple.firstIter() || moveMeshOuterCorrectors)
{
mesh.update();
if (mesh.changing())
{
MRF.update();
if (correctPhi)
{
// Calculate absolute flux
// from the mapped surface velocity
phi = mesh.Sf() & Uf();
#include "correctPhi.H"
// Make the flux relative to the mesh motion
fvc::makeRelative(phi, U);
}
if (checkMeshCourantNo)
{
#include "meshCourantNo.H"
}
}
}
#include "UEqn.H"
// --- Pressure corrector loop

View File

@ -1,5 +0,0 @@
#include "readTimeControls.H"
correctPhi = pimple.dict().lookupOrDefault("correctPhi", false);
checkMeshCourantNo = pimple.dict().lookupOrDefault("checkMeshCourantNo", false);

View File

@ -60,7 +60,7 @@ int main(int argc, char *argv[])
#include "setRootCase.H"
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "createControls.H"
#include "createDyMControls.H"
#include "createFields.H"
#include "createUcf.H"
#include "initContinuityErrs.H"
@ -69,7 +69,7 @@ int main(int argc, char *argv[])
while (runTime.run())
{
#include "readControls.H"
#include "readDyMControls.H"
#include "CourantNo.H"
#include "setDeltaT.H"

View File

@ -1,12 +0,0 @@
#include "createControl.H"
#include "createTimeControls.H"
bool correctPhi
(
pimple.dict().lookupOrDefault("correctPhi", false)
);
bool checkMeshCourantNo
(
pimple.dict().lookupOrDefault("checkMeshCourantNo", false)
);

View File

@ -1,5 +0,0 @@
#include "readTimeControls.H"
correctPhi = pimple.dict().lookupOrDefault("correctPhi", false);
checkMeshCourantNo = pimple.dict().lookupOrDefault("checkMeshCourantNo", false);

View File

@ -53,8 +53,7 @@ int main(int argc, char *argv[])
#include "setRootCase.H"
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "createControl.H"
#include "createTimeControls.H"
#include "createDyMControls.H"
#include "createFields.H"
#include "createFieldRefs.H"
#include "createRhoUf.H"
@ -70,7 +69,7 @@ int main(int argc, char *argv[])
while (runTime.run())
{
#include "readControls.H"
#include "readDyMControls.H"
{
// Store divrhoU from the previous time-step/mesh for the correctPhi

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2012-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -53,7 +53,6 @@ int main(int argc, char *argv[])
#include "setRootCase.H"
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "createControl.H"
#include "createControls.H"
#include "createFields.H"
#include "createUf.H"

View File

@ -1,9 +1,4 @@
#include "createTimeControls.H"
bool correctPhi
(
pimple.dict().lookupOrDefault("correctPhi", true)
);
#include "createDyMControls.H"
scalar maxAcousticCo
(

View File

@ -1,4 +1,4 @@
#include "readTimeControls.H"
#include "readDyMControls.H"
correctPhi = pimple.dict().lookupOrDefault("correctPhi", true);
maxAcousticCo = readScalar(runTime.controlDict().lookup("maxAcousticCo"));

View File

@ -60,10 +60,9 @@ int main(int argc, char *argv[])
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "initContinuityErrs.H"
#include "createControl.H"
#include "createDyMControls.H"
#include "createFields.H"
#include "createUf.H"
#include "createControls.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
@ -77,7 +76,7 @@ int main(int argc, char *argv[])
while (runTime.run())
{
#include "readControls.H"
#include "readDyMControls.H"
// Store divU from the previous mesh so that it can be mapped
// and used in correctPhi to ensure the corrected phi has the

View File

@ -1,9 +0,0 @@
#include "readTimeControls.H"
correctPhi = pimple.dict().lookupOrDefault("correctPhi", true);
checkMeshCourantNo =
pimple.dict().lookupOrDefault("checkMeshCourantNo", false);
moveMeshOuterCorrectors =
pimple.dict().lookupOrDefault("moveMeshOuterCorrectors", false);

View File

@ -3,7 +3,7 @@ CorrectPhi
U,
phi,
p_rgh,
surfaceScalarField("rAUf", fvc::interpolate(rAU)),
surfaceScalarField("rAUf", fvc::interpolate(rAU())),
geometricZeroField(),
pimple
);

View File

@ -1,14 +0,0 @@
bool correctPhi
(
pimple.dict().lookupOrDefault("correctPhi", true)
);
bool checkMeshCourantNo
(
pimple.dict().lookupOrDefault("checkMeshCourantNo", false)
);
bool moveMeshOuterCorrectors
(
pimple.dict().lookupOrDefault("moveMeshOuterCorrectors", false)
);

View File

@ -59,28 +59,31 @@ int main(int argc, char *argv[])
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "initContinuityErrs.H"
#include "createControl.H"
#include "createTimeControls.H"
#include "createDyMControls.H"
#include "createFields.H"
#include "createAlphaFluxes.H"
volScalarField rAU
tmp<volScalarField> rAU;
if (correctPhi)
(
IOobject
rAU = new volScalarField
(
"rAU",
runTime.timeName(),
IOobject
(
"rAU",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1.0)
dimensionedScalar("rAU", dimTime/dimDensity, 1)
)
);
#include "correctPhi.H"
#include "createUf.H"
#include "createUfIfPresent.H"
turbulence->validate();
@ -95,7 +98,7 @@ int main(int argc, char *argv[])
while (runTime.run())
{
#include "readControls.H"
#include "readDyMControls.H"
if (LTS)
{
@ -117,16 +120,10 @@ int main(int argc, char *argv[])
{
if (pimple.firstIter() || moveMeshOuterCorrectors)
{
scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
mesh.update();
if (mesh.changing())
{
Info<< "Execution time for mesh.update() = "
<< runTime.elapsedCpuTime() - timeBeforeMeshUpdate
<< " s" << endl;
// Do not apply previous time-step mesh compression flux
// if the mesh topology changed
if (mesh.topoChanging())

View File

@ -1,12 +1,12 @@
{
rAU = 1.0/UEqn.A();
surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU()));
volVectorField HbyA(constrainHbyA(rAU()*UEqn.H(), U, p_rgh));
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::flux(HbyA)
+ fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, Uf)
+ fvc::interpolate(rho*rAU())*fvc::ddtCorr(U, phi, Uf)
);
MRF.makeRelative(phiHbyA);
@ -47,7 +47,7 @@
p_rgh.relax();
U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
U = HbyA + rAU()*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
}
@ -55,11 +55,8 @@
#include "continuityErrs.H"
{
Uf = fvc::interpolate(U);
surfaceVectorField n(mesh.Sf()/mesh.magSf());
Uf += n*(phi/mesh.magSf() - (n & Uf));
}
// Correct Uf if the mesh is moving
fvc::correctUf(Uf, U, phi);
// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi, U);
@ -76,4 +73,9 @@
);
p_rgh = p - rho*gh;
}
if (!correctPhi)
{
rAU.clear();
}
}

View File

@ -1,9 +0,0 @@
#include "readTimeControls.H"
correctPhi = pimple.dict().lookupOrDefault("correctPhi", true);
checkMeshCourantNo =
pimple.dict().lookupOrDefault("checkMeshCourantNo", false);
moveMeshOuterCorrectors =
pimple.dict().lookupOrDefault("moveMeshOuterCorrectors", false);

View File

@ -64,9 +64,7 @@ int main(int argc, char *argv[])
#include "setRootCase.H"
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "createControl.H"
#include "createTimeControls.H"
#include "../interFoam/interDyMFoam/createDyMControls.H"
#include "createDyMControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
@ -96,7 +94,7 @@ int main(int argc, char *argv[])
while (runTime.run())
{
#include "../interFoam/interDyMFoam/readControls.H"
#include "readDyMControls.H"
// Store divU from the previous mesh so that it can be mapped
// and used in correctPhi to ensure the corrected phi has the

View File

@ -54,27 +54,31 @@ int main(int argc, char *argv[])
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "initContinuityErrs.H"
#include "createControl.H"
#include "createTimeControls.H"
#include "createDyMControls.H"
#include "createFields.H"
volScalarField rAU
tmp<volScalarField> rAU;
if (correctPhi)
(
IOobject
rAU = new volScalarField
(
"rAU",
runTime.timeName(),
IOobject
(
"rAU",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1.0)
dimensionedScalar("rAU", dimTime/dimDensity, 1)
)
);
#include "correctPhi.H"
#include "createUf.H"
#include "createUfIfPresent.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
@ -88,7 +92,7 @@ int main(int argc, char *argv[])
while (runTime.run())
{
#include "readControls.H"
#include "readDyMControls.H"
#include "CourantNo.H"
#include "alphaCourantNo.H"
@ -120,7 +124,7 @@ int main(int argc, char *argv[])
if (mesh.changing() && correctPhi)
{
// Calculate absolute flux from the mapped surface velocity
phi = mesh.Sf() & Uf;
phi = mesh.Sf() & Uf();
#include "correctPhi.H"

View File

@ -59,8 +59,6 @@ int main(int argc, char *argv[])
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "initContinuityErrs.H"
#include "createControl.H"
#include "createTimeControls.H"
#include "createDyMControls.H"
#include "createFields.H"
@ -89,7 +87,7 @@ int main(int argc, char *argv[])
while (runTime.run())
{
#include "readControls.H"
#include "readDyMControls.H"
#include "CourantNo.H"
#include "setDeltaT.H"

View File

@ -497,6 +497,12 @@ public:
// Mesh motion
//- Is mesh dynamic
virtual bool dynamic() const
{
return false;
}
//- Is mesh moving
bool moving() const
{

View File

@ -125,6 +125,12 @@ public:
// Member Functions
//- Is mesh dynamic
virtual bool dynamic() const
{
return true;
}
//- Update the mesh for both mesh motion and topology change
virtual bool update() = 0;
};
@ -135,10 +141,6 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "staticFvMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -23,7 +23,7 @@ License
\*---------------------------------------------------------------------------*/
#include "dynamicFvMesh.H"
#include "staticFvMesh.H"
// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //

View File

@ -1,8 +1,9 @@
#include "createControl.H"
#include "createTimeControls.H"
bool correctPhi
(
pimple.dict().lookupOrDefault("correctPhi", true)
pimple.dict().lookupOrDefault("correctPhi", mesh.dynamic())
);
bool checkMeshCourantNo

View File

@ -0,0 +1,19 @@
#include "readTimeControls.H"
correctPhi = pimple.dict().lookupOrDefault
(
"correctPhi",
correctPhi
);
checkMeshCourantNo = pimple.dict().lookupOrDefault
(
"checkMeshCourantNo",
checkMeshCourantNo
);
moveMeshOuterCorrectors = pimple.dict().lookupOrDefault
(
"moveMeshOuterCorrectors",
moveMeshOuterCorrectors
);

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -77,6 +77,12 @@ public:
// Member Functions
//- Is mesh dynamic
virtual bool dynamic() const
{
return false;
}
//- Dummy update function which does not change the mesh
virtual bool update();
};

View File

@ -25,7 +25,7 @@ Global
createRhoUf
Description
Creates and initialises the velocity field rhoUf if present.
Creates and initialises the velocity field rhoUf if required.
\*---------------------------------------------------------------------------*/
@ -33,19 +33,25 @@ Description
autoPtr<surfaceVectorField> rhoUf;
IOobject rhoUfHeader
(
"rhoUf",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
);
if (rhoUfHeader.typeHeaderOk<surfaceVectorField>(true))
if (mesh.dynamic())
{
Info<< "Reading face momentum rhoUf\n" << endl;
rhoUf.reset(new surfaceVectorField(rhoUfHeader, mesh));
Info<< "Constructing face momentum rhoUf" << endl;
rhoUf.reset
(
new surfaceVectorField
(
IOobject
(
"rhoUf",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
fvc::interpolate(rho*U)
)
);
}
// ************************************************************************* //

View File

@ -31,28 +31,6 @@ Description
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
if (mesh.changing())
{
if (!rhoUf.valid())
{
Info<< "Constructing face momentum rhoUf" << endl;
rhoUf.reset
(
new surfaceVectorField
(
IOobject
(
"rhoUf",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
fvc::interpolate(rho*U)
)
);
}
}
#error Remove rhoUpdateUf.H
// ************************************************************************* //

View File

@ -25,7 +25,7 @@ Global
createUf
Description
Creates and initialises the velocity field Uf if present.
Creates and initialises the velocity field Uf if required.
\*---------------------------------------------------------------------------*/
@ -33,19 +33,25 @@ Description
autoPtr<surfaceVectorField> Uf;
IOobject UfHeader
(
"Uf",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
);
if (UfHeader.typeHeaderOk<surfaceVectorField>(true))
if (mesh.dynamic())
{
Info<< "Reading face velocity Uf\n" << endl;
Uf.reset(new surfaceVectorField(UfHeader, mesh));
Info<< "Constructing face velocity Uf\n" << endl;
Uf.reset
(
new surfaceVectorField
(
IOobject
(
"Uf",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
fvc::interpolate(U)
)
);
}
// ************************************************************************* //

View File

@ -31,28 +31,6 @@ Description
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
if (mesh.changing())
{
if (!Uf.valid())
{
Info<< "Constructing face velocity Uf" << endl;
Uf.reset
(
new surfaceVectorField
(
IOobject
(
"Uf",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
fvc::interpolate(U)
)
);
}
}
#error Remove updateUf.H
// ************************************************************************* //

View File

@ -215,7 +215,7 @@ ddtCorr
const autoPtr<GeometricField<Type, fvsPatchField, surfaceMesh>>& Uf
)
{
if (U.mesh().changing())
if (U.mesh().dynamic())
{
return ddtCorr(U, Uf());
}
@ -280,7 +280,7 @@ ddtCorr
const autoPtr<GeometricField<Type, fvsPatchField, surfaceMesh>>& Uf
)
{
if (U.mesh().changing())
if (U.mesh().dynamic())
{
return ddtCorr(rho, U, Uf());
}

View File

@ -228,7 +228,7 @@ void Foam::fvc::correctUf
{
const fvMesh& mesh = U.mesh();
if (mesh.changing())
if (mesh.dynamic())
{
Uf() = fvc::interpolate(U);
surfaceVectorField n(mesh.Sf()/mesh.magSf());
@ -247,7 +247,7 @@ void Foam::fvc::correctRhoUf
{
const fvMesh& mesh = U.mesh();
if (mesh.changing())
if (mesh.dynamic())
{
rhoUf() = fvc::interpolate(rho*U);
surfaceVectorField n(mesh.Sf()/mesh.magSf());

View File

@ -47,7 +47,7 @@ correctFluxes
(phi none)
(nHatf none)
(rhoPhi none)
(alphaPhi none)
(alphaPhi0.water none)
(ghf none)
);