Gao's Drag bug solved

MFiX Bug report
#1

Hi,
With the mfix-21.1.0 update, I decided to test the new implemented drag models. When I tested the Gao drag model it gave me error in the fluid_bed_tfm_2d and fluid_bed_tfm_3d tutorials. To solve this bug, I propose to modify drag_gs.f: 1755
h3d = hInf * (ONE - exp(-alphaCoeff*(vSlipStar-VZero)**powerVal))
by the following

Pot = -alphaCoeff*(vSlipStar-VZero)**powerVal
IF (Pot>=73.6827d0)THEN
exp_temp = LARGE_NUMBER
ELSEIF (Pot<=-73.6827d0)THEN
exp_temp = -LARGE_NUMBER
ELSE
exp_temp = EXP(Pot) 
ENDIF
h3d =  hInf * (ONE - exp_temp)

Also, to add local variables at the top:
DOUBLE PRECISION :: Pot, exp_temp
This prevents the numbers from exceeding the fortran 1e32 limit.

I await your comments to solve this bug.
Best regards

#2

Thanks for testing the drag. Can you provide me more information on how the error happened? And can you also print out the values of some related variables when the error shows up: ep_g, alphaCoeff, vSlipStar, VZero, powerVal?

#3

Hi and thank you very much for answering. I want to print the values of some related variables so I can help development, I am trying (WRITE or PRINT) but still can’t get the values to show up in the stdout or mfix message. How could I do that?

I literally opened both tutorials (fluid_bed_tfm_2d or fluid_bed_tfm_3d), just changed the drag model and run with default solver.

#4

Copy the file drag_gs.f to your project folder;

Add the following lines to the subroutine of the Gao drag

if (Pot>=73.6827d0) then
WRITE(,’(/,“values of related variables”)’)
WRITE(,’(4(f18.8))’) alphaCoeff, vSlipStar, VZero, powerVal
endif

Recompile, run with only one core, then the values of related variables will show on your screen.

#5

I managed to extract data from the first iteration of the fluid_bed_tfm_2d tutorial and it came out 420 results where Pot> = 73.6827d0, out of a total of a 20x100 mesh.

image
It can be said that alphaCoeff is the cause of the divergence

#6

Thanks JPGD. I run without GUI, they are fine, both 2D and 3D run a long time without show any issue.

While When I run with GUI, indeed we reproduced the issue you found.

Jeff pointed out the reason:

The issue comes from a floating point exceptions (fpe), here specifically an overflow. The GUI solver from the GUI traps fpe. You can build the batch solver to track fpes with:

cmake … -DENABLE_MPI=0 -DCMAKE_Fortran_COMPILER=gfortran -DCMAKE_Fortran_FLAGS="-O0 -fcheck=all -ffpe-trap=invalid,zero,overflow -g"

And it should reproduce the error with the batch solver. The overflow occurs in the calculation of h3d in the drag law. The overflow occurs if the argument becomes larger than approximately 700.

We will remedy this issue based on the suggestion by you in the next release. Thanks.