View Issue Details
ID  Project  Category  View Status  Date Submitted  Last Update 

0003148  OpenFOAM  Bug  public  20190104 12:49  20190213 08:09 
Reporter  oertel59  Assigned To  henry  
Priority  normal  Severity  minor  Reproducibility  always 
Status  closed  Resolution  reopened  
Platform  GNU/Linux  OS  Ubuntu  OS Version  14.04 
Product Version  dev  
Fixed in Version  dev  
Summary  0003148: reactingMultiphaseEulerFoam: Dictionary entry order for the phases involved affects results  
Description  At least for reactingMultiphaseEulerFoam, it doesn't feel unreasonable to expect that the results are not affected by the order in which the phases involved are specified. I.e. in the phaseProperties dictionary, a change from phases (air water) to phases (water air) should yield identical results. This is only the case though if the simulation is incompressible, as evident from tests with a laminar bubbleColumn case using different entry orders and combinations of equationOfState. The original solver gives original/rhoConst/log.diff.alpha.airMean is empty original/rhoConstPerfectFluid/log.diff.alpha.airMean contains entries original/perfectGasRhoConst/log.diff.alpha.airMean contains entries original/perfectGasPerfectFluid/log.diff.alpha.airMean contains entries The attached patch separates the source term assembly and the solution of the alphaEqns. If only the changes to multiphaseSystem are considered, I get: multiphaseSystem/rhoConst/log.diff.alpha.airMean is empty multiphaseSystem/rhoConstPerfectFluid/log.diff.alpha.airMean is empty multiphaseSystem/perfectGasRhoConst/log.diff.alpha.airMean is empty multiphaseSystem/perfectGasPerfectFluid/log.diff.alpha.airMean contains entries Apparently, the order how the pEqnComps are added to the pEqn plays a role as well (roundoff behaviour?) If the pEqnComps are first summed up and then added to the pEqn as a whole, this yields: multiphaseSystemAndpEqn/rhoConst/log.diff.alpha.airMean is empty multiphaseSystemAndpEqn/rhoConstPerfectFluid/log.diff.alpha.airMean is empty multiphaseSystemAndpEqn/perfectGasRhoConst/log.diff.alpha.airMean is empty multiphaseSystemAndpEqn/perfectGasPerfectFluid/log.diff.alpha.airMean is empty. Admittedly, I am not so sure though whether the change to pEqn is convincing enough to be considered.  
Tags  No tags attached.  

0001reactingMultiphaseEulerFoamFixesabugwhereinput.patch (6,719 bytes)
From d00c51e3b73636fae3a85c9f81942098f6e07629 Mon Sep 17 00:00:00 2001 From: "Lehnigk, Ronald (FWDC)  10659" <r.lehnigk@hzdr.de> Date: Fri, 4 Jan 2019 12:54:34 +0100 Subject: [PATCH] reactingMultiphaseEulerFoam: Fixes a bug where input order affects results Prior to this commit, the order in which the phases involved are specified through the entry "phases" in the phaseProperties dictionary affected the results in case of an compressible simulation. Patch contributed by Institute of Fluid Dynamics, HelmholtzZentrum Dresden  Rossendorf (HZDR)  .../multiphaseSystem/multiphaseSystem.C  67 ++++++++++++ .../reactingMultiphaseEulerFoam/pU/pEqn.H  17 +++ 2 files changed, 55 insertions(+), 29 deletions() diff git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C index 62729c9b0..63f679db6 100644  a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C @@ 2,7 +2,7 @@ =========  \\ / F ield  OpenFOAM: The Open Source CFD Toolbox \\ / O peration  Website: https://openfoam.org  \\ / A nd  Copyright (C) 20132018 OpenFOAM Foundation + \\ / A nd  Copyright (C) 20132019 OpenFOAM Foundation \\/ M anipulation   License @@ 202,32 +202,38 @@ void Foam::multiphaseSystem::solveAlphas() } MULES::limitSum(alphafs, alphaPhiCorrs, fixedAlphaPhiCorrs);  // Solve for the moving phase alphas + // Compute alpha equation source terms for the moving phases + PtrList<volScalarField::Internal> Sps(movingPhases().size()); + PtrList<volScalarField::Internal> Sus(movingPhases().size()); forAll(movingPhases(), movingPhasei) { phaseModel& phase = movingPhases()[movingPhasei]; volScalarField& alpha = phase;  surfaceScalarField& alphaPhi = alphaPhiCorrs[phase.index()];  alphaPhi += upwind<scalar>(mesh_, phi_).flux(phase);  phase.correctInflowOutflow(alphaPhi);   volScalarField::Internal Sp + Sps.set (  IOobject + phase.index(), + new volScalarField (  "Sp",  mesh_.time().timeName(),  mesh_  ),  mesh_,  dimensionedScalar(dimless/dimTime, 0) + IOobject + ( + "Sp", + mesh_.time().timeName(), + mesh_ + ), + mesh_, + dimensionedScalar(dimless/dimTime, 0) + ) );  volScalarField::Internal Su + Sus.set (  "Su",  min(alpha, scalar(1))*fvc::div(fvc::absolute(phi_, phase.U())) + phase.index(), + new volScalarField + ( + "Su", + min(alpha, scalar(1))*fvc::div(fvc::absolute(phi_, phase.U())) + ) ); if (phase.divU().valid()) @@ 238,12 +244,12 @@ void Foam::multiphaseSystem::solveAlphas() { if (dgdt[celli] > 0.0) {  Sp[celli] = dgdt[celli];  Su[celli] += dgdt[celli]; + Sps[phase.index()][celli] = dgdt[celli]; + Sus[phase.index()][celli] += dgdt[celli]; } else if (dgdt[celli] < 0.0) {  Sp[celli] += + Sps[phase.index()][celli] += dgdt[celli] *(1  alpha[celli])/max(alpha[celli], 1e4); } @@ 265,29 +271,40 @@ void Foam::multiphaseSystem::solveAlphas() { if (dgdt2[celli] < 0.0) {  Sp[celli] += + Sps[phase.index()][celli] += dgdt2[celli] *(1  alpha2[celli])/max(alpha2[celli], 1e4);  Su[celli] = + Sus[phase.index()][celli] = dgdt2[celli] *alpha[celli]/max(alpha2[celli], 1e4); } else if (dgdt2[celli] > 0.0) {  Sp[celli] = dgdt2[celli]; + Sps[phase.index()][celli] = dgdt2[celli]; } } } } + } + + // Solve for the moving phase alphas + forAll(movingPhases(), movingPhasei) + { + phaseModel& phase = movingPhases()[movingPhasei]; + volScalarField& alpha = phase; + + surfaceScalarField& alphaPhi = alphaPhiCorrs[phase.index()]; + alphaPhi += upwind<scalar>(mesh_, phi_).flux(phase); + phase.correctInflowOutflow(alphaPhi); MULES::explicitSolve ( geometricOneField(), alpha, alphaPhi,  Sp,  Su + Sps[phase.index()], + Sus[phase.index()] ); phase.alphaPhiRef() = alphaPhi; diff git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H index f50de4b7b..a92bb2811 100644  a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H @@ 361,13 +361,22 @@ while (pimple.correct()) { fvScalarMatrix pEqn(pEqnIncomp);   forAll(phases, phasei) {  if (pEqnComps.set(phasei)) + fvScalarMatrix pEqnComp + ( + p_rgh, + pEqn.dimensions() + ); + + forAll(phases, phasei) {  pEqn += pEqnComps[phasei]; + if (pEqnComps.set(phasei)) + { + pEqnComp += pEqnComps[phasei]; + } } + + pEqn += pEqnComp; } pEqn.solve();  2.20.1 

Currently the phase fractions are solved in the order specified in the "phases" list and the latest solution of each is used so the order is likely to affect the solution unless an outer iteration is converged. It is possible to cache the intermediate results to ensure the order does not matter but it will be slightly slower (less convergent?), use more memory and add significant complexity, I guess this is what you have done in the patch you provide, but is it worth it? If you loop to convergence do you still get differences? How large are the difference? Do they matter if they are less than the numerical error of the solution? 

Pending feedback. 

I don't want to leave the question regarding the differences unanswered. I took the laminar bubble column tutorial again, with time averaging from 50100s. Perhaps this is not the best case. The differences are evident from the attached plot and vanish with the patch in place. The execution time is 1% higher with the patch. Using twice as many outer loops led to a crash of the case at about 5s physical time, so I couldn't produce the same plot again for the original solver with higher loop count. I am not sure whether this can be considered worth it and leave it up to your judgement and experience. 

If it is not possible to converge the case by increasing the number of outer loops then indeed the solution would be different by changing the solution order of any of the equations. Given that the case appears to be unstable this should be lookedat before substantial changes are made to the solution order based on the current results. 

Based on this case, I was unable to provide sufficient proof that this change is beneficial other than making the solution orderindependent at the cost of approximately 1% execution time and an unknown memory overhead. It doesn't converge properly  not for the original solver nor with the change. As a result, the differences between the solution for (air water) and (water air) with the original solver don't vanish either. Once I can come up with a case that provides the evidence, I will let you know. For now, you can close the issue. Thanks for the lesson. 

Suspended pending further analysis 
Date Modified  Username  Field  Change 

20190104 12:49  oertel59  New Issue  
20190104 12:49  oertel59  File Added: original.tar.gz  
20190104 12:49  oertel59  File Added: multiphaseSystemAndpEqn.tar.gz  
20190104 12:49  oertel59  File Added: multiphaseSystem.tar.gz  
20190104 12:49  oertel59  File Added: 0001reactingMultiphaseEulerFoamFixesabugwhereinput.patch  
20190104 13:56  henry  Note Added: 0010252  
20190123 11:10  henry  Assigned To  => henry 
20190123 11:10  henry  Status  new => closed 
20190123 11:10  henry  Resolution  open => suspended 
20190123 11:10  henry  Note Added: 0010261  
20190129 14:56  henry  Status  closed => feedback 
20190129 14:56  henry  Resolution  suspended => reopened 
20190129 16:00  oertel59  File Added: alpha.eps  
20190129 16:00  oertel59  Note Added: 0010269  
20190129 16:00  oertel59  Status  feedback => assigned 
20190129 16:19  henry  Note Added: 0010270  
20190213 08:03  oertel59  File Added: originalVsPatch.eps  
20190213 08:03  oertel59  File Added: original.eps  
20190213 08:03  oertel59  File Added: patch.eps  
20190213 08:03  oertel59  Note Added: 0010289  
20190213 08:09  henry  Status  assigned => closed 
20190213 08:09  henry  Fixed in Version  => dev 
20190213 08:09  henry  Note Added: 0010290 