interIsoFoam runs into an endless loop

@Henning86 (44a84d47), any particular suggestion (if possible)? (or would you prefer us to take a stab at it?)

Many thanks both.

Will look into it

@HappyKiter i cannot reproduce the error on commit 9b0f1b67 on the develop branch

What version are you using?

I used 2012

I cannot reproduce it on ubuntu 18.04.

my g++ --version:

g++ (Ubuntu 7.5.0-3ubuntu1~18.04) 7.5.0

Could you provide:

Operating system
Hardware info
Compiler

I've downloaded the latest dev version. Once it's compiled I'll test the tutorial on this version and see if I get the same error as in v2012

System: Kernel: 4.12.14-lp151.28.87-default x86_64 bits: 64 Desktop: Gnome 3.26.2 Distro: openSUSE Leap 15.1

gcc version 7.5.0 (SUSE Linux)

What kind of hardware information do you require? The output of lscpu is:

` Architecture: x86_64 CPU op-mode(s): 32-bit, 64-bit Byte Order: Little Endian Address sizes: 48 bits physical, 48 bits virtual CPU(s): 2 On-line CPU(s) list: 0,1 Thread(s) per core: 1 Core(s) per socket: 2 Socket(s): 1 NUMA node(s): 1 Vendor ID: AuthenticAMD CPU family: 21 Model: 112 Model name: AMD E2-9010 RADEON R2, 4 COMPUTE CORES 2C+2G Stepping: 0 CPU MHz: 2000.000 CPU max MHz: 2000.0000 CPU min MHz: 1200.0000 BogoMIPS: 3992.18 Virtualization: AMD-V L1d cache: 32K L1i cache: 64K L2 cache: 1024K NUMA node0 CPU(s): 0,1 Flags: fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush mmx fxsr sse sse2 ht syscall nx mmxext fxsr_opt pdpe1gb rdtscp lm constant_tsc rep_good acc_power nopl nonstop_tsc cpuid extd_apicid aperfmperf pni pclmulqdq monitor ssse3 fma cx16 sse4_1 sse4_2 movbe popcnt aes xsave avx f16c lahf_lm cmp_legacy svm extapic cr8_legacy abm sse4a misalignsse 3dnowprefetch osvw ibs xop skinit wdt lwp fma4 tce nodeid_msr tbm perfctr_core perfctr_nb bpext ptsc mwaitx cpb hw_pstate ssbd vmmcall fsgsbase bmi1 avx2 smep bmi2 xsaveopt arat npt lbrv svm_lock nrip_save tsc_scale vmcb_clean flushbyasid decodeassists pausefilter pfthreshold avic v_vmsave_vmload vgif overflow_recov

`

I tried with the newest dev verion and the issue is still there. The latest output of the log is:

--> FOAM Warning : From void Foam::cutFaceAdvect::cutPoints(Foam::label, Foam::scalar, Foam::DynamicList<Foam::Vector >&) in file cellCuts/cutFace/cutFaceAdvect.C at line 928 cutPoints = 3((0.0045956218192765 0.00020064918188055 0.018) (0.0045956218192765 0.00020064918188055 0.0182) (0.0047954314635929 0.00020937305935361 0.018)) for pts = 4((0.0045956218192765 0.00020064918188055 0.018) (0.0045956218192765 0.00020064918188055 0.0182) (0.0047954314635929 0.00020937305935361 0.0182) (0.0047954314635929 0.00020937305935361 0.018)), f - f0 = 4(9254 9356 9357 9255) and f0 = 0 --> FOAM Warning : From void Foam::cutFaceAdvect::cutPoints(Foam::label, Foam::scalar, Foam::DynamicList<Foam::Vector >&) in file cellCuts/cutFace/cutFaceAdvect.C at line 928 cutPoints = 3((0.0043958121749602 -0.00019192530440748 0.0182) (0.0045956218192765 -0.00020064918188055 0.0182) (0.0043958121749602 -0.00019192530440748 0.0184)) for pts = 4((0.0043958121749602 -0.00019192530440748 0.0182) (0.0045956218192765 -0.00020064918188055 0.0182) (0.0045956218192765 -0.00020064918188055 0.0184) (0.0043958121749602 -0.00019192530440748 0.0184)), f - f0 = 4(9304 9305 9407 9406) and f0 = 0 --> FOAM Warning : From void Foam::cutFaceAdvect::cutPoints(Foam::label, Foam::scalar, Foam::DynamicList<Foam::Vector >&) in file cellCuts/cutFace/cutFaceAdvect.C at line 928 cutPoints = 3((0.0023977157317965 -0.00010468652967681 0.0194) (0.0025975253761128 -0.00011341040714987 0.0194) (0.0023977157317965 -0.00010468652967681 0.0196)) for pts = 4((0.0023977157317965 -0.00010468652967681 0.0194) (0.0025975253761128 -0.00011341040714987 0.0194) (0.0025975253761128 -0.00011341040714987 0.0196) (0.0023977157317965 -0.00010468652967681 0.0196)), f - f0 = 4(9906 9907 10009 10008) and f0 = 0

Maybe this is of some help

@Henning86 @kuti : The "f0 = 0" means that isoAdvector it is trying to calculate the 0-contour isoface on a face (which it should never do). It is presumably doing this for a face surrounded by cells with alpha = 0, so the contour catches all the 4 face points and therefore cannot decide how to cut the face (this is only unambiguous when there are two cut points along the perimeter of a face). Focus from here should have been on finding out why isoAdvector is used on such a face rather than trying to solve the subsequent problems in the bounding step.

I could not reproduce the results on another machine but i'm looking into.

let me try with osuse as well.

for the following configurations, the attached mwe has been successfully completed (although setAlphaField returns a Fatal Error [cannot find file 0/alpha.air]).

    base0    = base
    base1    = develop
    api      = 2012
    patch    = 210210
    HEAD     = 9b0f1b67df
    version  = com
    compiler = Gcc (system)
             = gcc (SUSE Linux) 7.4.1 20190905 [gcc-7-branch revision 275407]
    mpi      = SYSTEMOPENMPI
             = mpirun (Open MPI) 1.10.7.0.5e373bf1fd
    OS       = openSUSE Leap 15.1
    opts     = linux64GccDPInt32Opt

    base0    = base
    base1    = develop
    api      = 2012
    patch    = 210210
    HEAD     = 9b0f1b67df
    version  = com
    compiler = Clang (system)
             = clang version 9.0.1
    mpi      = SYSTEMOPENMPI
             = mpirun (Open MPI) 1.10.7.0.5e373bf1fd
    OS       = openSUSE Leap 15.1
    opts     = linux64ClangDPInt32Opt

Having said that, @Henning86, is there any chance for us to avoid any infinite-loop pitfalls in the while loop in question? To my understanding, there is no break mechanism thereat if things go awry.

That is should be possible.

My approach would be to get rid of the loop to see if this really is the problem and if it is the case i would refactor the residual part of plicRDF to improve readability and performance (maps are quite slow).

I've just downloaded clang (clang version 7.0.1 (tags/RELEASE_701/final 349238). I'm currently compiling the dev version with this compiler. Let's see if I get it compiled and if this is successful i can reproduce the error again

By the way I noticed that in file $FOAM_SRC/transportModels/geometricVoF/surfaceIterators/surfaceIteratorIso.C line 85

there is another while statement:

while (L2 - L1 > 1)

I'm not familiar with the code, but there may be also in this loop be a potential infinite loop. I'm I wrong?

patch.inifiteloop

isoAdvectionTemplates.C plicRDF.C

Could you apply the patch or replace the isoAdvectionTemplates.C and plicRDF.C in your installation?

I think there are two options where the code might fail in:

boundFlux
while (resIter.found())

i replaced the while loop in plicRDF and added a counter to boundFlux so we are able to test where the error is.

Does the error occur if you use the isoAlpha reconstructionScheme?

ok with this files the simulation runs til the end.

Super thanks

The problem seems to be the while (resIter.found()).

I will rewrite that part

plicRDF.C

Could you test this implementation and see if the infinite loop is still fixed?

Hei hei. I got the following compilation error:

reconstructionSchemes/plicSchemes/plicRDF/plicRDF.C:236:10: error: ‘normalRes’ was not declared in this scope List& normalResidual ^~~~~~~~~ reconstructionSchemes/plicSchemes/plicRDF/plicRDF.C:236:10: note: suggested alternative: ‘normal’ List& normalResidual ^~~~~~~~~ normal reconstructionSchemes/plicSchemes/plicRDF/plicRDF.C:236:19: error: template argument 1 is invalid List& normalResidual

are some files missing?

Version below is tested and produces exactly the same results as before but without the while loop

plicRDF.H plicRDF.C

ok now the pressure solver crashes at time step Time = 0.625957. But no infinite loop

Hello. I found another parameter settings for which interIsoFaom runs in an infinite loop. I attached the case. I run it with the dev version with the two file send by @Henning86 nOuter16nInner2send.tar.gz

@Henning86 do you run also in an infinite loop?

Final patch (Both while loops are replaced with for loops):

fixInfinteLoop.patch

@johan_roenby do you approve the patch?

I modified the case so the bubbles will leave the domain by using another outlet pressure boundary condition:

outlet
{
     type  prghTotalPressure;
     p0    uniform 0;
     value uniform 0;
}

During testing i realized that checkMesh fails and this is probably the reason why alpha gets unbounded. Nevertheless, the code should never go into an infinite loop.

nOuter2nInner2nAlphaBounds6sendUpdate.zip

For small scale multiphase simulations where surface tension plays an important role there is a time step restriction which is not implemented. For surface tension driven flows CFL numbers smaller than 0.05 are not uncommon.

Details can be found in (Popinet describes it also very accurately see literature ): https://arxiv.org/abs/2103.00870

@HappyKiter Could you test the patch with nOuter16nInner2send.tar.gz?

can you send me the files? I somehow have errors when applying the patch

plicRDF.H

plicRDF.C

isoAdvectionTemplates.C

it runs still in the same infinite loop

Do you get the same results if you exclude #include "centerOfGravityAir"?

The case passed with the files above

it runs still in the same infinite loop if i comment out #include "centerOfGravityAir"

surfaceIteratorPLIC.C

Could you apply this file and send me the output?

which output? The one to the screen? Should i send you the log.interIsoFaom file?

The infinte loop should be gone and the solver crashed (hopefully) with an error message FatalError: Did the solver crash and if yes what was the FatalError message

I got the following warning when compiling: sportModels/geometricVoF/surfaceIterators/surfaceIteratorPLIC.o surfaceIterators/surfaceIteratorPLIC.C: In member function ‘Foam::label Foam::surfaceIteratorPLIC::vofCutCell(Foam::label, Foam::scalar, Foam::scalar, Foam::label, Foam::vector)’: surfaceIterators/surfaceIteratorPLIC.C:113:36: warning: ‘a3’ may be used uninitialized in this function [-Wmaybe-uninitialized] << abort(FatalError); ^ surfaceIterators/surfaceIteratorPLIC.C:113:36: warning: ‘L3’ may be used uninitialized in this function [-Wmaybe-uninitialized]

is this ok?

This is the last output of the solver since a while:

isoAdvection: After conservative bounding: min(alpha) = -22.597932862311, max(alpha) = 1 + 6.2961928257186 Phase-1 volume fraction = 0.99669030874796 Min(alpha.water) = -22.597932862311 Max(alpha.water) = 7.2961928257186 GAMG: Solving for p_rgh, Initial residual = 0.85858469200651, Final residual = 0.018618733552631, No Iterations 1 time step continuity errors : sum local = 0.0036636819726502, global = -2.0431449927149e-06, cumulative = -1.4898521524037e-05 GAMG: Solving for p_rgh, Initial residual = 0.63009111336733, Final residual = 9.5582449747527e-10, No Iterations 15 time step continuity errors : sum local = 2.0962552262343e-10, global = -8.4935885194242e-11, cumulative = -1.4898606459923e-05 PIMPLE: iteration 11 isoAdvection: After conservative bounding: min(alpha) = -0.083578656012244, max(alpha) = 1 + 0.32142333754742 Phase-1 volume fraction = 0.99669048968762 Min(alpha.water) = -0.083578656012244 Max(alpha.water) = 1.3214233375474 GAMG: Solving for p_rgh, Initial residual = 0.76214216166438, Final residual = 0.006416428785324, No Iterations 1 time step continuity errors : sum local = 0.0017186062739367, global = 3.4259321968989e-06, cumulative = -1.1472674263024e-05 GAMG: Solving for p_rgh, Initial residual = 0.8076337915325, Final residual = 5.0907395992102e-10, No Iterations 16 time step continuity errors : sum local = 2.1793387596923e-10, global = 1.0681875665576e-10, cumulative = -1.1472567444267e-05 PIMPLE: iteration 12 isoAdvection: After conservative bounding: min(alpha) = -2.1229786522706, max(alpha) = 1 + 0.98204523786559 Phase-1 volume fraction = 0.99669104421056 Min(alpha.water) = -2.1229786522706 Max(alpha.water) = 1.9820452378656 GAMG: Solving for p_rgh, Initial residual = 0.79459058571972, Final residual = 0.013872241894828, No Iterations 1 time step continuity errors : sum local = 0.0073797994261729, global = -1.5264576177153e-05, cumulative = -2.673714362142e-05 GAMG: Solving for p_rgh, Initial residual = 0.7503769404026, Final residual = 7.8975989236977e-10, No Iterations 17 time step continuity errors : sum local = 5.6044182904409e-10, global = -3.1986004039718e-10, cumulative = -2.673746348146e-05 PIMPLE: iteration 13 isoAdvection: After conservative bounding: min(alpha) = -35.21770078879, max(alpha) = 1 + 64.798805448994 Phase-1 volume fraction = 0.99669949850926 Min(alpha.water) = -35.21770078879 Max(alpha.water) = 65.798805448994 GAMG: Solving for p_rgh, Initial residual = 0.93654377416058, Final residual = 0.0341077161499, No Iterations 6 time step continuity errors : sum local = 0.037264486733739, global = 0.0086252756565533, cumulative = 0.0085985381930719 GAMG: Solving for p_rgh, Initial residual = 0.67397821930766, Final residual = 9.2137890088027e-10, No Iterations 33 time step continuity errors : sum local = 5.6787986865813e-09, global = 1.8501678804116e-09, cumulative = 0.0085985400432397 PIMPLE: iteration 14

@Henning86 @kuti : Note that "min(alpha) = -22.597932862311, max(alpha) = 1 + 6.2961928257186", i.e. the alpha field is already severely unbounded and there is no chance that the heuristic bounding step can deal with such severe unboundedness. Here the root to the boundedness should be sought earlier in the simulation.

no error message. Maybe can you send me the content of the whole folder geometricVof of you and a past it into mine just to be sure we have the same files

geometricVoF.zip

Hello @Henning86. With the latest files interIsoFoam crashes with the following message:

isoAdvection: After conservative bounding: min(alpha) = -563.23336830065, max(alpha) = 1 + 1044.0950152637 Phase-1 volume fraction = 0.99676141485956 Min(alpha.water) = -563.23336830065 Max(alpha.water) = 1045.0950152637 #0 Foam::error::printStack(Foam::Ostream&) at ??:? #1 (closed) Foam::sigFpe::sigHandler(int) at ??:? #2 (closed) ? in /lib64/libpthread.so.0 #3 (closed) Foam::GAMGSolver::scale(Foam::Field&, Foam::Field&, Foam::lduMatrix const&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, Foam::Field const&, unsigned char) const at ??:? #4 (closed) Foam::GAMGSolver::Vcycle(Foam::PtrListFoam::lduMatrix::smoother const&, Foam::Field&, Foam::Field const&, Foam::Field&, Foam::Field&, Foam::Field&, Foam::Field&, Foam::Field&, Foam::PtrList<Foam::Field >&, Foam::PtrList<Foam::Field >&, unsigned char) const at ??:? #5 (closed) Foam::GAMGSolver::solve(Foam::Field&, Foam::Field const&, unsigned char) const at ??:? #6 (closed) Foam::fvMatrix::solveSegregated(Foam::dictionary const&) at ??:? #7 (closed) Foam::fvMatrix::solveSegregatedOrCoupled(Foam::dictionary const&) at ??:? #8 (closed) Foam::fvMesh::solve(Foam::fvMatrix&, Foam::dictionary const&) const at ??:? #9 (closed) ? at ??:? #10 (closed) __libc_start_main in /lib64/libc.so.6 #11 (closed) ? at /home/abuild/rpmbuild/BUILD/glibc-2.26/csu/../sysdeps/x86_64/start.S:122

I had to recompile the solver too not only the library to get to this point. Strange.

the first case I sent crashes also with a floating point exception

Thats great news as the crash can be explained with the mesh quality and the inifite loop seems to be gone!

Hi @Henning86 - can you create a branch/merge request for your proposed fixes?

Sure.

First i want to confirm that the results of the fix and the previous results are identical. I hope to get it done tomorrow or Friday

mentioned in commit 497215cc

mentioned in merge request !436 (merged)

mentioned in commit 7cd3b270

closed

mentioned in merge request !455 (merged)

interIsoFoam runs into an endless loop

Child items

Activity

Admin message

interIsoFoam runs into an endless loop

Child items

Linked items

Related merge requests

Activity