Handling underflow errors in fortran 90

[meteo student]

I am using fortran 90 to solve a elliptic partial differential equation numerically(Poisson's equation) and I have encountered a underflow error. All the arrays involved are stored as real data types. The original expression is highlighted in the comments and the my attempt at refactoring into simpler expressions is shown as expr0, expr1, expr2,expr3. I am not sure about my attempt to handle the underflow error because while the FPE problem goes away the solution does not converge and I get all zero values.

Any suggestions to improve my coding attempts ?

    use,intrinsic :: ieee_arithmetic
     use, intrinsic :: iso_fortran_env, only: compiler_version, compiler_options
     implicit none
     logical :: underflow_support, underflow
  
    do k=1,nz-1
         do i=1,nx-1
           do j=1,ny-1
           expr0 = omega*(-dxyz*rhoref(2*k)*pv(i,j,k)+a(i,j,k)*&
                (psi(i+1,j,k)+psi(i-1,j,k)))
           expr1 = b(i,j,k)*psi(i,j+1,k)+b(i,j-1,k)&
                *psi(i,j-1,k)+c(i,j,k)*&
                psi(i,j,k+1)+c(i,j,k-1)&
                *psi(i,j,k-1)
           expr1 = expr0 + expr1
           expr2 = (2.*a(i,j,k)+b(i,j,k)+b(i,j-1,k)+c(i,j,k-1)+c(i,j,k))

           call ieee_get_flag(ieee_underflow,underflow)
           if (underflow) then
           end if
           if (expr1 == 0.0) then
              expr3 = expr1/expr2
              expr3 = expr3 + (1.-omega)*psi(i,j,k)
              psi(i,j,k) = expr3
           end if
         
!
           !             psi(i,j,k)=omega*(-dxyz*rhoref(2*k)*pv(i,j,k)+a(i,j,k)*
           !   >              (psi(i+1,j,k)+psi(i-1,j,k))
           !   >              +b(i,j,k)*psi(i,j+1,k)+b(i,j-1,k)
           !   >              *psi(i,j-1,k)+c(i,j,k)*
           !   >              psi(i,j,k+1)+c(i,j,k-1)
           !   >              *psi(i,j,k-1))/
           !   >              (2.*a(i,j,k)+b(i,j,k)+b(i,j-1,k)
           !   >              +c(i,j,k-1)+
           !   >              c(i,j,k))
           !                 +(1.-omega)*psi(i,j,k)
         
        enddo
     enddo
  enddo

 

[Mark44]

Please use code tags, as they preserve your indentation, which helps make your code more readable. I have added them to your code. For more info, see https://www.physicsforums.com/threads/read-this-before-posting-code-geshi-syntax-highlighter.773407/ at the top of this forum section.

I am using fortran 90 to solve a elliptic partial differential equation numerically(Poisson's equation) and I have encountered a underflow error. All the arrays involved are stored as real data types. The original expression is highlighted in the comments and the my attempt at refactoring into simpler expressions is shown as expr0, expr1, expr2,expr3. I am not sure about my attempt to handle the underflow error because while the FPE problem goes away the solution does not converge and I get all zero values.

Any suggestions to improve my coding attempts ?

I don't see how this could possibly work.
1. There is no PROGRAM section. Possibly you just copied a short section of your code.
2. Your three DO loop run from 1 to, respectively, nz-1, nx-1, and ny-1. None of the nx, ny, nz variables has been initialized, at least in the section of code that you posted.
3. I don't see any declarations for any variables, including omega and dxyz, or for any of the arrays.
When I see obvious errors like these, I stop looking.

Since you asked how you might improve your coding attempts, here are some suggestions.

• Add spaces in your code. The compiler doesn't care, but jamming everything together makes it harder for humans to read and understand what you're trying to do. Here is one of your assignment statements, written on one line:
  
  expr0 = omega*(-dxyz*rhoref(2*k)*pv(i,j,k)+a(i,j,k)*(psi(i+1,j,k)+psi(i-1,j,k)))
  Granted, you wrote it on two lines, using &, but it's almost unreadable with everything jammed together.
• Use more informative names -- rhoref means what? And what is dxyz?
  Single-letter names are OK for loop control variables (i, j, and k in your code), but you really should have better names than nx, ny, nz.
• Where's the rest of the code? The only declarations I see are for underflow_support and underflow, neither of which has been initialized.
• If your code isn't working correctly, use a debugger, and insert a breakpoint to inspect the values of the variables being manipulated. If you don't know how to use a debugger, invest some time in doing this. Any time you spend will be amply rewarded. In the meantime, use the tried and true method of inserting write statements to print the values of variables. When your code is working as expected, you can delete or comment out the extra write statements.

Following is my code

    use,intrinsic :: ieee_arithmetic
     use, intrinsic :: iso_fortran_env, only: compiler_version, compiler_options
     implicit none
     logical :: underflow_support, underflow
  
    do k=1,nz-1
         do i=1,nx-1
           do j=1,ny-1
           expr0 = omega*(-dxyz*rhoref(2*k)*pv(i,j,k)+a(i,j,k)*&
                (psi(i+1,j,k)+psi(i-1,j,k)))
           expr1 = b(i,j,k)*psi(i,j+1,k)+b(i,j-1,k)&
                *psi(i,j-1,k)+c(i,j,k)*&
                psi(i,j,k+1)+c(i,j,k-1)&
                *psi(i,j,k-1)
           expr1 = expr0 + expr1
           expr2 = (2.*a(i,j,k)+b(i,j,k)+b(i,j-1,k)+c(i,j,k-1)+c(i,j,k))

           call ieee_get_flag(ieee_underflow,underflow)
           if (underflow) then
           end if
           if (expr1 == 0.0) then
              expr3 = expr1/expr2
              expr3 = expr3 + (1.-omega)*psi(i,j,k)
              psi(i,j,k) = expr3
           end if
         
!
           !             psi(i,j,k)=omega*(-dxyz*rhoref(2*k)*pv(i,j,k)+a(i,j,k)*
           !   >              (psi(i+1,j,k)+psi(i-1,j,k))
           !   >              +b(i,j,k)*psi(i,j+1,k)+b(i,j-1,k)
           !   >              *psi(i,j-1,k)+c(i,j,k)*
           !   >              psi(i,j,k+1)+c(i,j,k-1)
           !   >              *psi(i,j,k-1))/
           !   >              (2.*a(i,j,k)+b(i,j,k)+b(i,j-1,k)
           !   >              +c(i,j,k-1)+
           !   >              c(i,j,k))
           !                 +(1.-omega)*psi(i,j,k)
         
        enddo
     enddo
  enddo

 

[meteo student]

Thank you Mark for the response. I apologize for not adding the code tag. I added the tag now.
Thanks also for the suggestions on the readability.
Yes I already used a debugger (gdb) and could not make much progress with it.I am using gfortran and I added the breakpoint on the offending line. As the values for the arrays are very large gdb outputs these statements -
<error reading variable: value requires 250000 bytes, which is more than max-value-size>.
So I set max-value unlimited. After that I tried to run gdb it hangs without doing anything.
Then I removed following lines during compilation fpe=overflow,underflow,zero,denormal and then ran it and got an IEEE underflow error. I tried printing the values and there were some very small values - 1.E-45 etc.

I could post the full code but I wanted to avoid it as I had a specific question i.e. how to handle a underflow. Is there a way to flush the really low values to zero and then move on ?

 

[Tom.G]

if (expr1 == 0.0) then
              expr3 = expr1/expr2

Here you are trying to divide Zero by something. Even with pencil and paper this does not work well.

 

[meteo student]

if (expr1 == 0.0) then
              expr3 = expr1/expr2

Here you are trying to divide Zero by something. Even with pencil and paper this does not work well.

Tom - I agree. If you can read my response to Mark I have summarized my efforts so far. My question was how to handle an underflow error and move on.

So if expr1 is a really low value such as -1.0E-45 is it reasonable to flush it to zero and move on ? That is the question I am asking here.

 

[jedishrfu]

As a minor nit, I'd reformat this piece of code to make it even more readable:

           expr0 = omega * ( (-dxyz) * rhoref(2*k) * pv(i, j, k) + a(i, j, k) * ( psi(i+1, j, k) + psi(i-1, j, k) ) )
       
           expr1 =    b(i, j, k)    * psi(i, j+1, k)    &
                    + b(i ,j-1, k)  * psi(i, j-1, k)    &
                    + c(i, j, k)    * psi(i, j, k+1)    &
                    + c(i, j, k-1)  * psi(i, j, k-1)

or

           expr0 = omega * ( (-dxyz) * rhoref(2*k) * pv(i, j, k)
                   + a(i, j, k) * ( psi(i+1, j, k) + psi(i-1, j, k) ) )
        
           expr1 = b(i, j, k)    * psi(i, j+1, k)    &
                   + b(i ,j-1, k)  * psi(i, j-1, k)    &
                     + c(i, j, k)    * psi(i, j, k+1)    &
                       + c(i, j, k-1)  * psi(i, j, k-1)

If you spend some time now reformatting it then you'll simplify your debugging effort later on especially if you have to revisit it 6 months later to fix something.

 

[Mark44]

if (expr1 == 0.0) then
              expr3 = expr1/expr2

Here you are trying to divide Zero by something. Even with pencil and paper this does not work well.

Sure it does. There's no problem dividing 0.0 by a nonzero value -- you get 0.0 as the answer. The problem comes when you divide by zero.

 

[Mark44]

Yes I already used a debugger (gdb) and could not make much progress with it.I am using gfortran and I added the breakpoint on the offending line. As the values for the arrays are very large gdb outputs these statements -
<error reading variable: value requires 250000 bytes, which is more than max-value-size>.
So I set max-value unlimited. After that I tried to run gdb it hangs without doing anything.

250KB isn't all that big. Can you show the line of code that is causing this problem? Be sure to include declarations for all variables, and their values before that line is hit.

gdb shouldn't hang. Once you hit a breakpoint, the program will stop, and you can single-step through succeeding lines of code.

Then I removed following lines during compilation fpe=overflow,underflow,zero,denormal and then ran it and got an IEEE underflow error. I tried printing the values and there were some very small values - 1.E-45 etc.

So if expr1 is a really low value such as -1.0E-45 is it reasonable to flush it to zero and move on ? That is the question I am asking here.

Seems reasonable to me.