MFiX Solver Crash: Invalid operation in wrap_pow

MFiX
#1

Hello,

I am attempting to model pyrolysis of a particle, and I am getting a solver crash on the first time step.
image

xpow.c leads me to this line:
image
I am unsure how to fix the error since xpow.c seems to be a “read only” function that is called.

I had previously got the program to work on mfix-22.4.3, but it stopped working. I updated to 23.1, and still couldn’t find a solution.

I have uploaded the .zip file below.
park_model_2023-04-21T135255.795423.zip (80.9 KB)

#2

Thanks for the bug report.

You are correct that xpow is not intended to be modified - it is an accelerated version of the pow function, which makes MFiX simulations run faster. However, even when disabling it I still get an invalid operation. xpow itself is not the problem.

I modified the code in xpow.c as follows:

double __wrap_pow(double x, double y) {
    double r, s, x0, y0;

    fprintf(stderr, "%f %f\n", x, y);
    fflush(stderr);
...

that is it simply logs its arguments to standard error on every call. And I see that the last call has

-1.223616 0.200000

that is, negative value to a fractional power which is not permitted.

This is coming from calc_gamma_des.f:142

139	         CASE (GUNN)
140	            Epg=Ep_g(IJK)
141	            N_Nu = (7.0d0-10.0*epg+5.0*epg*epg)*(1.0d0+0.7d0 *((N_Re)**0.2*(N_Pr)**THIRD)) &
142	                   + (1.33d0-2.40d0*epg+1.20d0*epg*epg)*(N_Re)**0.7*(N_Pr)**THIRD

That’s a big mess but I see the exponent 0.2 appears in the expression (N_Re)**0.2 so it appears that N_Re < 0. And the debugger (or print statements) confirms this:

(gdb) print N_re
$2 = -1.2236160342884503

So the question is - how did we get a negative Reynolds number? (It should always be > 0)

#3

Thanks for the response.

The Reynolds number should be calculated from calc_gamma_des.f:142
image

DP is the characteristic length, SLIP is the particle velocity calculated earlier in file, RO_g is the density, and mu is the viscosity.

I think RO_g and Mu_g are just grabbing the fluid properties and probably aren’t a source of problems.
At a glance, SLIP doesn’t seem like it can return negative values since it is calculating a magnitude and the operations there should handle any negative values.
image

DP should be calculated from line 121, and similarly should have positive values:
image

A “simple” solution might be to take the absolute value of the reynolds number before doing operations with it, but that solution seems like it could have the potential to hide other issues in the simulation.

#4

Since we don’t know what the source of the problem is we must be careful not to make assumptions. More print statements.

! Calculate the particle Reynolds Number
      write(*,*) "IJK", ijk, "MU", mu_g(ijk)
      IF(MU_G(IJK) > ZERO) THEN
           write(*,*)  "DP", dp, "SLIP", slip, "RO", ro_g(ijk)
           N_Re = (DP*SLIP*RO_g(IJK)) / MU_g(IJK)
      ELSE
           N_Re = LARGE_NUMBER
      ENDIF

And in the output we see:

 IJK 170                        MU 1.8000000000000000E-005
 DP  1.6000000000000000E-002    SLIP 1.4521806455795890E-004 
 RO -9.4793167968789156

So we’ve got a negative fluid density in cell 170.

I agree with your comment that simply taking an absolute value is not a good solution because it will merely hide other problems. We have to figure out where the negative ro (ρ) is coming from.

(Also I wish the original authors had spelled it rho but that’s too much code to change now!!)

– Charles

#5

Searching through the code (grep -ir 'negative density') I found many places where there are checks for negative density.

And searching the interface I find this control:

shot-2023-04-22_06-44-13
shot-2023-04-22_06-44-131190×743 163 KB

Enabling it results in a file ROgErr.log which contains:

  One or more cells have reported a negative/NaN gas density (RO_g(IJK)). If
  this is a persistent issue, lower keyword UR_FAC(1).

  Simulation time:  0.10000E-04

    IJK:      170       I:    4  J:    2  K:    4
      RO_g:  -9.4793      P_g: -0.72207E+06  T_g:   293.15

Unfortunately, lowering the underrelaxation factor UR_FAC(1) does not help unless I bring it all the way down to zero, in which case the run diverges.

#6

A few comments/suggestions:

  1. Please review your initial condition. Pressure is set to 10135 Pa (about 0.1 atmosphere).
  2. There is only a wall BC (no pressure outlet).
  3. Please try to turn off the reactions (or set the reaction rate to zero). If the reaction rate is too large, it can generate/consume a large amount of gas and make the solution unphysical.
1 Like
#7

Good catch. That almost certainly should be 101325 Pa.

#8

Thanks for the help. I was able to get the model to work. I added a pressure outlet boundary condition along the top of the domain. I also adjusted the particle input file to generate the particle higher up. I’m not sure if the particle can clip through a fixed wall boundary when it is generated by the input file and cause a crash.

Thank you again for all the assistance.