FPE at des_thermo_cond_mod.f:295, incomplete stack trace

Please select the most relevant MFiX category: | Installation | How to | Bug report | Share | for this topic.
Dear developers, hello. At present, I am simulating biomass pyrolysis (simply setting a working condition and adding a chemical reaction), but I have encountered some problems (in the Solids interface, turning on the Thermal Cond. function and running it for 2 steps will result in an error message: Solver crash! The MFix solvent has terminated unexpectedly. Error information: float in cal operation in (null) at (null): 0; If the thermal conduction function is not turned on, the calculation example can run normally. I cannot solve this problem.
I will appreciate it if this problem can be solved.
delete_2023-08-31T153552.188859.zip (68.4 MB)

Uploading: 2aca8ec8385a2fe01317d650e7d5b4f.png…

Hi. There are a few things going on here that I don’t completely understand.
First off, the stack trace says (null):0 where there should be a filename - when I run this case I get the following backtrace:

Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.

Backtrace for this error:
#0 radius
  at des/des_thermo_cond_mod.f:295
#1 des_thermo_cond_MOD_des_conduction
  at des/des_thermo_cond_mod.f:158
#2 calc_force_dem_mod_MOD_calc_force_dem
  at des/calc_force_dem.f:206
#3 des_time_march_MOD_des_time_step
  at des/des_time_march.f:194
#4 run_dem
  at mfix.f:211
#5 run_mfix
  at mfix.f:146
#6 main_MOD_run_mfix0
  at main.f:84

I’m not sure why you are not getting the full error report. You can try building the solver in Debug mode.

Here’s the offending code.

284	! Calculate
285	      VALUE = (R1**2 - R2**2 + CENTER_DIST_CORR**2)/(2.0d0 * R1 * CENTER_DIST_CORR)
286	! Check to ensure that VALUE is less than or equal to one. If VALUE is
287	! greater than one, the triangle inequality has been violated. Therefore
288	! there is no intersection between the fluid lens surrounding the larger
289	! particle and the surface of the smaller particle.
290	! Thus, there is no particle-fluid-particle heat transfer.
291	      IF( VALUE .GT. 1.0d0) THEN
292	         RADIUS = -1.0d0
293	      ELSE
294	! Calculate beta (Law of cosines)
295	         BETA = ACOS( VALUE )

The error is at line 295. The variable value is equal to -2.908 so it’s impossible to take the arccos of that.

Now, how did value get negative?

(gdb) print center_dist_corr 

This causes the denominator of the fraction at line 285 to be negative.

126	! Correct particle overlap to account for artificial softening
127	            OLAP_Sim = MAX_RAD + MIN_RAD - CENTER_DIST
129	               IF(Im.ge.Jm)THEN
130	                  OLAP_Actual = CORRECT_OLAP(OLAP_Sim,Jm,Im)
131	               ELSE
132	                  OLAP_Actual = CORRECT_OLAP(OLAP_Sim,Im,Jm)
133	               ENDIF
134	               CENTER_DIST_CORR = MAX_RAD + MIN_RAD - OLAP_Actual
(gdb) p max_rad
(gdb) p min_rad
(gdb) p center_dist

So center_dist is (slightly) less than max_rad + min_rad and we are executing the “Correct particle overlap” code, which is doing something wrong …

(gdb) p olap_sim
(gdb) p olap_actual
(gdb) p max_rad + min_rad - olap_actual

There should probably be some checks in the code to make sure this distance never gets negative.

First, make sure you use the grid based neighbor search, this will be faster than N-square, but this is not related to your issue.

The Young’s moduli you are using for the conduction are too small. You should use values that are representative of the actual material, probably in the order of 10 GPa to 100 GPa. This is different from the Youngs moduli you may use with the Hertzain collision model which are typically smaller that the actual material property (otherwise the DEM time step will likely be prohibitively small).

Thank you very much for your answer and feedback. This problem has been resolved.
Best wishes.

Thank you very much for your answer and feedback. This problem has been resolved.
Best wishes.

Dear cgw, excuse me, how did you obtain the error code? Could you please tell me the steps to obtain it? Recently, I have been trying to find a way to obtain error codes, but no matter how I operate, I cannot do it. If you have the time to provide detailed steps, I would really appreciate it.
Best wishes!

I’m running on Linux and using gdb to examine core dumps.

[quote=“GL1, post:10, topic:4829, full:true”]
Dear cgw, excuse me. Could you please provide specific debugging screenshots or processes? I really don’t understand this operation, and no one around me can operate this process. I really don’t have a solution now. If this problem can be solved, I would greatly appreciate it.
Best wishes!

Dear CGW, excuse me. According to the debugging results you provided, I tried on the Linux system, but it cannot be fully replicated. Due to my limited abilities, I really can’t come up with any solutions. The main question is why I can’t obtain the values shown in the figure below.

My own operation process is shown in the following figure:First, run Mfixsolver, and then type (gdb) p max_ Rad(), but there are no results or values. May I ask why this is? Or rather, how did you operate it? I am truly powerless and eager to learn this knowledge and content. Thank you again. If this issue can be resolved, I would greatly appreciate it.

Also, I have a question. Is the error code shown in the following figure obtained in the Mfix software or in the terminal interface?

Finally, thanks for your help again.
Best wishes.

Line 291 should be changed to

       IF(VALUE .LT. -1.0d0 .or. VALUE .GT. 1.0d0) THEN

I’ll follow up with some notes about debugging with GDB

Thank you very much for your reply. I am very eager to learn about gdb debugging, and I kindly request that you follow up and update as soon as possible. Thank you again.
Best wishes!

Hi. First off, It’s still not clear to me why the GUI is showing null:0 instead of the correct line number. I’d like to get to the bottom of that.

As promised, a few notes on GDB.

I suggest you try to find a GDB tutorial, I can offer some tips but I can’t teach you everything about GDB. The commands you need most are “where”, “up”, “down”, “list”, “print” and “whatis”. You can type “help” in GDB for the built-in help system which is not bad but presumes a bit of background knowledge.

Also note that debugging generally works better on Linux. On Linux you can do post-mortem debugging via “core” files, which don’t exist on Windows. On Windows you have to start the mfixsolver inside gdb, which it looks like you have figured out already.

The reason p min_rad is not working in your gdb session example is that you are on the wrong stack frame. If you type where at the (gdb) prompt it will show you the stack. Variables are only accessible in the stack frame where they live. In this example, the stack looks something like this:

(gdb) where
#0  0x00007f8a6861f3ae in feraiseexcept () from /lib64/libm.so.6
#1  0x00007f8a68621750 in acosf64 () from /lib64/libm.so.6
#2  0x00007f8a185d71cd in des_conduction::radius (r2=0.0020000000000000009, 
    r1=<optimized out>)
    at /home/cgw/mambaforge/envs/mfix-23.2/share/mfix/src/model/des/des_thermo_cond_mod.f:295
#3  0x00007f8a185d7a9f in des_thermo_cond::des_conduction (i=4, j=473, 
    center_dist=0.0039993857351258745, im=<optimized out>, 
    iijk=<optimized out>)
    at /home/cgw/mambaforge/envs/mfix-23.2/share/mfix/src/model/des/des_thermo_cond_mod.f:158

so the bottom of the stack is the library function feraiseexcept which is responsible for floating point exceptions. That function does not know about min_rad, etc. To see the variables which are part of des_thermo_conduction.f we need to go up 2 levels:

(gdb) up
#1  0x00007f8a68621750 in acosf64 () from /lib64/libm.so.6
(gdb) up
#2  0x00007f8a185d71cd in des_conduction::radius (r2=0.0020000000000000009, r1=<optimized out>) at /home/cgw/mambaforge/envs/mfix-23.2/share/mfix/src/model/des/des_thermo_cond_mod.f:295
295	         BETA = ACOS( VALUE )

Now we have a new problem because the term value seems to have a special meaning which blocks us accessing the variable:

(gdb) p value
Attempt to use a type name as an expression
(gdb) whatis value
type = Type, C_Union :: value   # ???

So now we look at more code to see how value is computed:

(gdb) list 270
265	!......................................................................!
266	! Subroutine Procedure: FUNCTION RADIUS                                !
267	!                                                                      !
268	! Propose: Calculate the center line radius between two particles.     !
269	! Depending on what is given as input, this function calculates:       !
270	!   1) radius of contact area between two particles                    !
271	!   2) radius delineating particle-fluid-particle region.              !
272	!......................................................................!
(gdb)   # hit return for more
277	! Radius values
279	! Distance between particle centers
281	! Generic value
284	! Calculate
(gdb)  # hit return for more
285	      VALUE = (R1**2 - R2**2 + CENTER_DIST_CORR**2)/(2.0d0 * R1 * CENTER_DIST_CORR)

Ok, here’s the computation of value. GDB can evaluate expressions, so let’s try to reproduce the computation:

(gdb) p (R1**2 - R2**2 + CENTER_DIST_CORR**2)/(2.0d0 * R1 * CENTER_DIST_CORR)
No symbol "CENTER_DIST_CORR" in current context.

Now we’re seeing a new issue, which is the fact that the compiler optimizer pass will eliminate some intermediate variables thereby making them inaccessible. So we rebuild the solver in Debug mode which turns the optimizer off:

Now, we’re in business!

(gdb) p (R1**2 - R2**2 + CENTER_DIST_CORR**2)/(2.0d0 * R1 * CENTER_DIST_CORR)
$2 = -2.9083823068238357

Hope this helps,

– Charles

Thank you for your reply. Your reply has been very helpful to me.
Best wishes!

1 Like