no progress of the solution 2

[ellipm]

Hi,
I'm working on a code for some months now and i have this problem: I can't forward the solution more than 2 steps than the initial steps. I think that the problem is that I set NYI+2 out of the boundaries,but I'm not sure. And as I increase the number of iterations for i, it seems that j reduces from nyi+2 to nyi-1.
This is a part of my code:

! initial conditions
X(1)=0.0
Z(1)=0.0
DX(1)=0.0
DZ(1)=0.0
THETA(1)=THETAO
DELTA(1)=RAD ! mixing layer thickness
DMIX(1)=0.0
YCORE(1)=RAD
!into the boundaries
DO J=1,NYI+1
U(1,J)=1.0
V(1,J)=0.0
DR(1,J)=20.0
EV(1,J)=1.0e-03
CV(1,J)=1.0e-03/0.7
AL1=2.0*EV(1,J)
AL2=2.0*EV(1,J)
ENDDO
!out of the boundaries
DO J=NYI+2,NJ
U(1,J)=UA
V(1,J)=0.0
DR(1,J)=0.0
EV(1,J)=1.0E-06
CV(1,J)=1.0E-06/0.7
AL1=2.0*KV
AL2=2.0*KV
ENDDO

!BOUNDARY CONDITION
DO I=1,NI
V(I,1)=0.0
ENDDO

DO 20 I=1,NI
DO 10 J=2,NJ
U(I+1,J)=U(I,J)-DS*V(I,J)*(U(I,J+1)-U(I,J))/(DY*U(I,J))+DS*((Y(J)+Y(J+1))/(4.0*U(I,J)*Y(J)*DY**2.0)*SMUP(I,J)*(U(I,J+1)-U(I,J))-(Y(J)+Y(J-1))/(4.0*U(I,J)*Y(J)*DY**2.0)*SMUM(I,J)*(U(I,J)-U(I,J-1)))+DS*G/U(I,J)*DR(I,J)/DENSO*SIN(THETA(I))

UMEAN(I+1)=(U(I+1,1)+UA)/2.0
IF(U(I+1,J).LE.UMEAN(I+1)) then
GOTO 32
else
goto 10
end if

32 YUHALF(I+1)=(J-2)*DY+((Y(J)-Y(J-1))/(U(I+1,J)-U(I+1,J-1)))*(UMEAN(I+1)-U(I+1,J-1))

!outer boundary condition
UTEST(I+1,J)=(U(I+1,1)-U(I+1,J))/(U(I+1,1)-UA)

IF(UTEST(I+1,J).GE.0.99) THEN
GOTO 27
ELSE
GOTO 10
END IF
10 CONTINUE

27 DELTA(I+1)=(j-1)*DY
DMIX(I+1)=DELTA(I+1)-YCORE(I+1)
write(*,*) 'delta(i+1)=',delta(i+1)
write(*,*) 'dmix(i+1)=',dmix(i+1)

NY=J
write(*,*) 'ny=',ny

DO J=NY+1,NJ
EV(I+1,J)=1.0e-06
U(I+1,J)=UA
ENDDO

If anyone can guess my mistake, because I have no more inspiration to find it, I'd really appreciate it!

 

[salgerman]

Quick feedback without looking beyond the first two loops.

What's the point of having
AL1=2.0*EV(I,J)
AL2=2.0*EV(I,J)
inside the first (into the boundaries) loop? It does not look like EV(I,J) is changing inside the loop at all...also, AL1 and AL2 are scalars, not arrays...what's the point of writing to them then? Just wait until you come out of the loop or make them scalars if that's what you need. How come they are two different variables, but get the exact same value?

What's the point of having
AL1=2.0*KV
AL2=2.0*KV
inside the second (out of the boundaries) loop? You just overwrote the values you set them to in the previous loops, etc, etc, same comments as before.

Also, please use the 'code' tags to surround code to ease reading it.

 

[ellipm]

I hope now it's better!

! initial conditions
X(1)=0.0           
Z(1)=0.0           
DX(1)=0.0
DZ(1)=0.0
THETA(1)=THETAO
DELTA(1)=RAD      ! mixing layer thickness
DMIX(1)=0.0   
YCORE(1)=RAD
!into the boundaries
  DO J=1,NYI+1            
     U(1,J)=1.0       
     V(1,J)=0.0       
     DR(1,J)=20.0       
     EV(1,J)=1.0e-03              
     CV(1,J)=1.0e-03/0.7
       AL1=2.0*EV(1,J)
     AL2=2.0*EV(1,J)
  ENDDO
!out of the boundaries
  DO J=NYI+2,NJ
     U(1,J)=UA         
     V(1,J)=0.0
     DR(1,J)=0.0                
     EV(1,J)=1.0E-06          
     CV(1,J)=1.0E-06/0.7
     AL1=2.0*KV
     AL2=2.0*KV
  ENDDO                  

!BOUNDARY CONDITION
  DO I=1,NI
       V(I,1)=0.0       
  ENDDO

DO 20 I=1,NI
DO 10 J=2,NJ
U(I+1,J)=U(I,J)-DS*V(I,J)*(U(I,J+1)-U(I,J))/(DY*U(I,J))+DS*((Y(J)+Y(J+1))/(4.0*U(I,J)*Y(J)*DY**2.0)*SMUP(I,J)*(U(I,J+1)-U(I,J))-(Y(J)+Y(J-1))/(4.0*U(I,J)*Y(J)*DY**2.0)*SMUM(I,J)*(U(I,J)-U(I,J-1)))+DS*G/U(I,J)*DR(I,J)/DENSO*SIN(THETA(I))
 
UMEAN(I+1)=(U(I+1,1)+UA)/2.0
IF(U(I+1,J).LE.UMEAN(I+1)) then
            GOTO 32
                else
               goto 10
               end if

32 YUHALF(I+1)=(J-2)*DY+((Y(J)-Y(J-1))/(U(I+1,J)-U(I+1,J-1)))*(UMEAN(I+1)-U(I+1,J-1))    
          
!outer boundary condition
       UTEST(I+1,J)=(U(I+1,1)-U(I+1,J))/(U(I+1,1)-UA)
       
IF(UTEST(I+1,J).GE.0.99) THEN
       GOTO 27
       ELSE
       GOTO 10
       END IF
    10 CONTINUE

 27  DELTA(I+1)=(j-1)*DY
       DMIX(I+1)=DELTA(I+1)-YCORE(I+1)
        write(*,*) 'delta(i+1)=',delta(i+1)
     write(*,*) 'dmix(i+1)=',dmix(i+1)
 
 NY=J
  write(*,*) 'ny=',ny
 
 DO J=NY+1,NJ
 EV(I+1,J)=1.0e-06
 U(I+1,J)=UA
ENDDO

Even without AL1,AL2 into the boundaries loop nothing changes!

 

[gummibaer]

because I have no more inspiration to find it

This would be the wrong approach. Chasing bugs is nothing of inspiration but of systematic work applying a certain method.

I must confess I am not very good in analysing foreign code when I do not know the theoretical background of it, especially when it is not so easy to be read when directly copied into the thread withoud using the code-tags.

But what would I do:

(1) get rid of the spaghetti code:

IF(UTEST(I+1,J).GE.0.99) THEN
       GOTO 27
       ELSE 
       GOTO 10
       END IF
    10 CONTINUE   ! this is the end of a do loop

 27  DELTA(I+1)=(j-1)*DY
       DMIX(I+1)=DELTA(I+1)-YCORE(I+1)
        write(*,*) 'delta(i+1)=',delta(i+1)
     write(*,*) 'dmix(i+1)=',dmix(i+1)

Structures like this are very hard to follow in debugging. It might be hard at first but you can code without ever using 'goto'.
And once you achieved this, you will find that you can follow execution your code much easier as before.

(2) use 'implicit none' in your code. Should be before the type declaration section. This forces you to declare the types of all the variables you use - but there is no better way to elimit typos from your code.

Then you should find answers to
-> Which is the loop that does the iteration (the do 20 or the do 10 ?)
-> which variable sets how often this loop is executed ? Or is execution terminated because some limit value is reached ?
-> where in your code is this variable / value determined ?
-> have this variable printed before the loop is started or before it is checked
-> If the value is not what you expect it to be, find the piece of code where it is determined and printout all the variables that contribute to it.
-> see which of those is outside of what to expect
-> find out why this is so ....

and so on and so on.

It might take some time - but you will sure find the reason why your code does not perform as it should

The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

 

[gummibaer]

Today I have had some more time to look at your code (it was late at night / early in the morning in my location when I entered my response).

I guess the loop that should do the iteration is the loop

do 10 j = 2, NJ
.
.
.

10 continue

If so, then you have two conditions that stop execution of the loop, (1) the loop counter J has counted up to NJ, (2) UTEST (I+1,J) exceeds 0.99. Can you make sure which one terminates execution ?

To get a better grip on the code's flow, I would add a third recommendation to the two ones of yesterday:

(3) break down the long formula for U(I+1,J) into a set of smaller formulas, that would bring some meaningful intermediate results that you can check on.

To rearrange your code I would propose

do 20 I = 1, NI                 ! I do not see the end of the loop, so ileave it as it is
j = 0                           ! I would use J as a counter, not as a loop control
do while (UTEST(I+1,J).LT.0.99) ! this terminates your loop once the criterion is reached
    j = j + 1                   ! count upwards
    if (j .gt. 1000000) then    ! to avoid endless execution set a reasonable number here
         write (*,*) 'loop terminated after 10000000 cycles'    
         exit
    endif
!
!   compute U(I+1,J) in reasonable steps with results you can check

    UMEAN = .....
!
!    I understand the first of your ifs within the loop, that you want to start a new cycle if U (I+1,J) is less than UMEAN (I+1)
!    otherwise to continue with the loop's execution
!
     if (U(I+1,J).GT.UMEAN(I+1)) cycle

     UHALF = ....
     UTEST = ....

     write (*,*) 'J, UTEST = ', J, UTEST    ! for debugging
!
!    I understand the second if, that you want to terminate the loop's execution when UTEST(I+1,J).GE.0.99. 
!    This is the effect of the do while (...) at the beginning of the loop
!
enddo               ! terminate the loop then

DELTA = .....
DMIX = ....

etc.
etc.

By this modification you have only one termination criterion for your loop and you can track how it does develop as a basis for hunting down this bug.

Norbert


The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

 

[ellipm]

thank you so much gummibaer!I think that your first post was more helpful than the 2nd one. It seems that there's some progress now.But I'll try the solution that you propose on your 2nd post,too.

 

[salgerman]

ellipm:

This is the part where you go click on the link below the post the helped you the most as a way of thanking and acknowledging the author.