[Alf90]

Hi everyone,

I have some problems using a self written Fortran user function for inlet total pressure in CFX.

I have written an iterative function, that is changing my inlet total pressure depending on downstream conditions. The script itself seems to work fine. I monitored the output (the value I set in CFX as inlet total pressure) of the script in CFX and it is exactly what I wanted it to be.

The problem is, that the total pressure at the inlet is not even close to this value. CFX is placing walls at 99%-100% of the inlet. I tried some different things until I'm now at a point, where I have no idea how to go on.

The last try was having a CEL function that has a constant total pressure for the first 50 iterations and after that switching to my function. In my function the total pressure was also hold constant for another 50 iterations (up until time step 100) with the same value as before. I have a monitor point for that CEL function that switches between the constant value and my function, that gives back the same constant value. The output of this monitor point is perfectly constant from timestep 50 to 51.
Up to timestep 50 the massFlowAve total pressure at the inlet is close to the constant value I set up (without any walls put at the inlet). At the point I switch to my script, the massFlowAve total pressure jumps within 2-3 timesteps to -10^20 Pa, which leads to a crash... I'm not wondering why it crashes, that is clear with that jump. But I have no idea why the input value remains constant but for the simulation this makes such a huge difference...

It is all done in a steady state simulation (Iteration is only done to match given conditions that are measured downstream) with Version 2024R1. This is my first time programming and using a user routine in CFX, so maybe there is something I completely misunderstood...
Unfortunately I cannot share a picture of the model or a full CCL. If it is not clear what I mean, let me know and I will build a simplified version/picture.

Thanks!

[Opaque]

You may need to post your "User Fortran". It would be impossible to know if there is an illegal memory access or an incorrect setting in your code w/o looking at it.

[Alf90]

The code is shown below.
Im running this in "Mode==1". I deleted the part for other modes to make it a bit shorter

Interval is 100

ValIni is the same value, I use before in the CEL equation.

FileNo is 1. This option is just there so I could use this script for different iterations within one run (if needed). In this simulation I'm using it only for my total inlet pressure.

When it crashes at timestep 53, I think all other input values should not be relevant.

#include "cfx5ext.h"
dllexport(usr_iterative_bc)
      SUBROUTINE USR_ITERATIVE_BC (
     & NLOC,NRET,NARG,RET,ARGS,CRESLT,CZ,DZ,IZ,LZ,RZ)
      IMPLICIT NONE
C
C ------------------------------
C        Argument list
C ------------------------------
C
      INTEGER NLOC, NRET, NARG
      CHARACTER*(4) CRESLT
      REAL    RET(NLOC,NRET), ARGS(NLOC,NARG)
C
      INTEGER IZ(*)
      CHARACTER CZ(*)*(1)
      DOUBLE PRECISION DZ(*)
      LOGICAL LZ(*)
      REAL RZ(*)
C
C ------------------------------
C        Local Variables
C ------------------------------
      CHARACTER (len=16) FNI
      INTEGER Iteration, Interval, Mode, Direc, FileNo
      DOUBLE PRECISION ValIni, ValDif, MaxStepAbs, MaxStepRel, NewVal
      REAL InData, InIter
      LOGICAL exists
C ---------------------------
C    Executable Statements
C ---------------------------
      Iteration=ARGS(1,1)
      Interval=ARGS(1,2)
      ValIni=ARGS(1,3)
      ValDif=ARGS(1,4)
      Mode=ARGS(1,5)
      MaxStepAbs=ARGS(1,6)
      MaxStepRel=ARGS(1,7)
      Direc=ARGS(1,8)
      FileNo=ARGS(1,9)
            
      IF (FileNo==1) THEN
         FNI='IterValLastA.txt'
      ELSEIF (FileNo==2) THEN
         FNI='IterValLastB.txt'
      ELSEIF (FileNo==3) THEN
         FNI='IterValLastC.txt'
      ELSE
         FNI='IterValLast.txt'
      END IF
          
      
      IF (Mode==1) THEN
         INQUIRE (FILE = FNI, EXIST = exists)
	 
	 IF (.NOT. exists) THEN
	    NewVal=ValIni
	    
            OPEN (UNIT=7, FILE = FNI, ACTION='WRITE')
            WRITE(7,*) (Iteration, NewVal)
            CLOSE(7)
         ELSEIF (Iteration<2) THEN
            NewVal=ValIni
	    
            OPEN (UNIT=7, FILE = FNI, ACTION='WRITE')
            WRITE(7,*) (Iteration, NewVal)
            CLOSE(7)
         ELSE
            OPEN (8, FILE = FNI, STATUS='OLD', ACTION='READ')
            READ(8,*) InIter, InData
            CLOSE(8)
            IF (InIter<Iteration) THEN
               IF (MODULO(Iteration,Interval)==0) THEN
	       
                  NewVal=InData+(ValDif/abs(ValDif))*
     & min(abs(ValDif)*MaxStepRel,MaxStepAbs)*Direc
	          
                  OPEN (UNIT=7, FILE = FNI, ACTION='WRITE')
                  WRITE(7,*) (Iteration, NewVal)
                  CLOSE(7)
               ELSE
                  NewVal=InData
               END IF
	    ELSE
	       NewVal=InData
            END IF
         END IF
      ELSE
         NewVal=ValIni
      END IF
      
      RET(1,1) = NewVal
      
      END

[ghorrocks]

I would do some basic debugging: Put lines in your code which output the values of key variables (either to the CFX output file or a separate output text file) and look at what all your variables are doing each iteration to check it is doing what you expect.

Your line: NewVal=InData+(ValDif/abs(ValDif))*
& min(abs(ValDif)*MaxStepRel,MaxStepAbs)*Direc

Looks very strange. What is this line trying to achieve? Is this going to make the value of NewVal jump around? If so this is always going to be hard to numerically converge. You might need some under-relaxation or other numerical techniques to get it to converge.

[Opaque]

Your code seems to assume to be called with an NLOC = 1 at all times.

Are you certain NLOC is 1 for all calls?

You need IF (NLOC .NE. 1) STOP, and see if that code works.

[Alf90]

Thank you for your responds. I was a week out of office and therefore had no chance to answer.

Originally Posted by ghorrocks

I would do some basic debugging: Put lines in your code which output the values of key variables (either to the CFX output file or a separate output text file) and look at what all your variables are doing each iteration to check it is doing what you expect.

Your line: NewVal=InData+(ValDif/abs(ValDif))*
& min(abs(ValDif)*MaxStepRel,MaxStepAbs)*Direc

Looks very strange. What is this line trying to achieve? Is this going to make the value of NewVal jump around? If so this is always going to be hard to numerically converge. You might need some under-relaxation or other numerical techniques to get it to converge.

I already tried that and wrote every value into external files. Every value was as I expected it to be. As shown in the code, "NewVal" is saved in an external file after every change and the value in this file is identical to my monitorpoint I set up in CFX, which also behaves how I want it to do.

The idea behind this line is to change my inlet condition regularly but not with every iteration. I want to set a constant total pressure at an interface downstream of the inlet. Therefore I wait several iterations so that my changes at the inlet can propagate to this interface and become (more or less) stable at this position. To avoid large jumps at the inlet, the change is limited by the given value "MaxStepAbs".
I had a similar code as a CEL function, which worked fine. There was just the problem, that I coudn't hold the value stable since I had no informations of prior timesteps. With the external file I tried to eliminate this problem.

Originally Posted by Opaque

Your code seems to assume to be called with an NLOC = 1 at all times.

Are you certain NLOC is 1 for all calls?

You need IF (NLOC .NE. 1) STOP, and see if that code works.

I thought that this would be called with NLOC=1 but it seems that this is not the case. NLOC is in the order of magnitude of the number of nodes at the inlet (where this function is used for inlet pressure) but it is also not exactly that number.
At this point I don't know how the solver interacts with the script and how it differs from just using a CEL expression.

What I want to get, is a single pressure value to use as inlet condition.

[Opaque]

Originally Posted by Alf90

I thought that this would be called with NLOC=1 but it seems that this is not the case. NLOC is in the order of magnitude of the number of nodes at the inlet (where this function is used for inlet pressure) but it is also not exactly that number.
At this point I don't know how the solver interacts with the script and how it differs from just using a CEL expression.

What I want to get, is a single pressure value to use as inlet condition.

The UserFortran function will be called as many times as the solver requires it. The inlet mesh is divided into groups of faces, and it will be called for each group in turn as needed.

As a debugging step, you should start by running the NULL case, i.e. send back EXACTLY the same value for all calls. Separately, run the case using the same value w/o UserFortran. Both simulations should be identical at every iteration; otherwise, the UserFortran code is introducing some error.

[Alf90]

I already read that the solver is calling the script as often as it needs it but I expected that this would be with NLOC=1, like "asking" for a single value as it would be when I just directly type in a single value for inlet pressure.


I just tried this by setting the value of "Mode" to 3. Therefore, it should give "ValIni" back. "ValIni" is a constant value without any change during the process. The behavior in the simulation is similar as before.
In parallel I also used the same value directly as input (without any CEL or FORTRAN code involved). It behaves different and in a way I expected it to do.

So there seems to be a problem with my script. For me, the problem seems to be in the way how my script is communicating with the solver, since the values I write to the external file are exactly as I want them to be...
I'll have a look into this when I have some time. If anyone has an idea what to change it would be great. Tanks a lot

[Opaque]

Originally Posted by Alf90

I already read that the solver is calling the script as often as it needs it but I expected that this would be with NLOC=1, like "asking" for a single value as it would be when I just directly type in a single value for inlet pressure.


I just tried this by setting the value of "Mode" to 3. Therefore, it should give "ValIni" back. "ValIni" is a constant value without any change during the process. The behavior in the simulation is similar as before.
In parallel I also used the same value directly as input (without any CEL or FORTRAN code involved). It behaves different and in a way I expected it to do.

So there seems to be a problem with my script. For me, the problem seems to be in the way how my script is communicating with the solver, since the values I write to the external file are exactly as I want them to be...
I'll have a look into this when I have some time. If anyone has an idea what to change it would be great. Tanks a lot

Perhaps there is a misunderstanding on your part on how the software work. When you specify a value in a boundary condition widget, you are not specifying the value on a single face, but several faces at once. That is why the "value" can an expression with spatial variation; otherwise, you would not be able to specify a distribution in space.

In your code, check how many locations are being provided -> NLOC. Your responsibility is to fill in ALL the entries; otherwise, the returned array is uninitialized.

[Alf90]

In general, I know that, but to be honest, in this case I perhaps didn't think about it properly... And that was the problem here. If I return not just a single value but a vector, everything now works the way I want it to.
I don't know what cfx does when the wrong number of values are returned, but if the missing entries are treated as zero, it would make perfect sense that walls were built on my input.

I feel a bit dumb now after seeing what went wrong, but many thanks for the help anyway.