ENH: Propagated caching of HbyA across solvers

This commit is contained in:
andy
2012-03-02 18:15:54 +00:00
parent 658f0a0680
commit 912a20b7a3
28 changed files with 321 additions and 259 deletions

View File

@ -2,19 +2,24 @@
volScalarField rAU("rAU", 1.0/UEqn.A());
surfaceScalarField rAUf("(1|A(U))", fvc::interpolate(rAU));
U = rAU*UEqn.H();
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
phi = (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rAU, U, phi);
surfaceScalarField phig(rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
surfaceScalarField buoyancyPhi(rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
phi -= buoyancyPhi;
surfaceScalarField phiHbyA
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
+ fvc::ddtPhiCorr(rAU, U, phi)
- phig
);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@ -24,14 +29,14 @@
if (pimple.finalNonOrthogonalIter())
{
// Calculate the conservative fluxes
phi -= p_rghEqn.flux();
phi = phiHbyA - p_rghEqn.flux();
// Explicitly relax pressure for momentum corrector
p_rgh.relax();
// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
U -= rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rAUf);
U = HbyA - rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
U.correctBoundaryConditions();
}
}

View File

@ -2,20 +2,27 @@
volScalarField rAU("rAU", 1.0/UEqn().A());
surfaceScalarField rAUf("(1|A(U))", fvc::interpolate(rAU));
U = rAU*UEqn().H();
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn().H();
UEqn.clear();
phi = fvc::interpolate(U) & mesh.Sf();
adjustPhi(phi, U, p_rgh);
surfaceScalarField phig(rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
surfaceScalarField buoyancyPhi(rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
phi -= buoyancyPhi;
surfaceScalarField phiHbyA
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
);
adjustPhi(phiHbyA, U, p_rgh);
phiHbyA -= phig;
while (simple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@ -25,14 +32,14 @@
if (simple.finalNonOrthogonalIter())
{
// Calculate the conservative fluxes
phi -= p_rghEqn.flux();
phi = phiHbyA - p_rghEqn.flux();
// Explicitly relax pressure for momentum corrector
p_rgh.relax();
// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
U -= rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rAUf);
U = HbyA - rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
U.correctBoundaryConditions();
}
}

View File

@ -8,21 +8,26 @@
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
U = rAU*UEqn.H();
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
phi = fvc::interpolate(rho)*
surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
surfaceScalarField phiHbyA
(
(fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rAU, rho, U, phi)
"phiHbyA",
fvc::interpolate(rho)
*(
(fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rAU, rho, U, phi)
)
+ phig
);
surfaceScalarField buoyancyPhi(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
phi += buoyancyPhi;
fvScalarMatrix p_rghDDtEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+ fvc::div(phi)
+ fvc::div(phiHbyA)
);
while (pimple.correctNonOrthogonal())
@ -38,14 +43,14 @@
if (pimple.finalNonOrthogonalIter())
{
// Calculate the conservative fluxes
phi += p_rghEqn.flux();
phi = phiHbyA + p_rghEqn.flux();
// Explicitly relax pressure for momentum corrector
p_rgh.relax();
// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
U += rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorAUf);
U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
U.correctBoundaryConditions();
K = 0.5*magSqr(U);
}

View File

@ -5,20 +5,27 @@
volScalarField rAU(1.0/UEqn().A());
surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
U = rAU*UEqn().H();
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn().H();
UEqn.clear();
phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
bool closedVolume = adjustPhi(phi, U, p_rgh);
surfaceScalarField phig(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
surfaceScalarField buoyancyPhi(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
phi -= buoyancyPhi;
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf())
);
bool closedVolume = adjustPhi(phiHbyA, U, p_rgh);
phiHbyA -= phig
while (simple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phi)
fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phiHbyA)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@ -27,14 +34,14 @@
if (simple.finalNonOrthogonalIter())
{
// Calculate the conservative fluxes
phi -= p_rghEqn.flux();
phi = phiHbyA - p_rghEqn.flux();
// Explicitly relax pressure for momentum corrector
p_rgh.relax();
// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
U -= rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorAUf);
U = HbyA - rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
U.correctBoundaryConditions();
}
}

View File

@ -8,24 +8,27 @@
volScalarField rAU(1.0/UEqn().A());
surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
U = rAU*UEqn().H();
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn().H();
surfaceScalarField phiU
surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::interpolate(rho)
*(
(fvc::interpolate(U) & mesh.Sf())
(fvc::interpolate(HbyA) & mesh.Sf())
+ fvc::ddtPhiCorr(rAU, rho, U, phi)
)
+ phig
);
phi = phiU - rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
{
fvScalarMatrix p_rghDDtEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+ fvc::div(phi)
+ fvc::div(phiHbyA)
);
// Thermodynamic density needs to be updated by psi*d(p) after the
@ -57,7 +60,11 @@
if (nonOrth == nNonOrthCorr)
{
phi += p_rghEqn.flux();
phi = phiHbyA + p_rghEqn.flux();
U = HbyA
+ rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
U.correctBoundaryConditions();
K = 0.5*magSqr(U);
}
}
@ -65,9 +72,6 @@
thermo.rho() += psi*p_rgh;
}
// Correct velocity field
U += rAU*fvc::reconstruct((phi - phiU)/rhorAUf);
U.correctBoundaryConditions();
p = p_rgh + rho*gh;
// Update pressure time derivative