<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Thu, Apr 30, 2015 at 5:35 AM, Danyang Su <span dir="ltr"><<a href="mailto:danyang.su@gmail.com" target="_blank">danyang.su@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">On 15-04-29 12:19 PM, Barry Smith wrote:<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
Ok, your code seems fine in terms of logic. But the vector x_vec_gbl i.e. b_react should not be in a read only state. Are you somewhere calling a VecGetArray() and not calling the VecRestoreArray()?<br>
<br>
Barry<br>
</blockquote>
I just made a triple check on VecGetArrayF90(). All the VecGetArrayF90() come with VecRestoreArrayF90(). No VecGetArray() is used in my codes.<br></blockquote><div><br></div><div>So you can check the lock status of the Vec:</div><div><br></div><div> <a href="http://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/Vec/VecLockGet.html">http://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/Vec/VecLockGet.html</a></div><div><br></div><div>Put these in your code until you can see who is locking the Vec before the call.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
Danyang<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
<br>
<br>
<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
On Apr 29, 2015, at 1:52 PM, Danyang Su <<a href="mailto:danyang.su@gmail.com" target="_blank">danyang.su@gmail.com</a>> wrote:<br>
<br>
On 15-04-29 11:30 AM, Barry Smith wrote:<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
On Apr 29, 2015, at 12:15 PM, Danyang Su <<a href="mailto:danyang.su@gmail.com" target="_blank">danyang.su@gmail.com</a>><br>
wrote:<br>
<br>
On 15-04-28 06:50 PM, Barry Smith wrote:<br>
<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
We started enforcing more checks on writing to vectors that you should not write to. Where are you calling DMLocalToGlobalBegin() ? It looks like you are trying to copy values into a global array that you should not because it is a read only input to a function, for example it is the right hand side of a linear system or the input into a SNESFormFunction(). So check the output to that call.<br>
<br>
Barry<br>
<br>
<br>
</blockquote>
This is used when local RHS vec is calculated and copied to global RHS.<br>
<br>
</blockquote>
Huhh, why are you copying anything into the RHS? Is this before you call the linear system solve? Send the code that calls the DMLocalToGlobalBegin()<br>
</blockquote>
Yes. This is called before the linear system solver. At present, we use PETSc KSP solver.<br>
The code is first developed in totally sequential version and then ported to the parallel version. So the subdomain has its own space for jacobi matrix and rhs. We do not use global RHS vector directly to set the RHS at this moment. All the RHS values of subdomain are copied to the global RHS.<br>
<br>
!> This is the function where RHS values of subdomain are copied to the global RHS values.<br>
subroutine compute_function(rank,da,x_array_loc,x_vec_loc, &<br>
x_vec_gbl,nngl,row_idx_l2pg,col_idx_l2pg, &<br>
b_non_interlaced)<br>
implicit none<br>
#include <finclude/petscsys.h><br>
#include <finclude/petscis.h><br>
#include <finclude/petscvec.h><br>
#include <finclude/petscvec.h90><br>
#include <finclude/petscmat.h><br>
#include <finclude/petscpc.h><br>
#include <finclude/petscksp.h><br>
#include <finclude/petscdmda.h><br>
#include <finclude/petscdmda.h90><br>
<br>
PetscInt :: rank<br>
DM :: da<br>
PetscReal, allocatable :: x_array_loc(:)<br>
Vec :: x_vec_loc<br>
Vec :: x_vec_gbl<br>
PetscInt :: nngl<br>
PetscInt, allocatable :: row_idx_l2pg(:)<br>
PetscInt, allocatable :: col_idx_l2pg(:)<br>
PetscBool :: b_non_interlaced<br>
<br>
PetscInt :: info_debug<br>
PetscInt :: i, j<br>
PetscErrorCode :: ierr<br>
PetscScalar, pointer :: vecpointer(:)<br>
!Zero entries<br>
call VecZeroEntries(x_vec_loc, ierr)<br>
!Get a pointer to vector data when you need access to the array<br>
call VecGetArrayF90(x_vec_loc, vecpointer, ierr)<br>
!Compute the function over the locally owned part of the grid<br>
if(b_non_interlaced) then<br>
j = nngl/2<br>
do i = 1, j<br>
vecpointer(2*i-1) = x_array_loc(i)<br>
vecpointer(2*i) = x_array_loc(i+j)<br>
end do<br>
else<br>
do i = 1, nngl<br>
vecpointer(i) = x_array_loc(i)<br>
end do<br>
end if<br>
!Restore the vector when you no longer need access to the array<br>
call VecRestoreArrayF90(x_vec_loc,vecpointer,ierr)<br>
!Insert values into global vector<br>
call DMLocalToGlobalBegin(da,x_vec_loc,INSERT_VALUES, &<br>
x_vec_gbl,ierr)<br>
!By placing code between these two statements, computations can be<br>
!done while messages are in transition.<br>
call DMLocalToGlobalEnd(da,x_vec_loc,INSERT_VALUES, &<br>
x_vec_gbl,ierr)<br>
return<br>
end subroutine<br>
<br>
<br>
<br>
!> This is the function where reactive transport equations are solved.<br>
subroutine solver_dd_snes_solve_react(ilog,idetail,a_in,b_in, &<br>
x_inout,ia_in,ja_in,nngl_in,itsolv, &<br>
over_flow,rnorm,row_idx_l2pg,col_idx_l2pg, &<br>
b_non_interlaced)<br>
use gen, only : rank, node_idx_l2lg, ittot_rt, &<br>
b_output_matrix, b_enable_output<br>
use solver_snes_function, only : form_initial_guess, &<br>
compute_function, &<br>
compute_jacobian<br>
use petsc_mpi_common, only : petsc_mpi_finalize<br>
implicit none<br>
#include <finclude/petscsys.h><br>
#include <finclude/petscvec.h><br>
#include <finclude/petscvec.h90><br>
#include <finclude/petscdmda.h><br>
<br>
PetscInt :: ilog<br>
PetscInt :: idetail<br>
PetscInt :: nngl_in<br>
PetscReal, allocatable :: a_in(:)<br>
PetscReal, allocatable :: b_in(:)<br>
PetscReal, allocatable :: x_inout(:)<br>
PetscInt, allocatable :: ia_in(:)<br>
PetscInt, allocatable :: ja_in(:)<br>
PetscInt, allocatable :: row_idx_l2pg(:)<br>
PetscInt, allocatable :: col_idx_l2pg(:)<br>
PetscInt :: itsolv<br>
PetscBool :: over_flow<br>
PetscBool :: b_non_interlaced<br>
PetscReal :: rnorm<br>
PetscViewer :: viewer<br>
PetscScalar,pointer :: vecpointer(:)<br>
<br>
PetscErrorCode :: ierr<br>
PetscInt :: info_debug<br>
character(72) :: strinum<br>
info_debug = 0<br>
if(b_output_matrix .and. b_enable_output) then<br>
write(strinum, *) ittot_rt<br>
strinum = "_"//trim(adjustl(strinum))<br>
end if<br>
<br>
!Form initial guess, only assemble the local owned part, without ghost nodes<br>
call form_initial_guess(rank,dmda_react%da,x_inout,x_react_loc,&<br>
x_react,nngl_in, row_idx_l2pg,col_idx_l2pg, &<br>
b_non_interlaced)<br>
<br>
!Compute function, only assemble the local part, without ghost nodes<br>
call compute_function(rank,dmda_react%da,b_in,b_react_loc, &<br>
b_react,nngl_in,row_idx_l2pg,col_idx_l2pg, &<br>
b_non_interlaced)<br>
!Compute jacobian matrix, only assemble the local part, without ghost nodes<br>
call compute_jacobian(rank,dmda_react%da, &<br>
a_react,a_in,ia_in,ja_in,nngl_in, &<br>
row_idx_l2pg,col_idx_l2pg, &<br>
b_non_interlaced)<br>
<br>
! Solve a x = b, where a is the Jacobian matrix.<br>
! - First, set the KSP linear operators. Here the matrix that<br>
! defines the linear system also serves as the preconditioning<br>
! matrix.<br>
! - Then solve the Newton system.<br>
#ifdef PETSC_V3_5_X<br>
call KSPSetOperators(ksp_react,a_react,a_react,ierr)<br>
#else<br>
call KSPSetOperators(ksp_react,a_react,a_react, &<br>
SAME_NONZERO_PATTERN,ierr)<br>
#endif<br>
call KSPSetDM(ksp_react,dmda_react%da,ierr)<br>
call KSPSetDMActive(ksp_react,PETSC_FALSE,ierr)<br>
<br>
call KSPSolve(ksp_react,b_react,x_react,ierr)<br>
!Get residual norm, by default, preconditioned residual norm is calculated.<br>
if (b_mykspconverged_react .and. &<br>
.not. b_use_petsc_default_react) then<br>
rnorm = rnorm_react<br>
else<br>
call KSPGetResidualNorm(ksp_react, rnorm, ierr)<br>
end if<br>
! Scatter ghost points to local vector, using the 2-step process<br>
! DMGlobalToLocalBegin(), DMGlobalToLocalEnd().<br>
! By placing code between these two statements, computations can be<br>
! done while messages are in transition.<br>
call DMGlobalToLocalBegin(dmda_react%da,x_react,INSERT_VALUES, &<br>
x_react_loc,ierr)<br>
call DMGlobalToLocalEnd(dmda_react%da,x_react,INSERT_VALUES, &<br>
x_react_loc,ierr)<br>
call VecGetArrayF90(x_react_loc,vecpointer,ierr)<br>
x_inout = vecpointer<br>
call VecRestoreArrayF90(x_react_loc,vecpointer,ierr)<br>
return<br>
end subroutine<br>
<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
Barry<br>
<br>
<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
It is strange why this error does not occur at the first timestep.<br>
<br>
Thanks,<br>
<br>
Danyang<br>
<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
On Apr 28, 2015, at 7:30 PM, Danyang Su <<a href="mailto:danyang.su@gmail.com" target="_blank">danyang.su@gmail.com</a>><br>
wrote:<br>
<br>
Hi Barry,<br>
<br>
There seems another bug (not pretty sure) in PETSc-dev, as shown below. The case I used is similar with the one I mentioned recently. I have no problem running this case using PETSc 3.5.2. But it give out the following error using PETSc-dev, even with only one processor.<br>
<br>
Thanks,<br>
<br>
Danyang<br>
<br>
timestep: 2048 time: 4.615E+00 years delt: 1.460E-20 years iter: 1 max.sia: 0.000E+00 tol.sia: 0.000E+00<br>
[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------<br>
[0]PETSC ERROR: Object is in wrong state<br>
[0]PETSC ERROR: Vec is locked read only, argument # 1<br>
[0]PETSC ERROR: See<br>
<a href="http://www.mcs.anl.gov/petsc/documentation/faq.html" target="_blank">http://www.mcs.anl.gov/petsc/documentation/faq.html</a><br>
for trouble shooting.<br>
[0]PETSC ERROR: Petsc Development GIT revision: v3.5.3-2796-g23b6c49 GIT Date: 2015-04-27 11:32:44 -0500<br>
[0]PETSC ERROR: ../min3p_thcm on a linux-gnu-dbg named nwmop by dsu Tue Apr 28 15:56:41 2015<br>
[0]PETSC ERROR: Configure options PETSC_ARCH=linux-gnu-dbg --with-cc=gcc --with-cxx=g++ --with-fc=gfortran --download-fblaslapack --download-mpich --download-mumps --download-hypre --download-superlu_dist --download-metis --download-parmetis --download-scalapack<br>
[0]PETSC ERROR: #349 VecGetArray() line 1646 in /home/dsu/Soft/PETSc/petsc-dev/src/vec/vec/interface/rvector.c<br>
[0]PETSC ERROR: #350 VecGetArrayPair() line 434 in /home/dsu/Soft/PETSc/petsc-dev/include/petscvec.h<br>
[0]PETSC ERROR: #351 VecScatterBegin_SSToSS() line 649 in /home/dsu/Soft/PETSc/petsc-dev/src/vec/vec/utils/vscat.c<br>
[0]PETSC ERROR: #352 VecScatterBegin() line 1694 in /home/dsu/Soft/PETSc/petsc-dev/src/vec/vec/utils/vscat.c<br>
[0]PETSC ERROR: #353 DMLocalToGlobalBegin_DA() line 56 in /home/dsu/Soft/PETSc/petsc-dev/src/dm/impls/da/dagtol.c<br>
[0]PETSC ERROR: #354 DMLocalToGlobalBegin() line 1944 in /home/dsu/Soft/PETSc/petsc-dev/src/dm/interface/dm.c<br>
Reduce time step for reactive transport<br>
no further time step reduction possible<br>
Optimal time increment computed by MIN3P 2.6644769214899893E-018<br>
Minimum time increment specified by user 3.6499999999999998E-018<br>
Please, reduce the MINIMUM TIME INCREMENT<br>
stop signal in time step reduction failed<br>
<br>
On 15-04-28 11:15 AM, Barry Smith wrote:<br>
<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
I am forwarding this on to the hypre support email list; hopefully they can have some advice on how to proceed.<br>
<br>
hypre folks, this is with version hypre-2.10.0b (the previous version had the same behavior). We are getting a divide by zero in gselim() line 4363 (this happens in a time dependent problem after many successful solves with time dependent matrix). I looked at the code, below, and note that there are some checks for the diagonal of A not being zero but not for the line that causes the divide by zero; is this perhaps an oversight in the hypre code? Any advice appreciated.<br>
<br>
Thanks<br>
<br>
Barry<br>
<br>
<br>
<br>
<br>
HYPRE_Int gselim(A,x,n)<br>
HYPRE_Real *A;<br>
HYPRE_Real *x;<br>
HYPRE_Int n;<br>
{<br>
HYPRE_Int err_flag = 0;<br>
HYPRE_Int j,k,m;<br>
HYPRE_Real factor;<br>
if (n==1) /* A is 1x1 */<br>
{<br>
if (A[0] != 0.0)<br>
{<br>
x[0] = x[0]/A[0];<br>
return(err_flag);<br>
}<br>
else<br>
{<br>
err_flag = 1;<br>
return(err_flag);<br>
}<br>
}<br>
else /* A is nxn. Forward elimination */<br>
{<br>
for (k = 0; k < n-1; k++)<br>
{<br>
if (A[k*n+k] != 0.0)<br>
{<br>
for (j = k+1; j < n; j++)<br>
{<br>
if (A[j*n+k] != 0.0)<br>
{<br>
factor = A[j*n+k]/A[k*n+k];<br>
for (m = k+1; m < n; m++)<br>
{<br>
A[j*n+m] -= factor * A[k*n+m];<br>
}<br>
/* Elimination step for rhs */<br>
x[j] -= factor * x[k];<br>
}<br>
}<br>
}<br>
}<br>
/* Back Substitution */<br>
for (k = n-1; k > 0; --k)<br>
{<br>
x[k] /= A[k*n+k];<br>
for (j = 0; j < k; j++)<br>
{<br>
if (A[j*n+k] != 0.0)<br>
{<br>
x[j] -= x[k] * A[j*n+k];<br>
}<br>
}<br>
}<br>
x[0] /= A[0];<br>
return(err_flag);<br>
}<br>
}<br>
<br>
<br>
<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
On Apr 28, 2015, at 12:55 PM, Danyang Su <<a href="mailto:danyang.su@gmail.com" target="_blank">danyang.su@gmail.com</a>><br>
wrote:<br>
<br>
Hi Barry,<br>
<br>
The development version of PETSc does not help to solve my problem. It still crashed due to the same error information.<br>
<br>
As Matthew mentioned, I checked the matrix and non of the diagonal entry is zero. But the diagonal values for different rows range from -1.0d-18 to 1.0d33 and the matrix is not strict diagonal dominant. When using 1 processor, the matrix is also not strict diagonal dominant, but the diagonal values have a much smaller range and the codes can run without error. Maybe there is something wrong in the matrix for this case and I need to check the matrix first.<br>
<br>
Thanks,<br>
<br>
Danyang<br>
<br>
Program received signal SIGFPE, Arithmetic exception.<br>
0x00007f668c1299d5 in gselim (A=0x2fe3330, x=0x2210680, n=9)<br>
at par_relax.c:4363<br>
4363 x[k] /= A[k*n+k];<br>
(gdb)<br>
<br>
[1]PETSC ERROR: *** unknown floating point error occurred ***<br>
[1]PETSC ERROR: The specific exception can be determined by running in a debugger. When the<br>
[1]PETSC ERROR: debugger traps the signal, the exception can be found with fetestexcept(0x3d)<br>
[1]PETSC ERROR: where the result is a bitwise OR of the following flags:<br>
[1]PETSC ERROR: FE_INVALID=0x1 FE_DIVBYZERO=0x4 FE_OVERFLOW=0x8 FE_UNDERFLOW=0x10 FE_INEXACT=0x20<br>
[1]PETSC ERROR: Try option -start_in_debugger<br>
[1]PETSC ERROR: likely location of problem given in stack below<br>
[1]PETSC ERROR: --------------------- Stack Frames ------------------------------------<br>
[1]PETSC ERROR: Note: The EXACT line numbers in the stack are not available,<br>
[1]PETSC ERROR: INSTEAD the line number of the start of the function<br>
[1]PETSC ERROR: is given.<br>
[1]PETSC ERROR: [1] PetscDefaultFPTrap line 379 /home/dsu/Soft/PETSc/petsc-dev/src/sys/error/fp.c<br>
[1]PETSC ERROR: [1] Hypre solve line 216 /home/dsu/Soft/PETSc/petsc-dev/src/ksp/pc/impls/hypre/hypre.c<br>
[1]PETSC ERROR: [1] PCApply_HYPRE line 203 /home/dsu/Soft/PETSc/petsc-dev/src/ksp/pc/impls/hypre/hypre.c<br>
[1]PETSC ERROR: [1] PCApply line 430 /home/dsu/Soft/PETSc/petsc-dev/src/ksp/pc/interface/precon.c<br>
[1]PETSC ERROR: [1] KSP_PCApply line 234 /home/dsu/Soft/PETSc/petsc-dev/include/petsc/private/kspimpl.h<br>
[1]PETSC ERROR: [1] KSPInitialResidual line 44 /home/dsu/Soft/PETSc/petsc-dev/src/ksp/ksp/interface/itres.c<br>
[1]PETSC ERROR: [1] KSPSolve_GMRES line 224 /home/dsu/Soft/PETSc/petsc-dev/src/ksp/ksp/impls/gmres/gmres.c<br>
[1]PETSC ERROR: [1] KSPSolve line 506 /home/dsu/Soft/PETSc/petsc-dev/src/ksp/ksp/interface/itfunc.c<br>
[1]PETSC ERROR: User provided function() line 0 in Unknown file trapped floating point error<br>
<br>
<br>
<br>
On 15-04-27 11:47 AM, Barry Smith wrote:<br>
<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
Danyang,<br>
<br>
I'm glad you were finally able to get to the problematic point; sorry it took so long (I want to strangle whoever turned on catching of underflows in PETSc).<br>
<br>
The hypre folks have done a great deal of work on their relaxation code since the PETSc 3.5.3 release. There is a good chance they have fixed the divide by zero problem you are hitting here. You will need to upgrade to the development version of PETSc (that uses the latest version of hypre), here are the instructions on how to obtain it<br>
<a href="http://www.mcs.anl.gov/petsc/developers/index.html" target="_blank">http://www.mcs.anl.gov/petsc/developers/index.html</a><br>
<br>
<br>
Please let us know if this resolves the problem with hypre failing.<br>
<br>
Barry<br>
<br>
<br>
<br>
<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
On Apr 27, 2015, at 11:44 AM, Danyang Su <<a href="mailto:danyang.su@gmail.com" target="_blank">danyang.su@gmail.com</a>><br>
wrote:<br>
<br>
Hi Barry,<br>
<br>
I got the following arithemetic exception after the previous bug is fixed.<br>
<br>
Loaded symbols for /lib/x86_64-linux-gnu/libnss_files.so.2<br>
0x00007f3b23b98f20 in __nanosleep_nocancel ()<br>
at ../sysdeps/unix/syscall-template.S:81<br>
81 ../sysdeps/unix/syscall-template.S: No such file or directory.<br>
(gdb) cont<br>
Continuing.<br>
<br>
Program received signal SIGFPE, Arithmetic exception.<br>
0x00007f3b260df449 in gselim (A=0x1e4c580, x=0x1d30c40, n=9)<br>
at par_relax.c:3442<br>
3442 x[k] /= A[k*n+k];<br>
(gdb)<br>
<br>
I tried both PETSc 3.5.2 and 3.5.3 and they return the same error as shown above. For 3.5.3, i edited fp.c file and then configure and make.<br>
<br>
Thanks,<br>
<br>
Danyang<br>
<br>
On 15-04-25 07:34 PM, Danyang Su wrote:<br>
<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
Hi All,<br>
<br>
The "floating point underflow" is caused by a small value divided by a very large value. This result is forced to zero and then it does not report any underflow problem. I just rerun this bad case to see if it still get stuck later. This will take a while.<br>
<br>
Thanks for all your kindly reply,<br>
<br>
Danyang<br>
<br>
On 15-04-25 07:02 PM, Barry Smith wrote:<br>
<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
Ok, you do have<br>
<br>
#ifndef PETSC_HAVE_XMMINTRIN_H<br>
#define PETSC_HAVE_XMMINTRIN_H 1<br>
#endif<br>
<br>
so the change you made should cause it to stop trapping underflow exceptions.<br>
<br>
Now in one email you reported a FPE within hypre, then I asked you to run with -start_in_debugger to determine where it happened exactly and then you reported the FPE happened in user code (what seemed to be an underflow issue). Why is this? Can you not run it where it generated the FPE in hypre using the -start_in_debugger option?<br>
<br>
Barry<br>
<br>
Perhaps you have multiple PETSC_ARCH or multiple PETSc installs to explain why you reported two different places where the exception occurred.<br>
<br>
<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
On Apr 25, 2015, at 8:31 PM, Danyang Su <<a href="mailto:danyang.su@gmail.com" target="_blank">danyang.su@gmail.com</a>><br>
wrote:<br>
<br>
<br>
<br>
On 15-04-25 06:26 PM, Matthew Knepley wrote:<br>
<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
On Sat, Apr 25, 2015 at 8:23 PM, Danyang Su <<a href="mailto:danyang.su@gmail.com" target="_blank">danyang.su@gmail.com</a>><br>
wrote:<br>
<br>
<br>
On 15-04-25 06:03 PM, Barry Smith wrote:<br>
If this is what you got in your last run<br>
<br>
at ../../gas_advection/velocity_g.F90:1344<br>
1344 cinfrt = cinfrt_dg(i1) * diff(ic,idim) !diff is a very small value, e.g., 1.0d-316<br>
then it is still catching floating point underflow, which we do not want. This means either the change I suggested you make in the fp.c code didn't work or it actually uses a different floating point trap than that one. BTW: absurd numbers like 1.0d-316 are often a symptom of uninitialized data; could that be a problem that diff is not filled correctly for all the ic, idim you are using?<br>
<br>
This going round and round is very frustrating and a waste of time. You need to be more proactive yourself and explore the code and poke around to figure out how to solve the problem.<br>
<br>
Please email $PETSC_DIR/$PETSC_ARCH/include/petscvariables.h so I can see what FP trap is being used on your machine.<br>
<br>
Barry<br>
Do you mean $PETSC_DIR/$PETSC_ARCH/conf/petscvariables? Otherwise I cannot find this file.<br>
<br>
Its include/petscconf.h<br>
Do I need to reconfigure PETSc after changing the code you mentioned?<br>
<br>
No, but you need to rebuild.<br>
<br>
</blockquote>
Yes, I have done 'gnumake'.<br>
<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
Matt<br>
Danyang<br>
<br>
<br>
<br>
On Apr 25, 2015, at 2:24 PM, Danyang Su<br>
<<a href="mailto:danyang.su@gmail.com" target="_blank">danyang.su@gmail.com</a>><br>
wrote:<br>
<br>
<br>
<br>
On 15-04-25 11:55 AM, Barry Smith wrote:<br>
On Apr 25, 2015, at 1:51 PM, Danyang Su<br>
<<a href="mailto:danyang.su@gmail.com" target="_blank">danyang.su@gmail.com</a>><br>
<br>
wrote:<br>
<br>
<br>
<br>
On 15-04-25 11:32 AM, Barry Smith wrote:<br>
<br>
I told you this yesterday.<br>
<br>
It is probably stopping here on a harmless underflow. You need to edit the PETSc code to not worry about underflow.<br>
<br>
Edit the file /home/dsu/Soft/PETSc/petsc-3.5.2/src/sys/error/fp.c and locate<br>
<br>
#elif defined PETSC_HAVE_XMMINTRIN_H<br>
_MM_SET_EXCEPTION_MASK(_MM_MASK_INEXACT);<br>
#else<br>
<br>
change it to<br>
<br>
#elif defined PETSC_HAVE_XMMINTRIN_H<br>
_MM_SET_EXCEPTION_MASK(_MM_MASK_INEXACT | _MM_MASK_UNDERFLOW);<br>
#else<br>
<br>
Then run make gnumake in the PETSc directory to compile the new version. Now link and run the program again with -fp_trap and see where it gets stuck this time.<br>
<br>
Did you do this?<br>
<br>
Barry<br>
<br>
Yes, I did change the code in fp.c and run 'make gnumake' in the PETSc directory. I just did a double check and ran make gnumake again and got the following information this time.<br>
<br>
<br>
<br>
dsu@nwmop:~/Soft/PETSc/petsc-3.5.2$<br>
<br>
make gnumake<br>
Building PETSc using GNU Make with 10 build threads<br>
==========================================<br>
make[1]: Entering directory `/home/dsu/Soft/PETSc/petsc-3.5.2'<br>
make[1]: Nothing to be done for `all'.<br>
make[1]: Leaving directory `/home/dsu/Soft/PETSc/petsc-3.5.2'<br>
=========================================<br>
<br>
<br>
Then I recompiled the codes, ran with -fp_trap and still got the following error<br>
<br>
Backtrace for this error:<br>
Note: The EXACT line numbers in the stack are not available,<br>
[2]PETSC ERROR: INSTEAD the line number of the start of the function<br>
[2]PETSC ERROR: is given.<br>
[2]PETSC ERROR: [2] PetscDefaultFPTrap line 379 /home/dsu/Soft/PETSc/petsc-3.5.2/src/sys/error/fp.c<br>
INSTEAD the line number of the start of the function<br>
[3]PETSC ERROR: is given.<br>
[3]PETSC ERROR: [3] PetscDefaultFPTrap line 379 /home/dsu/Soft/PETSc/petsc-3.5.2/src/sys/error/fp.c<br>
[2]PETSC ERROR: User provided function() line 0 in Unknown file trapped floating point error<br>
[3]PETSC ERROR: User provided function() line 0 in Unknown file trapped floating point error<br>
<br>
<br>
This is different then what you sent a few minutes ago where it crashed in hypre.<br>
<br>
Anyways you need to use the -start_in_debugger business I sent in the previous email to see the exact place the problem occurs.<br>
<br>
Here is the information shown on gdb screen<br>
<br>
Program received signal SIGFPE, Arithmetic exception.<br>
0x00000000006c2bef in velocity_g (l_sufx=1, suffix=..., nmax=12, njamxc=34,<br>
cinfradx=..., radial_coordx=.FALSE., _suffix=3)<br>
at ../../gas_advection/velocity_g.F90:1344<br>
1344 cinfrt = cinfrt_dg(i1) * diff(ic,idim) !diff is a very small value, e.g., 1.0d-316<br>
(gdb)<br>
<br>
After type cont on gdb screen, I got error information as below<br>
<br>
[1]PETSC ERROR: *** unknown floating point error occurred ***<br>
[1]PETSC ERROR: The specific exception can be determined by running in a debugger. When the<br>
[1]PETSC ERROR: debugger traps the signal, the exception can be found with fetestexcept(0x3d)<br>
[1]PETSC ERROR: where the result is a bitwise OR of the following flags:<br>
[1]PETSC ERROR: FE_INVALID=0x1 FE_DIVBYZERO=0x4 FE_OVERFLOW=0x8 FE_UNDERFLOW=0x10 FE_INEXACT=0x20<br>
[1]PETSC ERROR: Try option -start_in_debugger<br>
[1]PETSC ERROR: likely location of problem given in stack below<br>
[1]PETSC ERROR: --------------------- Stack Frames ------------------------------------<br>
[1]PETSC ERROR: Note: The EXACT line numbers in the stack are not available,<br>
[1]PETSC ERROR: INSTEAD the line number of the start of the function<br>
[1]PETSC ERROR: is given.<br>
[1]PETSC ERROR: [1] PetscDefaultFPTrap line 379 /home/dsu/Soft/PETSc/petsc-3.5.2/src/sys/error/fp.c<br>
[1]PETSC ERROR: User provided function() line 0 in Unknown file trapped floating point error<br>
[0]PETSC ERROR: *** unknown floating point error occurred ***<br>
[0]PETSC ERROR: The specific exception can be determined by running in a debugger. When the<br>
[0]PETSC ERROR: debugger traps the signal, the exception can be found with fetestexcept(0x3d)<br>
[0]PETSC ERROR: where the result is a bitwise OR of the following flags:<br>
[0]PETSC ERROR: FE_INVALID=0x1 FE_DIVBYZERO=0x4 FE_OVERFLOW=0x8 FE_UNDERFLOW=0x10 FE_INEXACT=0x20<br>
[0]PETSC ERROR: Try option -start_in_debugger<br>
[0]PETSC ERROR: likely location of problem given in stack below<br>
<br>
Thanks,<br>
<br>
Danyang<br>
Thanks,<br>
<br>
Danyang<br>
<br>
On Apr 25, 2015, at 1:05 AM, Danyang Su<br>
<<a href="mailto:danyang.su@gmail.com" target="_blank">danyang.su@gmail.com</a>><br>
<br>
wrote:<br>
<br>
Hi Barry and Satish,<br>
<br>
How can I get rid of unknown floating point error when a very small value is multiplied.<br>
<br>
e.g.,<br>
cinfrt_dg(i1) and diff(ic,idim) are 1.0250235986806329E-008 8.6178408169776945E-317 respectively,<br>
<br>
cinfrt = cinfrt_dg(i1) * diff(ic,idim)<br>
<br>
I get the following error when run with "-fp_trap -start_in_debugger".<br>
<br>
Backtrace for this error:<br>
*** unknown floating point error occurred ***<br>
[2]PETSC ERROR: The specific exception can be determined by running in a debugger. When the<br>
[2]PETSC ERROR: debugger traps the signal, the exception can be found with fetestexcept(0x3d)<br>
[2]PETSC ERROR: cinfrt_dg(i1),diff(ic,idim) 1.0250235986806329E-008 8.6178408169776945E-317<br>
where the result is a bitwise OR of the following flags:<br>
[2]PETSC ERROR: FE_INVALID=0x1 FE_DIVBYZERO=0x4 FE_OVERFLOW=0x8 FE_UNDERFLOW=0x10 FE_INEXACT=0x20<br>
[2]PETSC ERROR: Try option -start_in_debugger<br>
[2]PETSC ERROR: likely location of problem given in stack below<br>
<br>
Thanks,<br>
<br>
Danyang<br>
<br>
On 15-04-24 01:54 PM, Danyang Su wrote:<br>
<br>
On 15-04-24 01:23 PM, Satish Balay wrote:<br>
<br>
c 4 1.0976214263087059E-067<br>
<br>
I don't think this number can be stored in a real*4.<br>
<br>
Satish<br>
<br>
Thanks, Satish. It is caused by this number.<br>
<br>
On Fri, 24 Apr 2015, Danyang Su wrote:<br>
<br>
<br>
On 15-04-24 11:12 AM, Barry Smith wrote:<br>
<br>
On Apr 24, 2015, at 1:05 PM, Danyang Su<br>
<<a href="mailto:danyang.su@gmail.com" target="_blank">danyang.su@gmail.com</a>><br>
<br>
wrote:<br>
<br>
Hi All,<br>
<br>
One of my case crashes because of floating point exception when using 4<br>
processors, as shown below. But if I run this case with 1 processor, it<br>
works fine. I have tested the codes with around 100 cases up to 768<br>
processors, all other cases work fine. I just wonder if this kind of error<br>
is caused because of NaN in jacobi matrix, RHS or preconditioner?<br>
<br>
Yes, almost for sure it is one of these places.<br>
<br>
First run the bad case with -fp_trap if all goes well you'll see the<br>
function where the FPE is generated. Then run also with -start_in_debugger<br>
and<br>
type cont in all four debugger windows. When the FPE happens the debugger<br>
should stop showing exactly where the FPE happens.<br>
<br>
Barry<br>
<br>
Hi Barry,<br>
<br>
If run with -fp_trap -start_in_debugger, I got the following error<br>
<br>
[0]PETSC ERROR: *** unknown floating point error occurred ***<br>
[0]PETSC ERROR: The specific exception can be determined by running in a<br>
debugger. When the<br>
[0]PETSC ERROR: debugger traps the signal, the exception can be found with<br>
fetestexcept(0x3d)<br>
[0]PETSC ERROR: where the result is a bitwise OR of the following flags:<br>
[0]PETSC ERROR: FE_INVALID=0x1 FE_DIVBYZERO=0x4 FE_OVERFLOW=0x8<br>
FE_UNDERFLOW=0x10 FE_INEXACT=0x20<br>
[0]PETSC ERROR: Try option -start_in_debugger<br>
[0]PETSC ERROR: likely location of problem given in stack below<br>
[0]PETSC ERROR: --------------------- Stack Frames<br>
------------------------------------<br>
[0]PETSC ERROR: Note: The EXACT line numbers in the stack are not available,<br>
[0]PETSC ERROR: INSTEAD the line number of the start of the function<br>
[0]PETSC ERROR: is given.<br>
[0]PETSC ERROR: [0] PetscDefaultFPTrap line 379<br>
/home/dsu/Soft/PETSc/petsc-3.5.2/src/sys/error/fp.c<br>
[0]PETSC ERROR: User provided function() line 0 in Unknown file trapped<br>
floating point error<br>
<br>
Program received signal SIGABRT: Process abort signal.<br>
<br>
Backtrace for this error:<br>
#0 0x7F4FEAB1C7D7<br>
#1 0x7F4FEAB1CDDE<br>
#2 0x7F4FE9E1AD3F<br>
#3 0x7F4FE9E1ACC9<br>
#4 0x7F4FE9E1E0D7<br>
#5 0x7F4FEB0B6DCB<br>
#6 0x7F4FEB0B1825<br>
#7 0x7F4FEB0B817F<br>
#8 0x7F4FE9E1AD3F<br>
#9 0x6972C8 in tprfrtlc_ at tprfrtlc.F90:2393 (discriminator 3)<br>
#10 0x4C6C87 in gcreact_ at gcreact.F90:678<br>
#11 0x707E19 in initicrt_ at initicrt.F90:589<br>
#12 0x4F42D0 in initprob_ at initprob.F90:430<br>
#13 0x5AAF72 in driver_pc at driver_pc.F90:438<br>
<br>
I checked the code at tprfrtlc.F90:2393,<br>
<br>
realbuffer_gb(1:nvars) = (/time,(c(ic),ic=1,nc-1), &<br>
(cx(ix),ix=1,nxout)/)<br>
<br>
All the values (time, c, cx) are reasonable, as shown below. The only<br>
possibility is that realbuffer_gb is in declared as real*4 if using sing<br>
precision output while time, c, cx are declared in real*8. I have a lot of<br>
similar data conversion from real*8 to real*4 output, other code does not<br>
return error.<br>
<br>
time 0.0000000000000000<br>
c 1 9.9999999999999995E-008<br>
c 2 3.1555251077549618E-003<br>
c 3 7.1657814842179362E-008<br>
c 4 1.0976214263087059E-067<br>
c 5 5.2879822292305797E-004<br>
c 6 9.9999999999999964E-005<br>
c 7 6.4055731968811337E-005<br>
c 8 3.4607572892578404E-020<br>
cx 1 3.4376650636008101E-005<br>
cx 2 7.3989678854017763E-012<br>
cx 3 9.5317170613607207E-012<br>
cx 4 2.2344525794718353E-015<br>
cx 5 3.0624685689695889E-008<br>
cx 6 1.0046157902783967E-007<br>
cx 7 1.5320169154914984E-004<br>
cx 8 8.6930292776346176E-014<br>
cx 9 3.5944267559348721E-005<br>
cx 10 3.0072645866951157E-018<br>
cx 11 2.3592486321095017E-013<br>
<br>
Thanks,<br>
<br>
Danyang<br>
<br>
<br>
I can check all the entries of jacobi matrix to see if the value is valid,<br>
but this seems not a good idea as it takes a long time to reach this<br>
point. If I restart the simulation from a specified time (e.g., 7.685 in<br>
this case), then the error does not occur.<br>
<br>
Would you please give me any suggestion on debugging this case?<br>
<br>
Thanks and Regards,<br>
<br>
Danyang<br>
<br>
<br>
timestep: 2730 time: 7.665E+00 years delt: 1.000E-02 years iter: 1<br>
timestep: max.sia: 0.000E+00 tol.sia: 0.000E+00<br>
timestep: 2731 time: 7.675E+00 years delt: 1.000E-02 years iter: 1<br>
timestep: max.sia: 0.000E+00 tol.sia: 0.000E+00<br>
timestep: 2732 time: 7.685E+00 years delt: 1.000E-02 years iter: 1<br>
timestep: max.sia: 0.000E+00 tol.sia: 0.000E+00<br>
timestep: 2733 time: 7.695E+00 years delt: 1.000E-02 years iter: 1<br>
timestep: max.sia: 0.000E+00 tol.sia: 0.000E+00<br>
timestep: 2734 time: 7.705E+00 years delt: 1.000E-02 years iter: 1<br>
timestep: max.sia: 0.000E+00 tol.sia: 0.000E+00<br>
Reduce time step for reactive transport<br>
timestep: 2734 time: 7.700E+00 years delt: 5.000E-03 years iter: 1<br>
timestep: max.sia: 0.000E+00 tol.sia: 0.000E+00<br>
Reduce time step for reactive transport<br>
timestep: 2734 time: 7.697E+00 years delt: 2.500E-03 years iter: 1<br>
timestep: max.sia: 0.000E+00 tol.sia: 0.000E+00<br>
[1]PETSC ERROR: --------------------- Error Message<br>
--------------------------------------------------------------<br>
[1]PETSC ERROR: Floating point exception<br>
[2]PETSC ERROR: --------------------- Error Message<br>
--------------------------------------------------------------<br>
[2]PETSC ERROR: Floating point exception<br>
[2]PETSC ERROR: Vec entry at local location 0 is not-a-number or infinite<br>
at end of function: Parameter number 3<br>
[2]PETSC ERROR: See<br>
<br>
<a href="http://www.mcs.anl.gov/petsc/documentation/faq.html" target="_blank">http://www.mcs.anl.gov/petsc/documentation/faq.html</a><br>
<br>
<br>
for trouble shooting.<br>
[2]PETSC ERROR: Petsc Release Version 3.5.2, Sep, 08, 2014<br>
[2]PETSC ERROR: [1]PETSC ERROR: Vec entry at local location 0 is<br>
not-a-number or infinite at end of function: Parameter number 3<br>
[1]PETSC ERROR: See<br>
<br>
<a href="http://www.mcs.anl.gov/petsc/documentation/faq.html" target="_blank">http://www.mcs.anl.gov/petsc/documentation/faq.html</a><br>
<br>
<br>
for trouble shooting.<br>
[1]PETSC ERROR: Petsc Release Version 3.5.2, Sep, 08, 2014<br>
[1]PETSC ERROR: ../min3p_thcm_petsc_dbg on a linux-gnu-dbg named nwmop by<br>
dsu Thu Apr 23 15:38:52 2015<br>
[1]PETSC ERROR: Configure options PETSC_ARCH=linux-gnu-dbg --with-cc=gcc<br>
--with-cxx=g++ --with-fc=gfortran --download-fblaslapack --download-mpich<br>
--download-mumps --download-hypre --download-superlu_dist --download-metis<br>
--download-parmetis --download-scalapack<br>
[1]PETSC ERROR: #1 VecValidValues() line 34 in<br>
/home/dsu/Soft/PETSc/petsc-3.5.2/src/vec/vec/interface/rvector.c<br>
../min3p_thcm_petsc_dbg on a linux-gnu-dbg named nwmop by dsu Thu Apr 23<br>
15:38:52 2015<br>
[2]PETSC ERROR: Configure options PETSC_ARCH=linux-gnu-dbg --with-cc=gcc<br>
--with-cxx=g++ --with-fc=gfortran --download-fblaslapack --download-mpich<br>
--download-mumps --download-hypre --download-superlu_dist --download-metis<br>
--download-parmetis --download-scalapack<br>
[2]PETSC ERROR: #1 VecValidValues() line 34 in<br>
/home/dsu/Soft/PETSc/petsc-3.5.2/src/vec/vec/interface/rvector.c<br>
[2]PETSC ERROR: [1]PETSC ERROR: #2 PCApply() line 442 in<br>
/home/dsu/Soft/PETSc/petsc-3.5.2/src/ksp/pc/interface/precon.c<br>
[1]PETSC ERROR: #2 PCApply() line 442 in<br>
/home/dsu/Soft/PETSc/petsc-3.5.2/src/ksp/pc/interface/precon.c<br>
[2]PETSC ERROR: #3 KSP_PCApply() line 230 in<br>
/home/dsu/Soft/PETSc/petsc-3.5.2/include/petsc-private/kspimpl.h<br>
#3 KSP_PCApply() line 230 in<br>
/home/dsu/Soft/PETSc/petsc-3.5.2/include/petsc-private/kspimpl.h<br>
[1]PETSC ERROR: #4 KSPInitialResidual() line 63 in<br>
/home/dsu/Soft/PETSc/petsc-3.5.2/src/ksp/ksp/interface/itres.c<br>
[2]PETSC ERROR: #4 KSPInitialResidual() line 63 in<br>
/home/dsu/Soft/PETSc/petsc-3.5.2/src/ksp/ksp/interface/itres.c<br>
[1]PETSC ERROR: #5 KSPSolve_GMRES() line 234 in<br>
/home/dsu/Soft/PETSc/petsc-3.5.2/src/ksp/ksp/impls/gmres/gmres.c<br>
[2]PETSC ERROR: #5 KSPSolve_GMRES() line 234 in<br>
/home/dsu/Soft/PETSc/petsc-3.5.2/src/ksp/ksp/impls/gmres/gmres.c<br>
[2]PETSC ERROR: #6 KSPSolve() line 459 in<br>
/home/dsu/Soft/PETSc/petsc-3.5.2/src/ksp/ksp/interface/itfunc.c<br>
[1]PETSC ERROR: #6 KSPSolve() line 459 in<br>
/home/dsu/Soft/PETSc/petsc-3.5.2/src/ksp/ksp/interface/itfunc.c<br>
^C[mpiexec@nwmop] Sending Ctrl-C to processes as requested<br>
[mpiexec@nwmop] Press Ctrl-C again to force abort<br>
<br>
<br>
<br>
<br><span class=""><font color="#888888">
<br>
-- <br>
What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener<br>
<br>
</font></span></blockquote><span class=""><font color="#888888">
<petscconf.h><br>
<br>
</font></span></blockquote></blockquote></blockquote></blockquote></blockquote></blockquote></blockquote></blockquote></blockquote></blockquote></blockquote></blockquote></blockquote>
<br>
</blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div>
</div></div>