Skip to content
Snippets Groups Projects
Commit 942a2f79 authored by Dave Whipp's avatar Dave Whipp
Browse files

Added extra output for Cholesky factorization errors. It should be commented out in the future

parent e1c451ad
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/solve_with_pwssmp.f90
+ 64
− 20
View file @ 942a2f79
......@@ -64,7 +64,7 @@ double precision aux, diag
integer mrp,naux
double precision elem_dz
integer badnode,badproc,cnt,idof,rbadnode
integer badnode,badproc,cnt,idof,rbadnode,badindex
integer,dimension(:),allocatable :: bad_elems
!==============================================================================!
......@@ -134,6 +134,14 @@ iparm(1)=0
iparm(2)=0
iparm(3)=0
! Enable scaling of the WSMP matrices for LL^T factorization
if (params%wsmp_scaling) then
iparm(10)=1
endif
! Skip ordering
iparm(16)=-1
!iparm(20)=2
call pwssmp(n_iproc,ia,ja,avals,diag,perm,invp,b,ldb,nrhs,aux,naux,mrp,iparm,dparm)
......@@ -146,6 +154,7 @@ if (iproc.eq.0) then
call stop_run ('WSMP pb$')
else
write(*,'(a,f15.7,a)') shift//'wsmp: initialisation time: ',time2-time1,' s'
write(*,'(a,e15.7)') shift//'wsmp: signularity threshold (dparm(10)): ',dparm(10)
end if
end if
......@@ -192,6 +201,7 @@ if (iproc.eq.0) then
stop
else
write(*,'(a,f15.7,a)') shift//'wsmp: symbolic fact. time: ',time2-time1,' s'
write(*,'(a,e15.7,a,e15.7)') shift//'wsmp: range of avals: ',minval(avals),' - ',maxval(avals)
end if
end if
......@@ -212,10 +222,13 @@ call pwssmp(n_iproc,ia,ja,avals,diag,perm,invp,b,ldb,nrhs,aux,naux,mrp,iparm,dpa
time2=mpi_wtime()
write(*,'(a,e15.7)') shift//'wsmp: avals min on main diagonal (dparm(5)): ',dparm(5)
write(*,'(a,e15.7)') shift//'wsmp: avals max on main diagonal (dparm(4)): ',dparm(4)
if (iproc.eq.0) then
if (iparm(64) .gt. 0) then
call write_streepjes(6,1)
write (*,*) 'Cholesky factorization problem at location ',iparm(64),' in WSMP matrix'
write (*,*) 'Cholesky factorization problem at location ',iparm(64),' in WSMP matrix (iparm(64))'
call write_streepjes(6,1)
elseif (iparm(64).lt.0) then
call write_streepjes(6,1)
......@@ -228,16 +241,47 @@ badproc=-1
if (iparm(64).gt.0) then
if (iparm(64).ge.n_iproc_st .and. iparm(64).le.n_iproc_end) then
badproc=iproc
write (*,'(a36,i12,a1)') 'WSMP failed on process ',badproc,'.'
write (*,'(a36,i15,a1)') 'WSMP failed on MPI process ',badproc,'.'
endif
! List values +/- 50 rows/columns around bad value
if (iproc == 0) write(*,'(a,i15)') '+/- 50 values surrounding bad index ',iparm(64)
do i=iparm(64)-50,iparm(64)+50
if (i.ge.n_iproc_st .and. i.le.n_iproc_end) badproc=iproc
if (iproc.eq.badproc) then
write (*,'(a,i15,a,e15.7,a,i15)') 'avals at row/col ',i,': ',avals(ia(i-n_iproc_st+1)),' from process ',badproc
endif
enddo
if (iparm(64).ge.n_iproc_st .and. iparm(64).le.n_iproc_end) badproc=iproc
if (iproc.eq.badproc) then
if (params%debug.gt.0) then
write (*,'(a36,i12,a1)') 'Debug output from process ',badproc,':'
write (*,'(a36,e12.5)') 'Bad main diagonal value (avals): ',avals(ia(iparm(64)-n_iproc_st+1))
idof=mod(iparm(64),ndof)
write (*,'(a36,i15,a1)') 'Debug output from process ',badproc,':'
write (*,'(a36,e15.7)') 'Bad main diagonal value (avals): ',avals(ia(iparm(64)-n_iproc_st+1))
!write (*,'(a36,e15.7)') 'Permuted bad main diagonal value (avals): ',avals(ia(iparm(64)-n_iproc_st+1))
!write (*,'(a36,e15.7)') 'Invp bad main diagonal value (avals): ',avals(ia(iparm(64)-n_iproc_st+1))
! Do inverse permutation to find original location of bad value in octree
!badindex=invp(iparm(64))
badindex=iparm(64)
idof=mod(badindex,ndof)
if (idof.eq.0) idof=3
rbadnode = ((iparm(64)-idof)/ndof)+1
rbadnode = ((badindex-idof)/ndof)+1
badnode = vo%rtf(rbadnode)
write (*,'(a36,i15)') 'Node in DOUAR octree with bad value: ',badnode
allocate (bad_elems(8))
cnt=0
do i=1,osolve%nleaves
......@@ -251,20 +295,20 @@ if (iparm(64).gt.0) then
call write_streepjes(6,1)
do i=1,cnt
write (*,*) ''
write (*,'(a,i12,a)') '-- Output for element ',bad_elems(i),' --'
write (*,'(a32,e12.5)') 'Pressure: ',osolve%pressure(bad_elems(i))
write (*,'(a32,e12.5)') 'Smoothed pressure: ',osolve%spressure(bad_elems(i))
write (*,'(a32,e12.5)') 'Effective viscosity: ',osolve%eviscosity(bad_elems(i))
write (*,'(a32,e12.5)') 'Effective strain rate: ',osolve%e2d(bad_elems(i))
write (*,'(a,i15,a)') '-- Output for element ',bad_elems(i),' --'
write (*,'(a32,e15.7)') 'Pressure: ',osolve%pressure(bad_elems(i))
write (*,'(a32,e15.7)') 'Smoothed pressure: ',osolve%spressure(bad_elems(i))
write (*,'(a32,e15.7)') 'Effective viscosity: ',osolve%eviscosity(bad_elems(i))
write (*,'(a32,e15.7)') 'Effective strain rate: ',osolve%e2d(bad_elems(i))
elem_dz=(osolve%z(osolve%icon(5,bad_elems(i)))-osolve%z(osolve%icon(1,bad_elems(i))))*params%vex
write (*,'(a32,e12.5)') 'Surface damping contribution: ',params%damp_factor*(params%dt*(-weightel(bad_elems(i))))/(2.d0*elem_dz)
write (*,'(a32,e15.7)') 'Surface damping contribution: ',params%damp_factor*(params%dt*(-weightel(bad_elems(i))))/(2.d0*elem_dz)
write (*,'(a)') '-- Nodal output --'
write (*,'(a16,8i14)') 'Node number: ',osolve%icon(1,bad_elems(i)), &
write (*,'(a16,8i15)') 'Node number: ',osolve%icon(1,bad_elems(i)), &
osolve%icon(2,bad_elems(i)),osolve%icon(3,bad_elems(i)), &
osolve%icon(4,bad_elems(i)),osolve%icon(5,bad_elems(i)), &
osolve%icon(6,bad_elems(i)),osolve%icon(7,bad_elems(i)), &
osolve%icon(8,bad_elems(i))
write (*,'(a16,8e14.5)') 'x-coordinate: ', &
write (*,'(a16,8e15.7)') 'x-coordinate: ', &
osolve%x(osolve%icon(1,bad_elems(i))), &
osolve%x(osolve%icon(2,bad_elems(i))), &
osolve%x(osolve%icon(3,bad_elems(i))), &
......@@ -273,7 +317,7 @@ if (iparm(64).gt.0) then
osolve%x(osolve%icon(6,bad_elems(i))), &
osolve%x(osolve%icon(7,bad_elems(i))), &
osolve%x(osolve%icon(8,bad_elems(i)))
write (*,'(a16,8e14.5)') 'y-coordinate: ', &
write (*,'(a16,8e15.7)') 'y-coordinate: ', &
osolve%y(osolve%icon(1,bad_elems(i))), &
osolve%y(osolve%icon(2,bad_elems(i))), &
osolve%y(osolve%icon(3,bad_elems(i))), &
......@@ -282,7 +326,7 @@ if (iparm(64).gt.0) then
osolve%y(osolve%icon(6,bad_elems(i))), &
osolve%y(osolve%icon(7,bad_elems(i))), &
osolve%y(osolve%icon(8,bad_elems(i)))
write (*,'(a16,8e14.5)') 'z-coordinate: ', &
write (*,'(a16,8e15.7)') 'z-coordinate: ', &
osolve%z(osolve%icon(1,bad_elems(i))), &
osolve%z(osolve%icon(2,bad_elems(i))), &
osolve%z(osolve%icon(3,bad_elems(i))), &
......@@ -291,7 +335,7 @@ if (iparm(64).gt.0) then
osolve%z(osolve%icon(6,bad_elems(i))), &
osolve%z(osolve%icon(7,bad_elems(i))), &
osolve%z(osolve%icon(8,bad_elems(i)))
write (*,'(a16,8e14.5)') 'u-velocity: ', &
write (*,'(a16,8e15.7)') 'u-velocity: ', &
osolve%u(osolve%icon(1,bad_elems(i))), &
osolve%u(osolve%icon(2,bad_elems(i))), &
osolve%u(osolve%icon(3,bad_elems(i))), &
......@@ -300,7 +344,7 @@ if (iparm(64).gt.0) then
osolve%u(osolve%icon(6,bad_elems(i))), &
osolve%u(osolve%icon(7,bad_elems(i))), &
osolve%u(osolve%icon(8,bad_elems(i)))
write (*,'(a16,8e14.5)') 'v-velocity: ', &
write (*,'(a16,8e15.7)') 'v-velocity: ', &
osolve%v(osolve%icon(1,bad_elems(i))), &
osolve%v(osolve%icon(2,bad_elems(i))), &
osolve%v(osolve%icon(3,bad_elems(i))), &
......@@ -309,7 +353,7 @@ if (iparm(64).gt.0) then
osolve%v(osolve%icon(6,bad_elems(i))), &
osolve%v(osolve%icon(7,bad_elems(i))), &
osolve%v(osolve%icon(8,bad_elems(i)))
write (*,'(a16,8e14.5)') 'w-velocity: ', &
write (*,'(a16,8e15.7)') 'w-velocity: ', &
osolve%w(osolve%icon(1,bad_elems(i))), &
osolve%w(osolve%icon(2,bad_elems(i))), &
osolve%w(osolve%icon(3,bad_elems(i))), &
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment