Hello,everyone!
Add code to the vtk_out.f file to output the granular temperature in the grid, but the operation failed. Strangely, the same code can be successfully output in version 21.3.2. Do you know why?
fb_sqp_2023-04-27T195758.338235.zip (201.2 KB)
I’m not sure I understand your intention here.
There is an MFIX keyword VTK_THETA_M which will turn on saving granular temp in a region. But per the documentation:
- Enable writing solids phase granular temperature
- Requires TFM solids and KT_TYPE /= “ALGEBRAIC”
- Sets keyword VTK_THETA_M(#,#)
- DEFAULT value .FALSE.
Your case is a DEM case so VTK_THETA_M is unavailable.
In vtk_out.f you have:
DO M = MMAX+1,DES_MMAX+MMAX
IF(VTK_Theta_m(VTK_REGION,M).AND.ALLOCATED(Theta_m)) THEN
!#############################################################################################
CALL DES_GRANULAR_TEMPERATURE
!#############################################################################################
WRITE(SUBM,*)M
CALL WRITE_SCALAR_IN_VTU_BIN('Granular_temperature_'//ADJUSTL(SUBM),Theta_m(:,M),PASS)
ENDIF
END DO
That is inside an IF(VTK_THETA_M so unless that keyword is set you will not trigger that code.
But in your .mfx file you somehow have:
vtk_theta_m(1,1) = .True.
vtk_theta_m(1,2) = .True.
how did you set those for a DEM case? Did you edit the file outside of MFiX?
Here’s a summary of your changes to vtk_out.f
diff -u ~/Work/NETL/mfix/model/cartesian_grid/vtk_out.f /tmp/fb_sqp_2023-04-27T195758.338235
--- /home/cgw/Work/NETL/mfix/model/cartesian_grid/vtk_out.f 2023-03-19 08:40:05.706646717 -0500
+++ /tmp/fb_sqp_2023-04-27T195758.338235/vtk_out.f 2023-04-27 11:37:57.718384657 -0500
@@ -52,6 +52,9 @@
USE vtp, only: write_vtp_file, update_frames
USE, intrinsic :: iso_c_binding
USE WRITE_OUT0_MOD, only:location
+!###########################YIN################################
+ USE DES_GRANULAR_TEMPERATURE_MOD, only: DES_GRANULAR_TEMPERATURE
+!##############################################################
CONTAINS
@@ -190,7 +193,7 @@
IF(VTK_P_star(VTK_REGION).AND.ALLOCATED(P_STAR)) &
CALL WRITE_SCALAR_IN_VTU_BIN('P_STAR',P_STAR,PASS)
- DO M = 1,MMAX
+ DO M = MMAX+1,DES_MMAX+MMAX
IF(VTK_P_s(VTK_REGION,M).AND.ALLOCATED(P_S)) THEN
WRITE(SUBM,*)M
CALL WRITE_SCALAR_IN_VTU_BIN('P_S_'//ADJUSTL(SUBM),P_S(:,M),PASS)
@@ -213,14 +216,14 @@
ENDIF
END DO
- DO M = 1,MMAX
+ DO M = MMAX+1,DES_MMAX+MMAX
IF(VTK_ROP_s(VTK_REGION,M).AND.ALLOCATED(ROP_S)) THEN
WRITE(SUBM,*)M
CALL WRITE_SCALAR_IN_VTU_BIN('Solids_bulk_density_'//ADJUSTL(SUBM),ROP_S(:,M),PASS)
ENDIF
END DO
- DO M = 1,MMAX
+ DO M = MMAX+1,DES_MMAX+MMAX
IF(VTK_RO_s(VTK_REGION,M).AND.ALLOCATED(RO_S)) THEN
WRITE(SUBM,*)M
CALL WRITE_SCALAR_IN_VTU_BIN('Solids_density_'//ADJUSTL(SUBM),RO_S(:,M),PASS)
@@ -308,8 +311,11 @@
END DO
END DO
- DO M = 1,MMAX
+ DO M = MMAX+1,DES_MMAX+MMAX
IF(VTK_Theta_m(VTK_REGION,M).AND.ALLOCATED(Theta_m)) THEN
+!#############################################################################################
+ CALL DES_GRANULAR_TEMPERATURE
+!#############################################################################################
WRITE(SUBM,*)M
CALL WRITE_SCALAR_IN_VTU_BIN('Granular_temperature_'//ADJUSTL(SUBM),Theta_m(:,M),PASS)
ENDIF
I don’t understand your changes to the loop indices.
When I run this I get a floating-point exception … I assume this is what you meant by “The operation fails” (please try to be a more explicit about errors and include any stack traces you get)
Here’s the backtrace:
Error: Solver crash!
Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.
Backtrace for this error:
#0 vtk_out_mod_MOD_write_scalar_in_vtu_bin
at fb_sqp_2023-04-27T195758.338235/vtk_out.f:1289
#1 vtk_out_mod_MOD_write_vtu_file
at fb_sqp_2023-04-27T195758.338235/vtk_out.f:324
As always, write statements are invaluable when debugging:
WRITE(SUBM,*)M
write(*,*) "HERE", theta_m(:,M)
CALL WRITE_SCALAR_IN_VTU_BIN('Granular_temperature_'//ADJUSTL(SUBM),Theta_m(:,M),PASS)
results in
0.0000000000000000 0.0000000000000000 0.0000000000000000 9.7546105778997083E+063 9.7546105778997112E+063 9.7546105778997083E+063 9.7546105778997083E+063 9.7546105778997083E+063 9.7546105778997112E+063
9.7546105778997083E+063 9.7546105778997083E+063 9.7546105778997083E+063 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 9.7546105778997083E+063 9.7546105778997112E+063
9.7546105778997083E+063 9.7546105778997083E+063 9.7546105778997083E+063 9.7546105778997112E+063 9.7546105778997083E+063 9.7546105778997083E+063 9.7546105778997083E+063 0.0000000000000000
0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.00000
This looks like garbage data, note all the values are either 0 or 9.7546105778997112E+063
And I suspect VTK is trying to stick this into a 32-bit variable:
>>> import numpy as np
>>> x=9.7546105778997112E+063
>>> np.float64(x)
9.754610577899711e+63
>>> np.float32(x)
inf
the value is too big for 32 bits, triggering the exception you are seeing.
The subroutine DES_GRANULAR_TEMPERATURE is not supported. It is turned off in the code with a ! Invoke at own risk comment so my guess is it is not usable as is (it was originally written in 2004!). You may be better off writing you own routine.
Thank you for your help. @cgw
Let me take a closer look at the code. The same code can successfully output particle temperature in MFIX-21.3.2, which may be due to incompatible code versions.
Hello, @cgw
I found that it is related to the interpolation algorithm. Using the Garg2012 algorithm can output the granular temperature of particles, while using the Satellite DPVM algorithm will generate an error when building. Do you know which modification can be applied to the Satellite DPVM algorithm? Thank you！