Hello,everyone!
I use the SQP method to simulate irregular particles and write particle_input.data file imports particle location information, but it cannot be successfully read and there is no particle display. Do you know where the mistake is? Thank you！
n2d-lsfb_2023-04-21T220448.184259.zip (581.3 KB)
1 Like
Thank you for the bug report.
When I run your case I get the following error popup:
Do you see something similar?
The message says: “Cannot open file sq_input.dat”
I do not find any mention of such a file in the current MFiX documentation. This is a fairly new feature and it looks like the documentation has not caught up. Sorry for the difficulty.
Examining the code at des/red_particle_input.f at line 1557 (you can click the links in the error backtrace) I find this:
open(unit=1091,file='sq_input.dat', action='read')
do lcurpar = 1, particles
read(1091,*) dpar_super_r(lcurpar,1),dpar_super_r(lcurpar,2),&
dpar_super_r(lcurpar,3),dpar_super_mn(lcurpar,1),&
dpar_super_mn(lcurpar,2), dpar_super_q(lcurpar,1),&
dpar_super_q(lcurpar,2),dpar_super_q(lcurpar,3),&
dpar_super_q(lcurpar,4)
enddo
close(1091)
There is no check for presence of sq_input.dat or error handling if the file is not present, so the Fortran runtime error is raised.
Looking at the read statement, we see that 9 fields are read - the 3 semiaxes, followed by the m and n shape exponents, followed by 4 components of the quaternion describing the particle orientation.
This means that each particle could be a different size and shape. I am assuming you want them all to be the shape that you set up in the case file. (It would be nice to make the sq_input.dat optional. Clearly there is some room for improvement here.)
I tried the following:
for x in `seq 4500`; do echo 0.00635 0.00635 0.00175 2.0 8.0 1.0 0.0 0.0 0.0 ; done > sq_input.dat
which creates a sq_input.dat with 4500 identical lines - but running this case I get:
Error: Solver crash!
Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.
Backtrace for this error:
#0 sq_properties_mod_MOD_sq_volume
at des/sq_properties.f:675
and at line 675 I find:
670 SUBROUTINE SQ_VOLUME(axi,m,n,svolume)
671 IMPLICIT NONE
672
673 DOUBLE PRECISION :: axi(3),m,n,svolume,bt5, bt6
674 DOUBLE PRECISION :: e1, e2
675 e1 = 2.0D0/n
676 e2 = 2.0D0/m
somehow despite putting n=8 in the sq_input.dat file we are getting n=0 here resulting in a floating-point error. Strange.
Will follow up when I know what’s going on.
Thanks again for the report and for your patience. SQP is a new feature and it looks like there are still a few issues with it. Your feedback helps us improve MFiX.
– Charles
You are using a “Version 1” particle_input.dat which is just a list of numbers. The code in read_particle_input.f looks at the first line of that file to determine what format is used. A new format was defined in MFiX 20, which is indicated with Version 2.0 on the first line. This is described in the manual at:
The advantage of this format is that it allows to either read variables or set a constant value for all particles. That saves having to create files like the sq_input.dat file I constructed above.
read_particle_input.f actually defines 3 implementations, Version 1, Version 2, and Version 3, which doesn’t seem to be documented yet.
The Version 1 code always tries to read the superquadric parameters from sq_input.dat after loading particle_input.dat.
Version 2 does not care at all about superquadric parameters.
And Version 3 has support for SQP parameters being either constant or specified per-particle.
Let’s try upgrading your V1 file to V2:
Version 2 requires a long header, including the instructions, and expects data to start on line 35. So the easiest way is to start with an existing Version 2 file and modify it.
The DEM 3D hopper tutorial uses V2 format so I copied the header from there, and then modified it as follows:
Version 2.0
========================================================================
Instructions:
Dimension: Enter "2" for 2D, "3" for 3D (Integer)
Particles: Number of particles to read from this file (Integer)
For each variable, enter whether it will be read from the file
("T" to be read (True), "F" to not be read (False)). When not read, enter the
default values assign to all particles.
Coordinates are always read (X,Y in 2D, X,Y,Z in 3D)
Phase_ID, Diameter, Density, Temperature are scalars
Velocity requires U,V in 2D, U,V,W in 3D
Temperature is only read/set if the energy equation is solved
Species are only read/set if the species equations are solved, and requires
all species to be set (if phase 1 has 2 species and phase 2 has 4 species,
all particles must have 4 entries, even phase 1 particles (pad list with zeros)
User scalars need DES_USR_VAR_SIZE values
Data starts at line 35. Each column correspond to a variable.
========================================================================
Dimension: 3
Particles: 4500
========================================================================
Variable Read(T/F) Default value (when Read=F)
========================================================================
Coordinates: T Must always be T
Phase_ID: F 1
Diameter: T
Density: T
Velocity: T
Temperature: F (Ignored if energy eq. is not solved)
Species: F (Ignored if species eq. are not solved)
User_Scalar: F (Ignored if no user scalars are defined)
============================================================================================
X (m) Y (m) Z (m) Diameter (m) Density Velocity
============================================================================================
1.34940000E-02 1.10180000E-02 1.55820000E-02 6.42583409E-03 2.65000000E+03 0.00000000E+00 0.00000000E+00 0.00000000E+00
2.69890000E-02 1.10180000E-02 1.55820000E-02 6.42583409E-03 2.65000000E+03 0.00000000E+00 0.00000000E+00 0.00000000E+00
4.04830000E-02 1.10180000E-02 1.55820000E-02 6.42583409E-03 2.65000000E+03 0.00000000E+00 0.00000000E+00 0.00000000E+00
Note that I set Diameter, Density and Velocity enabled (T) because your file contains these. Alternately, you could remove all but the first 3 columns of data and specify the constant density and velocity values in the header.
After doing this I still get the same floating-point exception at
Floating-point exception - erroneous arithmetic operation.
Backtrace for this error:
#0 sq_properties_mod_MOD_sq_volume
at des/sq_properties.f:675
due to shape parameter N being 0.
It seems that when using particle_input.dat the values for particle parameters like diameter, density, etc which are set in the MFiX GUI are ignored completely. (Potentially confusing)
We need to use a Version 3 particle input file. Although the Version 3 format is not well documented yet, by looking at the functionWRITE_PART_OUT_V3P0 in read_particle_input.f you can see what the header and instructions look like. Now the data starts on line 39 instead of 35 and there are additional fields for SuperDEM parameters.
Here’s your file in version 3 format. I went ahead and removed the constant colums.
Version 3.0
========================================================================
Instructions:
Dimension: Enter "2" for 2D, "3" for 3D (Integer)
Particles: Number of particles to read from this file (Integer)
For each variable, enter whether it will be read from the file
("T" to be read (True), "F" to not be read (False)). When not read, enter the
default values assign to all particles.
Coordinates are always read (X,Y in 2D, X,Y,Z in 3D)
Phase_ID, Diameter, Density, Temperature are scalars
Velocity requires U,V in 2D, U,V,W in 3D
Temperature is only read/set if the energy equation is solved
Species are only read/set if the species equations are solved, and requires
all species to be set (if phase 1 has 2 species and phase 2 has 4 species,
all particles must have 4 entries, even phase 1 particles (pad list with zeros)
User scalars need DES_USR_VAR_SIZE values
Data starts at line 39. Each column correspond to a variable.
SQP_semi-axes, SQP_roundness, and SQP_quaternion are only needed when using SuperDEM.
========================================================================
Dimension: 3
Particles: 4500
========================================================================
Variable Read(T/F) Default value (when Read=F)
========================================================================
Coordinates: T Must always be T
Phase_ID: F 1
Diameter: F 6.42583409E-03
Density: F 2.65000000E+03
Velocity: F 0 0 0
Temperature: F (Ignored if energy eq. is not solved)
Species: F (Ignored if species eq. are not solved)
User_Scalar: F (Ignored if no user scalars are defined)
SuperDEM_semiaxis F 0.00635 0.00635 0.00175
SuperDEM_roundness F 2 8
SuperDEM_quaternion F 1 0 0 0
===========================================
X (m) Y (m) Z (m)
===========================================
1.34940000E-02 1.10180000E-02 1.55820000E-02
2.69890000E-02 1.10180000E-02 1.55820000E-02
4.04830000E-02 1.10180000E-02 1.55820000E-02
...
particle_input.dat (199.7 KB)
Using this input file, the case runs! And here’s a picture to prove it:
Look at those pretty non-spherical particles! 
But alas, after a while there is a new error:
Backtrace for this error:
#0 wrap_pow
at xpow.c:31
#1 sq_equivalent_radius_mod_MOD_sq_gradient
at des/sq_equivalent_radius.f:53
#2 sq_contact_newton_dpmethod_mod_MOD_func_dp_a
at des/sq_contact_detection_newton.f:949
#3 sq_contact_newton_dpmethod_mod_MOD_sq_contact_newton_dp_a
at des/sq_contact_detection_newton.f:267
#4 sq_calc_force_mod_MOD_calc_force_superdem
at des/sq_calc_force_superdem.f:341
#5 des_time_march_MOD_des_time_step
at des/des_time_march.f:192
#6 run_dem
at mfix.f:211
#7 run_mfix
at mfix.f:146
#8 main_MOD_run_mfix0
at main.f:79
we are trying to compute x*x where x=6.6205546496795669e+154 and it’s an overflow. It’s also Friday, I will return to this topic next week
– Charles
Thank you for your reply! @cgw
The fourth column of data in the particle_input.data file is radius rather than diameter. I’ll multiply the radius by two to see if it works.
For superquadric particles the radius/diameter is not used - the three semiaxes describe the particle shape.
Thank you for your reply! @cgw
According to your method, an error occurred when calculating to 0.34s on the Windows system
, and an error was reported at the beginning of the Linux system.
Yes, I see the same float overflow error.
If you try 23.1 you will get a better error popup on Windows (it will point to file and line, like it does on Linux).
Still not sure what the underlying cause of this overflow is.
Thank you for your reply! @cgw
I use automatic particle generation to calculate irregular particles, but after a period of time, there will be errors and the prompt “overlap between two super square particles is too large!” will continue throughout the calculation process. Can you help me see where the problem lies? Thank you very much！
n2d-lsfb_2023-04-24T085636.440765.zip (626.0 KB)
This might be caused by m=8; normally, we should set m<=6 for cylindrical particles or box particles.
If you use particle_input.dat for SQP, please use version 3. In above figure, it should be diameter (minimum bounding diameter), not radius. Note a, b, a are semi-axil of a superquadri particle.
Thanks for the comments @gaoxi.
We are using a Version3 file.
You are correct that the diameter needs to be doubled. But I think there is some room for improvement here - we have a bounding sphere calculation in both the GUI and the solver, why do we make the user specify it here?
And why do we make users specify constant values in particle_input.dat instead of deferring to the keyword values? That would be a much nicer fit with the GUI - user could set up particle shape, etc in the GUI and then that would be used. Having the particle_input.dat override the keywords seems error-prone and not very user-friendly.
Also, if m=8 is not stable, we should not allow users to specify it in the GUI. But rather than telling users not to use this shape (or change Young’s modulus, etc) I’d like to investigate why the solver is not stable with these values. I believe that, rather than adjusting the parameters to make the solver happy, we should try to get a better idea why the solver crashes and make the calculations more robust.
With the correct diameter value in the file (0.012787997374897409, as computed by the superquadric designer) the simulation crashes almost immediately with:
Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.
Backtrace for this error:
#0 sq_calc_force_mod_MOD_calc_force_superdem
at des/sq_calc_force_superdem.f:377
#1 des_time_march_MOD_des_time_step
at des/des_time_march.f:192
#2 run_dem
at mfix.f:211
#3 run_mfix
at mfix.f:146
#4 main_MOD_run_mfix0
at main.f:79
I’ve been experimenting with a version of xpow which temporarily disables overflow checking, so that instead of a value like 1e240 ^ 2 leading to an overflow, it returns inf.
With that modifcation, the simulation starts off by printing a huge number (191742) of overlap between two superquadric particles is too large! messages, and after that proceeds extremely slowly:
But it’s running. I’ll let it go overnight and see what happens.
I’ve made a little progress on this case and found some interesting things.
First, thanks again to @gaoxi for pointing out the radius/diameter mixup.
Comment: I think it would be best for this input to be optional. It can be computed from the semiaxes, and for non-particle_input.dat cases (automatic particle generation), if the diameter is not specified it is computed from the semiaxes. So why make users specify it here?
Once I corrected the value in the file, I noticed something interesting. In the VTK window:
when we look closely they do seem to be overlapping a bit.
This is quite strange because yesterday when I looked at the particles they were not touching at all:
I think the only thing that has changed is the diamater specified in particle_input.dat. But that should not affect how the particles are displayed in the VTK window, because, again, the superquadrics are completely described by the semiaxes and m,n exponents. Perhaps there is a bug in the visualization code, I will follow up on this.
So I decided to take seriously the thousands of
overlap between two superquadric particles is too large!
messages on the screen.
Looking at your data file, the minimum spacing between the particles in the y-direction is 0.0131989 which is very close to the bounding sphere diameter, it looks like it is close enough to cause trouble:
0.0131989 - 0.01278799737489740 = 0.000410
too close for comfort.
Here’s a little Python script that adds a little breathing room:
#!/usr/bin/env python
import numpy as np
input = open('particle_input.dat')
output = open('particle_input.dat.new', 'w')
# copy header
for lno in range(38):
output.write(input.readline())
# load numeric data
data = np.loadtxt(input)
# rescale
data *= [1, 1.1, 1]
# write output
for row in data:
row.tofile(output, sep=' ')
output.write('\n')
(Rename particle_input.dat.new to particle_input.dat after running this.)
After doing this the overlap is too large! messages went away, at least for a while. I am still getting a crash, but it runs a little further before crashing, and goes a lot faster (I am also running 8-thread SMP which helps a lot).
Here’s the particles with a little more space:
and here’s a little movie showing how far we’ve come:
More to come …
Disabling the FPE trap in xpow and allowing it to return inf:
Now we’ve got a new backtrace:
Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.
Backtrace for this error:
#0 sq_rotation_mod_MOD_qmatrixrotatevector
at des/sq_rotation.f:345
#1 sq_equivalent_radius_mod_MOD_sq_gradient
at des/sq_equivalent_radius.f:80
#2 sq_contact_newton_dpmethod_mod_MOD_func_dp_a
at des/sq_contact_detection_newton.f:949
#3 sq_contact_newton_dpmethod_mod_MOD_sq_contact_newton_dp_a
at des/sq_contact_detection_newton.f:267
#4 sq_calc_force_mod_MOD_calc_force_superdem._omp_fn.0
at des/sq_calc_force_superdem.f:341
Core was generated by `/home/cgw/mambaforge/envs/mfix-git/bin/python -m mfix_solver.pymfix -d /tmp/n2d'.
Program terminated with signal SIGFPE, Arithmetic exception.
#0 0x00007f4b2a8fcd3d in sq_rotation_mod::qmatrixrotatevector (qrot=..., x_frame_1=..., x_frame_2=..., ixg=2)
at /home/cgw/mambaforge/envs/mfix-git/share/mfix/src/model/des/sq_rotation.f:345
345 & *X_frame_2(2)+Qrot(3,1)*X_frame_2(3))
[Current thread is 1 (Thread 0x7f4b26b046c0 (LWP 17346))]
(gdb) list
340 ! its transpose
341 !
342 ! Rotate the vector "X_frame_2" from the local frame (2) to the
343 ! global frame (1)
344 X_frame_1(1) = 2.D0*(Qrot(1,1)*X_frame_2(1)+Qrot(2,1) &
345 & *X_frame_2(2)+Qrot(3,1)*X_frame_2(3))
346 X_frame_1(2) = 2.D0*(Qrot(1,2)*X_frame_2(1)+Qrot(2,2) &
347 & *X_frame_2(2)+Qrot(3,2)*X_frame_2(3))
348 X_frame_1(3) = 2.D0*(Qrot(1,3)*X_frame_2(1)+Qrot(2,3) &
349 & *X_frame_2(2)+Qrot(3,3)*X_frame_2(3))
(gdb) p qrot
$1 = ((0.49903154015346318, -0.0007951998096093725, 0.031094848276286433) (-0.028471787774220762, 0.18959288376215444, 0.46178338615382208) (-0.012525144031196185, -0.46265960066123457, 0.18918037605149429))
(gdb) p x_frame_2
$2 = (-inf, -inf, inf)
Reducing the number of particles from 4500 to 45 - this is probably a good idea for debugging because it speeds things up considerably. Here’s another movie, with the colors corrected (negative y-velocity, downward, is blue/purple, positive y-velocity is yellow). I used a shorter VTK write interval to capture more detail. It starts off slowly but gets pretty wild when the particles hit the bottom and start spinning around like crazy!
The error now is in a new place:
Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.
Backtrace for this error:
#0 0x7fc530f8f5cf in ???
#1 0x7fc4e4688e64 in __sq_calc_force_mod_MOD_calc_force_superdem._omp_fn.0
at model/des/sq_calc_force_superdem.f:377
Next things to try:
- disable ALL overflow checks (not just in
xpow.c) - try m=6 as Xi suggests.
With all FPE (overflow/zero-division) checks disabled (invalid operations will result in NaN or +/-Inf as appropriate, per IEEE754)
This is quite interesting - it took about 20 minutes to simulate 2 seconds, and there were no crashes, just some overlap is too big! messages.
Thes disc-shaped particles are spinning very fast. As entertaining as it is to watch, I’m not sure how physical it is. It might be correct. Particles tend to find their happy place and hover there, spinning madly, but eventually they drift off or are hit by another particle.
If you look closely, going frame-by-frame, you can find a few things like this:
two particles have completely intersected in a physically impossible manner. This is perhaps to be expected, given the high spin rate - we might have just missed the collision between them. dt is 0.01 sec, it’s possible that reducing max_dt would help with these overlaps (at the cost of slowing down the simulation). But since FPE checks are disabled, this was non-fatal - the simulation continued regardless.
This tells me that perhaps we need to reconsider the way we handle FPE’s are handled in MFiX. Disabling them altogether seems unwise - but for sections of code like the SQP code where we know large values are common (due to corners and other singularities) it might be easier to disable FPE checks than to rewrite the code to avoid doing things like (huge_number) ** 2 - checking all the inputs to all calculations would be a lot of work and could also slow down the code substantially.
Or perhaps we could add a control to the Numerics pane to enable/disable the various type of FPE check.
Next things to try:
- Change shape exponents
- Add more particles
Here’s a plot of simulation time vs. real time. Note that when particles are colliding, the slope is much lower. With fewer collisions to check, it speeds up.
Also I found that unfortunately there was a mistake in my previous runs - I had inadvertently put in the bounding diameter for m=2, n=6rather than m=2, n=8. Oops. (As I said, having to specify this is error-prone). The numbers are very close, 0.012851 vs. 0.012788 so perhaps this was not too big of an error. Running again with the correct value gives fairly similar results:
(note that the range of the colorbar for y-velocity has been expanded, so that small velocities are no longer over-emphasized)
and there are still a few overlaps:
Next: change exponent from 8 to 6