- #1

- 72

- 0

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 ?

Any suggestions to improve my coding attempts ?

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: