SQP performance

MFiX
,
5

Unfortunately the covergence criteria in the Numerics pane are for the fluid and DEM solvers. The convergence criteria for the SQP contact detection algorithm are not adjustible as of now. (I’ll see about adding this to a future release). The SQP code uses its own GMRES implementation in the file sq_gmres.f and the convergence criteria are all hard-coded there.

I did a little experiment and was able to get your case to run in half the time:

comparison
comparison1189×577 40.5 KB

There are still some flat spots in the “modified” run but they are shorter, as can be seen more clearly in this closeup of the plot:

comparison-detail
comparison-detail1189×577 30.8 KB

Here are the changes I made.

  1. Compile sq_gmres with higher level of compiler optimization. In order to do this you will have to edit the file CMakeLists.txt in the mfix source directory. This file cannot be copied to the work directory like Fortran files, you have to edit the master copy. It is in $CONDA_PREFIX/share/mfix/src
    Add these lines:
set_source_files_properties(${MFIX_SRC}/des/sq_gmres.f PROPERTIES
  COMPILE_OPTIONS "${CFLAGS};-O3;-ggdb;${MARCH}")

You can add these right below the similar setting for xpow.c (around line 295)

  1. Make the following changes in model/des/sq_contact_detection_newton.f . You can either edit the master copy or copy this file to your project directory. If you copy to project directory you will still have the original code for comparison.
diff --git a/model/des/sq_contact_detection_newton.f b/model/des/sq_contact_detection_newton.f
index 1f09d3f5d..1addf15d3 100644
--- a/model/des/sq_contact_detection_newton.f
+++ b/model/des/sq_contact_detection_newton.f
@@ -84,7 +84,7 @@ CONTAINS
 
           INTEGER,PARAMETER :: MAX_ITER = 100
           DOUBLE PRECISION,PARAMETER :: g = (sqrt(5.0) - 1.0) / 2.0 ! golden ratio
-          DOUBLE PRECISION,PARAMETER :: tol2 = 1e-7**2  !squared tolerance for golden mean
+          DOUBLE PRECISION,PARAMETER :: tol2 = 1e-5**2  !squared tolerance for golden mean
           DOUBLE PRECISION :: AAA(4), BBB(4), X1(4), X2(4), F1, F2, DELX(4)
           DOUBLE PRECISION :: funca1(4), funca2(4), funca01(4), funca02(4), FUNCA03(4)
           DOUBLE PRECISION :: X_A01(4), X_A02(4), X_A03(4)
@@ -249,7 +249,7 @@ CONTAINS
               else
                   ! solve linear equation using gmres_general
                   m = 1000
-                  threshold = 1E-10
+                  threshold = 1E-8
                   X_A_DELTA(:) = 0.0
                   call gmres(FJACOBIA, B, X_A_DELTA, 4,m, threshold, X_AA,c)
               endif
@@ -424,7 +424,7 @@ CONTAINS
 
           INTEGER,PARAMETER :: MAX_ITER=100
           DOUBLE PRECISION,PARAMETER :: G = (sqrt(5.0) - 1.0) / 2.0 ! golden ratio
-          DOUBLE PRECISION,PARAMETER :: TOL2 = 1E-7 ** 2 !squared tolerance for golden mean
+          DOUBLE PRECISION,PARAMETER :: TOL2 = 1E-5 ** 2 !squared tolerance for golden mean
           DOUBLE PRECISION :: AAA(4), BBB(4), X1(4), X2(4), F1, F2, DELX(4)
           DOUBLE PRECISION :: funcb1(4), funcb2(4), FUNCB01(4), FUNCB02(4), funcb03(4)
           DOUBLE PRECISION :: X_B01(4), X_B02(4), X_B03(4)

Then build a custom solver for your project.

We have an MFiX release coming out soon (25.3) but I do not feel comfortable adding these changes to the release without more testing. You can help us test this :slight_smile: I suspect that there are further optimizations possible, we may not need to compute the inverse matrix at all.

And this condition on the determinant > 1 has never been clear to me: (@gaoxi , any comment)?

              if (dabs(detA) > 1.0d0)  then
                  call matinv4(FJACOBIA, FJACOBIA_inv)
                  call matvec4(FJACOBIA_inv, B, X_AA)
              else
                  ! solve linear equation using gmres_general
                  m = 1000
                  threshold = 1E-8
                  X_A_DELTA(:) = 0.0
                  call gmres(FJACOBIA, B, X_A_DELTA, 4,m, threshold, X_AA,c)
              endif

why are we using the gmres method when |det(A)| <= 1 ? Why is the cutoff equal to 1?

It’s also possible there’s a better way to solve the linear equation (BiCGSTAB, etc). We could try allowing different linear solvers here like we do for the fluid solver, etc.

So there’s probably quite a bit we could to to speed this code up. But we have to be careful to not introduce excessive errors.

With the modified code, I notice a tendency for the particles to continue to move after they have settled - this seems a bit non-physical (kinetic energy seemingly coming from nowhere) - you can see this at around the 6s mark in this movie clip:

However I see very similar things happening with the unmodified code, so I’m not convinced that this is due to my changes:

Thanks for your interest in MFiX and for helping us test this code. Please let us know if you find anything interesting. We may be able to find some improvements for the next release.

– Charles