photoreactor_1 - Copia (3).zip (2.3 MB)
Hello everyone, I am currently working on photocatalysis in fluidized bed reactors and I was wondering if anyone ever tried to implement its modeling with MFiX. I was looking into how it could be implemented, but I have no idea on how to access user subroutines or anything like that. Can somebody help me?
photoreactor_kinetics user.zip (32.5 MB)
Ok, I have started writing some code with an equation for light penetration in the bed. I was planning on using the Scalar to introduce initial condition for illumination, by setting the intensity of radiation equal to the Scalar value at the start of simulation with usr0.f, then calculating the actual intial condition also in usr0.f, in usr1.f I wanted to write the calculation of attenuation factor and irradiation and output them using the dummy reaction rate matrix. The lamp is simulated at the moment by the region filled with gas and wall on the side of z=0 (I can’t manage to simulate it simply with wall conditions). I doesn’t seem to calculate any value for the attenutaion factor and the light intensity in the vtk, they are always zero. What am I doing wrong?
I can’t comment on the attenuation factor calculation, but there are a few issues I can see:
ReactionRates array. Before using it, you need to allocate by specifying its size (bottom of the Model setup pane, set ReactionRates array size to 2 if you plan to store 2 quantities per cell). To output these arrays in the vtk files, go to Output>VTK>Reactions tab for a cell vtk file, and check the corresponding check boxes. This is why you don’t see anything in the vtk files.IJK loop in usr0.f, local_in loop in usr1.f)local_in_prev is not assigned any value so it could be zero or any garbage value. This should trigger an out of bound array error.Thank you for your answer. I did some changes to the code and I made some progress, but now I get this overflow error that I don’t know how to fix.
In usr_1, it should calculate the attenuation factor, and then calculate the light in the cell based on the light on the preceding cell in the z direction (the lamp is the space before the impermeable wall I set as internal surface (wall boundary conditions did not work, so I included it as an internal surface). The problem here is that it doesn’t seem to like the operation I_temp=Light_in*att_ijk in usr1, where Light_in should be a temporary variable that I use to store the value of light in the previous cell in the k direction and att_ijk is the value that it should calculate for the attnuation factor in the cell at the current time.
photoreactor_kinetics user (2).zip (33.0 MB)
If you are on Linux and you get a core file you can try to use GDB to see what values are overflowing. But from your bug report I see you are on Windows where unfortunately the debugging support is lagging a bit. Your best bet is to add some “write” statements before the line that fails. The error message is telling you the overflow is at line 112 in usr1.f:
107 IF (Ep_g(local_in)==1) THEN
108 I_temp=Light_in*1
109 Inten(local_in,1)=I_temp
110 Light_in=I_temp
111 ELSE
112 I_temp=Light_in*att_ijk
113 Inten(local_in,1)=I_temp
114 Light_in=I_temp
So before line 112, add
write(*,*) "Light_in=", light_in, "att_ijk=", att_ijk
I also ran this by Claude (Opus 5) which found other problems. I cannot guarantee this is all correct but it looks right to me. (You don’t have to worry about item 3 if you’re not planning on running in DMP mode).
IF (K_OF(IJK_2)<=nint(lamp_z)) THEN
att(local_in,1) = 1 ! <-- local_in, not IJK_2
local_in has no value at this point in the routine (it’s first assigned in the triple loop below, and it’s not SAVEd, so it’s indeterminate on every call). You’re writing 1 into att at a garbage index — out-of-bounds write, or silent clobbering of an unrelated cell.
Consequence twofold: att(IJK_2,1) is never set for the lamp region, so the following ReactionRates(IJK_2,1) = att(IJK_2,1) reads stale/uninitialized data.
Light_in is a running product carried across the entire triple loop, but the nesting is K outer, I, then J inner. For K > lamp_z the value arriving at cell (i,j,k) has been attenuated by every cell visited before it in scan order, not by the column above it. Beyond a few levels it underflows to zero.
The march has to be innermost in K, with the reset per column:
DO I_usr2 = Istart3, Iend3
DO J_usr2 = Jstart3, Jend3
Light_in = I_lamp
DO K_usr2 = Kstart3, Kend3
IF (DEAD_CELL_AT(I_usr2,J_usr2,K_usr2)) CYCLE
local_in = FUNIJK(I_usr2,J_usr2,K_usr2)
IF (K_usr2 > nint(lamp_z)) Light_in = Light_in * att(local_in,1)
Inten(local_in,1) = Light_in
ReactionRates(local_in,2) = Light_in
END DO
END DO
END DO
(The Ep_g==1 branch collapses out — it’s Light_in*1, and att is already 1.0 there.)
Kstart3..Kend3 are rank-local. With any decomposition in z, each rank restarts its column at I_lamp and the attenuation accumulated upstream is lost — results depend on NODESK. There’s no communication here to fix that; the march either needs a serialized sweep along the z-pencil (send the exit intensity plane to the next rank up) or the decomposition has to be constrained to NODESK=1. Also, the loop includes ghost layers, so boundary cells get counted twice into the product.
att is a per-cell transmittance with no length in itatt = -6.0D-3 + Ep_g*ep_new2 is applied as a bare multiplicative factor per cell. That makes the total attenuation att**(number of cells) — refine DZ and the light profile changes completely. Beer–Lambert wants exp(-mu*DZ(K)) with mu an absorption coefficient in 1/m. As written, it’s not grid-convergent.
Related: nothing bounds att. Ep_g*ep_new2 - 0.006 can be negative (sign-flipping intensity) or, if ep_new2 > 1, greater than 1 (light amplifying downstream). Clamp to [0,1] at minimum.
IF (J_OF(IJK_1)<=0.02)
J_OF returns the integer J index. Compared against 0.02, this is true only for J ≤ 0 — i.e. the ghost layer, never a real cell. Whoever wrote it meant a y-coordinate of 2 cm. You want the cell-center y (from YN/DY in geometry), not the index.
nint(lamp_z) is the same ambiguity in the other direction: if lamp_z is a height in meters, rounding it and comparing to K_OF is meaningless; if it’s meant to be a cell index it should be an integer in usr_mod. Worth resolving which it is before anything else here matters.
IF (Ep_g(IJK_1) /= Ep_g(IJK_1)) is a NaN test. Substituting a plausible-looking value for a NaN void fraction hides a solver failure rather than reporting it — the run continues producing garbage. Also, under -ffast-math / -ffinite-math-only gfortran is entitled to fold this to .FALSE., so the guard may not even exist in the optimized build. If you want the check, use IEEE_IS_NAN from IEEE_ARITHMETIC and call MFIX_EXIT on it.
IF (Ep_g(IJK_2)==1) — float equality. Use Ep_g >= ONE - SMALL_NUMBER.DOUBLE PRECISION, DIMENSION(DIMENSION_3,1) :: my_epg — automatic array on the stack, sized by the full local grid, reallocated every time step. Make it rank-1, and either SAVEd/allocatable or module-level. The trailing ,1 extent serves no purpose.ReactionRates(:,2) requires NRR >= 2 in the project file; silently out of bounds otherwise.local_in_prev, bjbjbj, I_temp (once the branch collapses), and the USE of ic, rxns, param.IJKSTART3..IJKEND3 loops can be fused; my_epg doesn’t need to exist as an array at all if att is computed in one pass.The two that will actually bite you today are #1 and #2. #3 is what will make the answer change when you switch node counts.
Final comments (from me, not Claude)
Error from check_data/check_bc_geometry.f:350
Error 1100: Invalid location specified for BC 14.
X: 0.35000E-02, 0.66500E-01 I: 6, 68
Y: -0.50000E-02, 0.75000E-01 J: 2, -1
Z: 0.0000 , 0.0000 K: 1, 1
This usually occurs when the BC region is outside the fluid region,
or when the BC region is smaller than the grid spacing.
Fatal error reported on one or more processes. The .LOG file
may contain additional information about the failure.
Hope this helps,
– Charles
Thank you for your answer, I am very much not used to writing in Fortran. For the I, J, K loop I indeed suspected that I was doing it in the wrong order, thanks for confirming it. And thank you for signaling that I forgot to bind the value of the attenuation between 0 and 1, I indeed already corrected it on my own in a later version, I am testing the other modifications. I don’t know why it doesn’t let you even start the run, I did not have that kind of problem.
I made the changes you suggested and now it seems to be working correctly. Now I’ll try connecting the light intensity to the kinetic calculations. Thank you again for your help, and sorry for the terrible indentation in the code, I fixed that too now.