00001 if (pressureImplicitPorosity) 00002 { 00003 U = trTU()&UEqn().H(); 00004 } 00005 else 00006 { 00007 U = trAU()*UEqn().H(); 00008 } 00009 00010 UEqn.clear(); 00011 00012 phi = fvc::interpolate(rho*U) & mesh.Sf(); 00013 bool closedVolume = adjustPhi(phi, U, p); 00014 00015 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) 00016 { 00017 tmp<fvScalarMatrix> tpEqn; 00018 00019 if (pressureImplicitPorosity) 00020 { 00021 tpEqn = (fvm::laplacian(rho*trTU(), p) == fvc::div(phi)); 00022 } 00023 else 00024 { 00025 tpEqn = (fvm::laplacian(rho*trAU(), p) == fvc::div(phi)); 00026 } 00027 00028 tpEqn().setReference(pRefCell, pRefValue); 00029 // retain the residual from the first iteration 00030 if (nonOrth == 0) 00031 { 00032 eqnResidual = tpEqn().solve().initialResidual(); 00033 maxResidual = max(eqnResidual, maxResidual); 00034 } 00035 else 00036 { 00037 tpEqn().solve(); 00038 } 00039 00040 if (nonOrth == nNonOrthCorr) 00041 { 00042 phi -= tpEqn().flux(); 00043 } 00044 } 00045 00046 #include <finiteVolume/continuityErrs.H> 00047 00048 // Explicitly relax pressure for momentum corrector 00049 p.relax(); 00050 00051 if (pressureImplicitPorosity) 00052 { 00053 U -= trTU()&fvc::grad(p); 00054 } 00055 else 00056 { 00057 U -= trAU()*fvc::grad(p); 00058 } 00059 00060 U.correctBoundaryConditions(); 00061 00062 bound(p, pMin); 00063 00064 // For closed-volume cases adjust the pressure and density levels 00065 // to obey overall mass continuity 00066 if (closedVolume) 00067 { 00068 p += (initialMass - fvc::domainIntegrate(psi*p)) 00069 /fvc::domainIntegrate(psi); 00070 } 00071 00072 rho = thermo.rho(); 00073 rho.relax(); 00074 Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl; 00075 00076 // ************************ vim: set sw=4 sts=4 et: ************************ //