I think I have found a bug in the implementation of the cyclic boundary condition using GSP particles. I have attached a case where the bug if obvious if you run the case for an iteration. Although there are .f files in the project folder, you can just use the default (native) solver to see the bug. I am running the new mfix 26.1.2.
The problem is when the cyclic boundary condition is used then the estimated drag coefficient is 0 for particles next to the cyclic boundary condition. I have illustrated in the figure below. The BC is cyclic in the x-direction. As can be seen the outer most particles have 0 C_D. The result is this if I run with the project solver, or the default (native) solver using serial or SMP.
The problem is the same for cyclic in the z-direction. Either I am forgetting something or maybe it can be fixed in a future version. Would be nice for cases like this where SMP is preferred since it allows me to use more cores since DMP is limited to 3 computational cells in each direction.
Further experimentation showed it is only a problem when:
The problem is with the evaluation of the EP_G field. If I print various values from a usr_drag.f function I get the following where it is clear that the 0 drag coef lDgA is caused by EPg being evaluated as ~10^31 which I assume is since its undefined or something.
2. Note this problem only occurs if field-to-particle interpolation is used. i.e. des_interp_on = True. It happens for both options “Field-to-particle” only and “Field-to-particle and particle-to-field”
I assume it means it is a problem in the comp_mean_fields.f and maybe sendrecvnode_mod.f. I am a bit in over my head here but I will try to look in to it. I noticed this code in comp_mean_fields1.f
! Sum data interpolated into ghost cells into physical cells
!---------------------------------------------------------------------//
IF(DES_INTERP_ON) THEN
CALL DES_COLLECT_gDATA(SOLVOLINC(:,mLB:mUB))
ENDIF
And when I find the DES_COLLECT_gDATA there seems to be a lot of stuff about node count i.e. the dmp solver where again the EP_G problem does NOT happen. Maybe that is a clue
I think I got somewhere by using claude for a while. It suggests that the bug is in line 115 of calc_interp_weights.f. See image below. If I use the original line I get the bug and if I use the new code with MAP_C it seems to work. I need to do some further validation. From reading the code I gather that the EP_G evaluation goes wrong sinec the wrapping for cyclic boundaries is not present in FUNIJK. And that when using DMP then the problem is fixed for nodesi>1 since then the interpolation is handled elsewhere in the sendrecv file.
I included a .zip with a minimal example. You can compile with the old line and the new line to see the difference. It is also important to add FUNIJK_MAP_C to the openmp parallel clause which is also done in this file.
