Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Handling underflow errors in fortran 90

  1. Jul 10, 2017 #1
    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 ?
    Code (Fortran):


        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
     
    Last edited: Jul 11, 2017
  2. jcsd
  3. Jul 11, 2017 #2

    Mark44

    Staff: Mentor

    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 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:
      Code (Fortran):
      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.

     
  4. Jul 11, 2017 #3
    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 ?
     
    Last edited: Jul 11, 2017
  5. Jul 11, 2017 #4

    Tom.G

    User Avatar
    Science Advisor

    Code (Fortran):

    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.
     
  6. Jul 11, 2017 #5
    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.
     
  7. Jul 11, 2017 #6

    jedishrfu

    Staff: Mentor

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

    Code (Fortran):

               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

    Code (Fortran):

               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.
     
  8. Jul 11, 2017 #7

    Mark44

    Staff: Mentor

    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.
     
  9. Jul 11, 2017 #8

    Mark44

    Staff: Mentor

    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.
    Seems reasonable to me.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Handling underflow errors in fortran 90
Loading...