Issue in Handling Unsmooth Outputs After AMR

[MontelukastNa]

Dear FoamScience and kind friend M. Elwardi Fadeli,
Thank you for your great work on this AMR and Loadbalance library.

Recently, I’ve encountered an issue that may merit further clarification or improvement,

Specifically, there may be sometime that when the adjustableRunTime writeInterval occurs immediately after a refineInterval (e.g., writeInterval 0.001s, refineInterval 30 timesteps, unrefineInterval 90 timesteps), the output fields sometimes appear unphysical or unsmooth (e.g., P result as follows Fig).
I suspect this might be due to the fact that the solution fields require several time steps to fully adjust and stabilize after mesh refinement and rebalancing.

It’s unclear if this is the exact cause, but the results imply that output written immediately after refinement might be unsmooth because the mesh and interpolation have just changed.

At present, I have some preliminary thoughts regarding how output might be better handled after mesh refinement in AMR simulations, and I would greatly appreciate any guidance:

• One idea is to set OpenFOAM’s writeInterval to None, then instead manage output through a controlDict function triggered after each refineInterval, which would output by delaying several time steps (e.g., about 10) to allow the fields to smooth. However, this makes it challenging to maintain outputs at fixed adjustableRunTime intervals (e.g., every 0.001 seconds).

• Another idea is to keep writeInterval based on adjustableRunTime, but introduce a condition so that when refineInterval is close to the output time, smaller time steps are used for several steps to smooth the fields before writing output.

I do not have experience implementing such approaches myself, so any suggestions, best practices, or insights on addressing this issue would be most welcome. 🥂Thank you very much for your time and the continued development of this great tool.

Best regards,
Na

[FoamScience]

Hello; This seems like a more serious issue (more than an "output time thing"). Looks like you need:

1. More nBufferLayers in dynamicMeshDict so that you avoid losing accuracy around the interface you're refining.
2. Use correctFluxes keyword in dyanmicMeshDict (something like correctFluxes ((phi none) (... any more special flux fields));). You might experiment with correcting the flux with the velocity field (instead of none) as this depends on the solver you are running.

Now, let me explain few notes;

I suspect this might be due to the fact that the solution fields require several time steps to fully adjust and stabilize after mesh refinement and rebalancing.
Generally, we cannot run a CFD simulation for even a single step with a degrading field - it will mess up the accuracy from that point onward; so, this is not allowed to happen

As the "unrefinement" (coarsening) happens, face fluxes generally need to be reconstructed from velocity (because that's what we solved the PDE for) as this is the best option to recover the fluxes. Some solvers do this already (so you don't need the correctFluxes for these).

Try these two changes and let me know what comes out. Also, might be good to show the mesh and maybe dynamicMeshDict content :)

[MontelukastNa]

Dear FoamScience and kind friend M. Elwardi Fadeli,

Thanks for your kind reply.

You're right, many thanks for your clarification. Althrough the unsmooth solution fields still "exist" even after running reconstructPar, the bubble behavior shows good agreement with the experimental data reported in the reference. I also apologize for my earlier doubt. As you mentioned, even a single degraded field in a CFD simulation can compromise the accuracy from that point onward.

Below, I present both the original and modified dynamicMeshDict settings, along with the observed AMR test results. For comparison, SHM results are also included.

I have observed the same phenomenon with both blastAMR and SHM: either P or U shows smooth results (e.g. smooth P with unsmooth U, or vice versa), but not both simultaneously. So I believe the unsmooth outcome may be due to solver settings that are not appropriate, as you mentioned.

For the correctFluxes, i don't know whether the interFlow solvers do this already, so i try these two changes (nBufferLayers and correctFluxes), as you suggested.

My enviroument info:
OpenFOAM: v2312
Solver: interFlow (TwoPhaseFlow v2312)
BlastAMR: V2212 

💡1️⃣ Original dynamicMeshDict without correctFluxes keyword, and increasing nBufferLayers from 3 to 4 appears to have no further effect on the solution fields, the unsmooth regions still persist.

/·························original dynamicMeshDict ·························/
dynamicFvMesh   adaptiveFvMesh;
refiner         polyRefiner;
errorEstimator  delta;
deltaField      alpha.water;
balance         yes;
refineInterval  30;
unrefineInterval 90;
maxRefinement   4;
nBufferLayers   3;  // 3 or 4;
nRefinementBufferLayers 3;
nUnrefinementBufferLayers 6;
lowerRefineLevel 0.01;
upperRefineLevel 0.99;
unrefineLevel   0.01;
upperUnrefineLevel 0.99;
dumpLevel       no;
refineProbes    no;

💡2️⃣Then I use correctFluxes keyword in dynamicMeshDict, but it seems can't work with the interFlow solver.

/·························correctFluxes dynamicMeshDict ·························/
dynamicFvMesh   adaptiveFvMesh;
refiner         polyRefiner;
errorEstimator  delta;
deltaField      alpha.water;
balance         yes;
refineInterval  30;
unrefineInterval 90;
maxRefinement   4;
nBufferLayers   3;
nRefinementBufferLayers 3;
nUnrefinementBufferLayers 6;
lowerRefineLevel 0.01;
upperRefineLevel 0.99;
unrefineLevel   0.01;
upperUnrefineLevel 0.99;
dumpLevel       no;
refineProbes    no;
correctFluxes   ( ( phi U ) );  // adding to my simulation

It runs for several time steps, but crashes before reaching the step of calculating the maximum imbalance. Increasing nBufferLayers form 3 to 4 get the same error.

/························error output ( correctFluxes dynamicMeshDict )·························/
PIMPLE: iteration 1
Protecting 0 cells next to requested boundary patches.
Selected 8 cells for refinement out of 115540.
[10] [stack trace]
[10] =============
[11] [stack trace]
[11] =============
[12] [stack trace]

💡3️⃣Furthermore, I use mapSurfaceFields keyword in dynamicMeshDict, but its outcomes same as the case No. 1️⃣ in this issue (Original dynamicMeshDict).

/·························mapSurfaceFields  dynamicMeshDict ·························/
dynamicFvMesh   adaptiveFvMesh;
refiner         polyRefiner;
errorEstimator  delta;
deltaField      alpha.water;
balance         yes;
refineInterval  30;
unrefineInterval 90;
maxRefinement   4;
nBufferLayers   3;
nRefinementBufferLayers 3;
nUnrefinementBufferLayers 6;
lowerRefineLevel 0.01;
upperRefineLevel 0.99;
unrefineLevel   0.01;
upperUnrefineLevel 0.99;
dumpLevel       no;
refineProbes    no;
mapSurfaceFields  ( Uf  Uf_0 ); // adding to my simulation

💡4️⃣ Lastly, the SHM results is given for comparison.

❗❗❗However, whether using SHM or AMR, I’ve encountered a common issue: discontinuities appear in the solution fields at locations where the mesh refinement level changes.
From my comparison, I noticed that either p is smooth and U is not, or vice versa — but never both simultaneously.
✅Since this behavior also occurs in the SHM case, where AMR is not involved at all. Thus the issue may not be related to AMR itself, but rather to certain aspects of the simulate settings or the interFlow solver.

Below are the contents of my fvSchemes and fvSolution files for your reference.
I’d appreciate any suggestions you might have regarding possible adjustments or optimizations.

⏩info of fvScheme

/·························fvScheme info·························/
ddtSchemes
{
    default         CrankNicolson 0.9;
}
 
gradSchemes
{
    default         pointCellsLeastSquares;//Gauss linear;//
}
 
divSchemes
{
    default             none;
    div(rhoPhi,U)	Gauss limitedLinearV 1;
    div(phi,alpha)	Gauss vanLeer;
    div(phirb,alpha)	Gauss interfaceCompression;
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}
 
laplacianSchemes
{
    default	Gauss linear corrected;
}
 
interpolationSchemes
{
    default	linear;
}
 
snGradSchemes
{
    default	corrected;
}
 
fluxRequired
{
    default         no;
    p_rgh;
    pcorr;
    alpha.water;
}

⏩info of fvSolution

/·························fvSolution info·························/
solvers
{
    "alpha.water.*"
    {
        advectionScheme isoAdvection;//geoAdvection;//
        reconstructionScheme plicRDF; //RDFadvect isoInverseDistance

        // advectionScheme MULESScheme;//geoAdvection;//
        // reconstructionScheme isoAlpha; //plicRDF; //RDFadvect isoInverseDistance

        // isoAdvector parameters
        vof2IsoTol      1e-6;
        surfCellTol     1e-6;
        nAlphaBounds    3;
        snapTol         1e-12;
        clip            true;
        gradAlphaNormal false;
        writeVTK        false;

        // MULES parameters
        nAlphaCorr      2;
        cAlpha          1;
        icAlpha         0;
        MULESCorr       yes;
        nLimiterIter    3;
        alphaApplyPrevCorr yes;
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-8;
        relTol          0;
        minIter         1;

        // parameters read by both isoAdvector and MULES
        nAlphaSubCycles 1;
    }

    pcorr
    {
        solver          GAMG;
        smoother        DIC;
        tolerance       1e-10;
        relTol          0;
    }

    pcorrFinal
    {
        $pcorr;
        relTol          0;
    }

    p_rgh
    {
        solver          GAMG;
        smoother        DIC;
        tolerance       1e-07;
        relTol          0.001;
    }

    p_rghFinal
    {
        $p_rgh;
        tolerance       1e-8; /* Should be 1 to 2 orders of 
                                 magnitude smaller than 
                                 alpha.phase1.*.surfCellTol */
        relTol          0;
    }

    "(U|k|epsilon).*"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-07;
        relTol          0;
    }
}

PIMPLE
{
    momentumPredictor   yes;
    nOuterCorrectors    1;
    nCorrectors         3;
    nNonOrthogonalCorrectors 1;
    correctPhi true;
    pRefCell            0;
    pRefValue           0;
}

relaxationFactors
{
    equations
    {
        ".*"            1;
    }
}

I would greatly appreciate any suggestions, best practices, or insights you may have on how to deal with it.
🥂Thank you once again for your time and for the ongoing development of such a valuable AMR tool.

Best regards,
Na

[FoamScience]

Thanks for the extensive report; but this, most likely, is not a library issue as you discovered. what you have to do:

• Refine less often: unrefineInterval #calc "5*$refineInterval"; (you may want to make it even larger)
• Have a decent amount of nBuffers: nBufferLayers 8; - to 10 (This is valid for both blastAMR and SHM)
  This should fix the issues you are seeing. interFlow should correct the fluxes on its own.

Also mapSurfaceFields has no effect; it's never read by the library. Guess I need better docs of possible keywords...

[MontelukastNa]

Dear FoamScience and kind friend M. Elwardi Fadeli,

Thanks for your kind suggestions.

I’ve followed your advice and modified the settings:

unrefineInterval 150;  // originally set via calc "5*$refineInterval"
nBufferLayers 10;      // increased from the previous value

However, after re-running the simulation twice with these AMR settings, I didn’t observe any noticeable changes. The mesh size remained the same, suggesting that nBufferLayers might not be taking effect. Unfortunately, the discontinuity in the U field still persists, whereas the pressure field remains smooth, as previously summarized.

I'm using the interFlow solver in OpenFOAM v2312, and I’ve confirmed that the dictionary changes are being correctly loaded (for example, I even conducted a third test where I changed nBufferLayers from 10 to 12 during runtime, and the log file showed that the update was applied).

For now, this is the conclusion based on careful double-checks and repeated testing (I plan to rerun it twice more shortly). I’ll continue investigating and will keep you posted if I find anything new.

Any further insights and suggestions would be greatly appreciated.

Best regards,
Na

[FoamScience]

You may also try using none instead of "U" in correctFluxes to tell the AMR engine not to interpolate fluxes on new faces, but I don't think It'll make a big difference.

It's strange that nBufferLayers appears inactive; will investigate it when i get time.

Instead, we can fix this by including those cells into the refinement; however for this I need the case as only a coded refinement can take care of this. If you are going through with this, please:

1. Fork this repository
2. Clone locally, create a new branch, add your case to tutorials folder, commit the changes and push to your fork
3. Open a Pull Request into v2112 branch. Please also mark the PR as "Draft"

Make sure to only commit the clean state of the case, no result files

We can easily collaborate on this from there; as I feel this is a general issue many users will face...

[MontelukastNa]

Dear FoamScience and kind friend M. Elwardi Fadeli,

Thank you very much for your kind reply and for offering to help, I truly appreciate it.

Apologies for the slight delay in getting back to you; I had a fever on one of the days recently and couldn’t check messages in time.

I’d be glad to follow your instructions — forking the repository, adding the case to the tutorials folder, and opening a draft PR to the v2112 branch with a clean setup. I’ll get started on that as soon as possible.

Before I begin, I’d just like to mention that I’ve been using Henning’s (twoPhaseFlow-v2312+) interFlow solver, not the ones that come from OpenFOAM v2312. I wonder if that might make things more complicated or inconvenient for you?

If that’s not an issue, I’ll proceed right away.

Thanks again for your generous support!

Best regards,
Na

[FoamScience]

No need to apologize for anything, this is the deal with Open Source projects, things are expected to go a little slow as not everyone will have time to do everything at once...

Effects of working with interFlow

1. Your case won't be included in this repo since it needs an external dependency
2. The fix I have in mind has nothing to do with the solver; only the way we select cells for refinement is important so I should be able to help without actually running the case.
3. You will be doing the actual testing (I don't expect this to be much work)

If you're fine with this, feel free to add your case

[MontelukastNa]

Dear FoamScience and kind friend M. Elwardi Fadeli,

Thank you very much for your time and helpful suggestions regarding my previous question about the unsmooth outputs.

I'm using the interFlow solver in OpenFOAM v2312, and I’ve confirmed that the dictionary changes are being correctly loaded (for example, I even conducted a third test where I changed nBufferLayers from 10 to 12 during runtime, and the log file showed that the update was applied).

📑 Fig. 1. A bubble formation under an additional body force

After further investigation, I found that one source of the issue was related to the way I explicitly introduced an additional body force into the simulation. This implementation caused discontinuities in the velocity field U, which confounded the diagnosis and amplified the apparent non-smoothness near refinement-level transitions.
I sincerely apologize for the confusion caused by this.

However, I should also mention that even in cases without any additional body force, such as the example shown below (a bubble rising in a quiescent liquid driven only by gravity), a certain level of non-smoothness in the solution still seems to exist near refinement level transitions. So the problem may not be completely understood yet. If I obtain any further insights or updates, I will report them here.

📑 Fig. 2. A bubble rising in a quiescent liquid driven only by gravity (without any additional body force)

Thank you again for your patience and help.

Best regards,
Na