223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Merge branch 'less-terminallyIll' into 'development' 

removed terminallyIll as global variable

See merge request damask/DAMASK!895
This commit is contained in:
Philip Eisenlohr 2024-01-08 15:22:21 +00:00
parent fdcea09133 1f9a372aa7
commit 9c796001e7
10 changed files with 72 additions and 68 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/Marc/DAMASK_Marc.f90
View File

@ -280,7 +280,6 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
elseif (lovl == 6) then ! stress requested by marc
computationMode = materialpoint_CALCRESULTS
if (cptim > theTime .or. inc /= theInc) then ! reached "convergence"
terminallyIll = .false.
cycleCounter = -1 ! first calc step increments this to cycle = 0
if (inc == 0) then ! >> start of analysis <<
lastIncConverged = .false.
@ -299,7 +298,6 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
else if ( timinc < theDelta ) then ! >> cutBack <<
lastIncConverged = .false.
outdatedByNewInc = .false.
terminallyIll = .false.
cycleCounter = -1 ! first calc step increments this to cycle = 0
print'(a,i6,1x,i2)', '<< HYPELA2 >> cutback detected..! ',m(1),nn
end if ! convergence treatment end

15
src/Marc/materialpoint_Marc.f90
View File

@ -37,6 +37,8 @@ module materialpoint_Marc
integer, public :: &
cycleCounter = 0 !< needs description
logical, public :: &
broken = .false. !< needs description
integer, parameter, public :: &
materialpoint_CALCRESULTS = 2**0, &
@ -129,9 +131,8 @@ subroutine materialpoint_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip,
integer elCP, & ! crystal plasticity element number
i, j, k, l, m, n, ph, homog, mySource,ce
real(pREAL), parameter :: ODD_STRESS = 1e15_pREAL, & !< return value for stress if terminallyIll
ODD_JACOBIAN = 1e50_pREAL !< return value for jacobian if terminallyIll
real(pREAL), parameter :: ODD_STRESS = 1e15_pREAL, & !< return value for stress if broken
ODD_JACOBIAN = 1e50_pREAL !< return value for jacobian if broken
elCP = discretization_Marc_FEM2DAMASK_elem(elFE)
ce = discretization_Marc_FEM2DAMASK_cell(ip,elFE)
@ -149,17 +150,17 @@ subroutine materialpoint_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip,
if (iand(mode, materialpoint_CALCRESULTS) /= 0) then
validCalculation: if (terminallyIll) then
validCalculation: if (broken) then
call random_number(rnd)
if (rnd < 0.5_pREAL) rnd = rnd - 1.0_pREAL
materialpoint_cs(1:6,ip,elCP) = ODD_STRESS * rnd
materialpoint_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_eye(6)
else validCalculation
call homogenization_mechanical_response(dt,(elCP-1)*discretization_nIPs + ip, &
(elCP-1)*discretization_nIPs + ip)
call homogenization_mechanical_response(broken, &
dt,(elCP-1)*discretization_nIPs + ip, (elCP-1)*discretization_nIPs + ip)
terminalIllness: if (terminallyIll) then
terminalIllness: if (broken) then
call random_number(rnd)
if (rnd < 0.5_pREAL) rnd = rnd - 1.0_pREAL

17
src/grid/grid_mech_FEM.f90
View File

@ -83,8 +83,8 @@ module grid_mechanical_FEM
real(pREAL) :: &
err_BC !< deviation from stress BC
integer :: &
totalIter = 0 !< total iteration in current increment
integer :: totalIter = 0 !< total iteration in current increment
logical :: broken
public :: &
grid_mechanical_FEM_init, &
@ -271,7 +271,7 @@ subroutine grid_mechanical_FEM_init(num_grid)
end if restartRead
call utilities_updateCoords(F)
call utilities_constitutiveResponse(P_current,P_av,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2
call utilities_constitutiveResponse(broken,P_current,P_av,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2
F, & ! target F
0.0_pREAL) ! time increment
call DMDAVecRestoreArrayF90(DM_mech,u_PETSc,u,err_PETSc)
@ -325,8 +325,7 @@ function grid_mechanical_FEM_solution(incInfoIn) result(solution)
solution%converged = reason > 0
solution%iterationsNeeded = totalIter
solution%termIll = terminallyIll
terminallyIll = .false.
solution%termIll = broken
P_aim = merge(P_av,P_aim,params%stress_mask)
end function grid_mechanical_FEM_solution
@ -336,7 +335,7 @@ end function grid_mechanical_FEM_solution
!> @brief forwarding routine
!> @details find new boundary conditions and best F estimate for end of current timestep
!--------------------------------------------------------------------------------------------------
subroutine grid_mechanical_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,&
subroutine grid_mechanical_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining, &
deformation_BC,stress_BC,rotation_BC)
logical, intent(in) :: &
@ -493,7 +492,7 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,e
BCTol = max(maxval(abs(P_av))*num%eps_stress_rtol, num%eps_stress_atol)
if ((totalIter >= num%itmin .and. all([err_div/divTol, err_BC/BCTol] < 1.0_pREAL)) &
.or. terminallyIll) then
.or. broken) then
reason = 1
elseif (totalIter >= num%itmax) then
reason = -1
@ -571,10 +570,10 @@ subroutine formResidual(da_local,x_local, &
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
call utilities_constitutiveResponse(P_current,&
call utilities_constitutiveResponse(broken,P_current,&
P_av,C_volAvg,devNull, &
F,params%Delta_t,params%rotation_BC)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
call MPI_Allreduce(MPI_IN_PLACE,broken,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
call parallelization_chkerr(err_MPI)
!--------------------------------------------------------------------------------------------------

15
src/grid/grid_mech_spectral_basic.f90
View File

@ -83,8 +83,8 @@ module grid_mechanical_spectral_basic
err_BC, & !< deviation from stress BC
err_div !< RMS of div of P
integer :: &
totalIter = 0 !< total iteration in current increment
integer :: totalIter = 0 !< total iteration in current increment
logical :: broken
public :: &
grid_mechanical_spectral_basic_init, &
@ -228,7 +228,7 @@ subroutine grid_mechanical_spectral_basic_init(num_grid)
end if restartRead
call utilities_updateCoords(reshape(F,shape(F_lastInc)))
call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
call utilities_constitutiveResponse(broken,P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
reshape(F,shape(F_lastInc)), & ! target F
0.0_pREAL) ! time increment
call DMDAVecRestoreArrayF90(DM_mech,F_PETSc,F,err_PETSc) ! deassociate pointer
@ -287,8 +287,7 @@ function grid_mechanical_spectral_basic_solution(incInfoIn) result(solution)
solution%converged = reason > 0
solution%iterationsNeeded = totalIter
solution%termIll = terminallyIll
terminallyIll = .false.
solution%termIll = broken
P_aim = merge(P_av,P_aim,params%stress_mask)
end function grid_mechanical_spectral_basic_solution
@ -453,7 +452,7 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
BCTol = max(maxval(abs(P_av))*num%eps_stress_rtol, num%eps_stress_atol)
if ((totalIter >= num%itmin .and. all([err_div/divTol, err_BC/BCTol] < 1.0_pREAL)) &
.or. terminallyIll) then
.or. broken) then
reason = 1
elseif (totalIter >= num%itmax) then
reason = -1
@ -515,10 +514,10 @@ subroutine formResidual(residual_subdomain, F, &
end if newIteration
associate (P => r)
call utilities_constitutiveResponse(P, &
call utilities_constitutiveResponse(broken,P, &
P_av,C_volAvg,C_minMaxAvg, &
F,params%Delta_t,params%rotation_BC)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
call MPI_Allreduce(MPI_IN_PLACE,broken,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
call parallelization_chkerr(err_MPI)
err_div = utilities_divergenceRMS(P)
end associate

17
src/grid/grid_mech_spectral_polarization.f90
View File

@ -95,8 +95,8 @@ module grid_mechanical_spectral_polarization
err_curl, & !< RMS of curl of F
err_div !< RMS of div of P
integer :: &
totalIter = 0 !< total iteration in current increment
integer :: totalIter = 0 !< total iteration in current increment
logical :: broken
public :: &
grid_mechanical_spectral_polarization_init, &
@ -257,7 +257,7 @@ subroutine grid_mechanical_spectral_polarization_init(num_grid)
end if restartRead
call utilities_updateCoords(reshape(F,shape(F_lastInc)))
call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
call utilities_constitutiveResponse(broken,P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
reshape(F,shape(F_lastInc)), & ! target F
0.0_pREAL) ! time increment
call DMDAVecRestoreArrayF90(DM_mech,FandF_tau_PETSc,FandF_tau,err_PETSc) ! deassociate pointer
@ -322,8 +322,7 @@ function grid_mechanical_spectral_polarization_solution(incInfoIn) result(soluti
solution%converged = reason > 0
solution%iterationsNeeded = totalIter
solution%termIll = terminallyIll
terminallyIll = .false.
solution%termIll = broken
P_aim = merge(P_av,P_aim,params%stress_mask)
end function grid_mechanical_spectral_polarization_solution
@ -517,7 +516,7 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
BCTol = max(maxval(abs(P_av))*num%eps_stress_rtol, num%eps_stress_atol)
if ((totalIter >= num%itmin .and. all([err_div/divTol, err_curl/curlTol, err_BC/BCTol] < 1.0_pREAL)) &
.or. terminallyIll) then
.or. broken) then
reason = 1
elseif (totalIter >= num%itmax) then
reason = -1
@ -605,14 +604,14 @@ subroutine formResidual(residual_subdomain, FandF_tau, &
err_curl = utilities_curlRMS(F)
#ifdef __GFORTRAN__
call utilities_constitutiveResponse(r_F, &
call utilities_constitutiveResponse(broken,r_F, &
#else
associate (P => r_F)
call utilities_constitutiveResponse(P, &
call utilities_constitutiveResponse(broken, P, &
#endif
P_av,C_volAvg,C_minMaxAvg, &
F - r_F_tau/num%beta,params%Delta_t,params%rotation_BC)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
call MPI_Allreduce(MPI_IN_PLACE,broken,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
#ifdef __GFORTRAN__
err_div = utilities_divergenceRMS(r_F)
#else

7
src/grid/grid_mech_utilities.f90
View File

@ -113,9 +113,10 @@ end function utilities_maskedCompliance
!--------------------------------------------------------------------------------------------------
!> @brief Calculate constitutive response.
!--------------------------------------------------------------------------------------------------
subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
subroutine utilities_constitutiveResponse(broken, P,P_av,C_volAvg,C_minmaxAvg,&
F,Delta_t,rotation_BC)
logical, intent(out) :: broken
real(pREAL), intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness
real(pREAL), intent(out), dimension(3,3) :: P_av !< average PK stress
real(pREAL), intent(out), dimension(3,3,cells(1),cells(2),cells3) :: P !< PK stress
@ -123,19 +124,19 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
real(pREAL), intent(in) :: Delta_t !< loading time
type(tRotation), intent(in), optional :: rotation_BC !< rotation of load frame
integer :: i
integer(MPI_INTEGER_KIND) :: err_MPI
real(pREAL), dimension(3,3,3,3) :: dPdF_max, dPdF_min
real(pREAL) :: dPdF_norm_max, dPdF_norm_min
real(pREAL), dimension(2) :: valueAndRank !< pair of min/max norm of dPdF to synchronize min/max of dPdF
print'(/,1x,a)', '... evaluating constitutive response ......................................'
flush(IO_STDOUT)
homogenization_F = reshape(F,[3,3,product(cells(1:2))*cells3]) ! set materialpoint target F to estimated field
call homogenization_mechanical_response(Delta_t,1,product(cells(1:2))*cells3) ! calculate P field
call homogenization_mechanical_response(broken,Delta_t,1,product(cells(1:2))*cells3) ! calculate P field
P = reshape(homogenization_P, [3,3,cells(1),cells(2),cells3])
P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt

5
src/grid/grid_thermal_spectral.f90
View File

@ -56,6 +56,7 @@ module grid_thermal_spectral
integer :: totalIter = 0 !< total iteration in current increment
real(pREAL), dimension(3,3) :: K_ref
real(pREAL) :: mu_ref, Delta_t_
logical :: broken
public :: &
grid_thermal_spectral_init, &
@ -206,7 +207,7 @@ function grid_thermal_spectral_solution(Delta_t) result(solution)
call SNESGetConvergedReason(SNES_thermal,reason,err_PETSc)
CHKERRQ(err_PETSc)
solution%converged = reason > 0
solution%converged = reason > 0 .and. .not. broken
solution%iterationsNeeded = merge(totalIter,num%itmax,solution%converged)
call SNESGetDM(SNES_thermal,DM_thermal,err_PETSc)
@ -323,7 +324,7 @@ subroutine formResidual(residual_subdomain,x_scal,r,dummy,err_PETSc)
real(pREAL), dimension(3,cells(1),cells(2),cells3) :: vectorField
call homogenization_thermal_response(Delta_t_,1,product(cells(1:2))*cells3)
call homogenization_thermal_response(broken,Delta_t_,1,product(cells(1:2))*cells3)
associate(T => x_scal)
vectorField = utilities_ScalarGradient(T)

42
src/homogenization.f90
View File

@ -46,9 +46,6 @@ module homogenization
thermal_active, &
damage_active
logical, public :: &
terminallyIll = .false. !< at least one material point is terminally ill
!--------------------------------------------------------------------------------------------------
! General variables for the homogenization at a material point
real(pREAL), dimension(:,:,:), allocatable, public :: &
@ -213,8 +210,9 @@ end subroutine homogenization_init
!--------------------------------------------------------------------------------------------------
!> @brief
!--------------------------------------------------------------------------------------------------
subroutine homogenization_mechanical_response(Delta_t,cell_start,cell_end)
subroutine homogenization_mechanical_response(broken,Delta_t,cell_start,cell_end)
logical, intent(out) :: broken
real(pREAL), intent(in) :: Delta_t !< time increment
integer, intent(in) :: &
cell_start, cell_end
@ -226,8 +224,8 @@ subroutine homogenization_mechanical_response(Delta_t,cell_start,cell_end)
doneAndHappy
!$OMP PARALLEL
!$OMP DO PRIVATE(en,ho,co,converged,doneAndHappy)
broken = .false.
!$OMP PARALLEL DO PRIVATE(en,ho,co,converged,doneAndHappy)
do ce = cell_start, cell_end
en = material_entry_homogenization(ce)
@ -241,7 +239,7 @@ subroutine homogenization_mechanical_response(Delta_t,cell_start,cell_end)
doneAndHappy = [.false.,.true.]
convergenceLooping: do while (.not. (terminallyIll .or. doneAndHappy(1)))
convergenceLooping: do while (.not. (broken .or. doneAndHappy(1)))
call mechanical_partition(homogenization_F(1:3,1:3,ce),ce)
converged = all([(phase_mechanical_constitutive(Delta_t,co,ce),co=1,homogenization_Nconstituents(ho))])
@ -252,17 +250,22 @@ subroutine homogenization_mechanical_response(Delta_t,cell_start,cell_end)
doneAndHappy = [.true.,.false.]
end if
end do convergenceLooping
if (.not. converged) then
if (.not. broken) print*, ' Cell ', ce, ' failed (mechanics)'
broken = .true.
end if
converged = converged .and. all([(phase_damage_constitutive(Delta_t,co,ce),co=1,homogenization_Nconstituents(ho))])
if (.not. converged) then
if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill'
terminallyIll = .true.
if (.not. broken) print*, ' Cell ', ce, ' failed (damage)'
broken = .true.
end if
end do
!$OMP END DO
!$OMP END PARALLEL DO
!$OMP DO PRIVATE(ho)
if (broken) return
!$OMP PARALLEL DO PRIVATE(ho)
do ce = cell_start, cell_end
ho = material_ID_homogenization(ce)
do co = 1, homogenization_Nconstituents(ho)
@ -270,8 +273,7 @@ subroutine homogenization_mechanical_response(Delta_t,cell_start,cell_end)
end do
call mechanical_homogenize(Delta_t,ce)
end do
!$OMP END DO
!$OMP END PARALLEL
!$OMP END PARALLEL DO
end subroutine homogenization_mechanical_response
@ -279,8 +281,10 @@ end subroutine homogenization_mechanical_response
!--------------------------------------------------------------------------------------------------
!> @brief
!--------------------------------------------------------------------------------------------------
subroutine homogenization_thermal_response(Delta_t,cell_start,cell_end)
subroutine homogenization_thermal_response(broken, &
Delta_t,cell_start,cell_end)
logical, intent(out) :: broken
real(pREAL), intent(in) :: Delta_t !< time increment
integer, intent(in) :: &
cell_start, cell_end
@ -289,18 +293,20 @@ subroutine homogenization_thermal_response(Delta_t,cell_start,cell_end)
co, ce, ho
broken = .false.
!$OMP PARALLEL DO PRIVATE(ho)
do ce = cell_start, cell_end
if (terminallyIll) continue
if (broken) continue
ho = material_ID_homogenization(ce)
do co = 1, homogenization_Nconstituents(ho)
if (.not. phase_thermal_constitutive(Delta_t,material_ID_phase(co,ce),material_entry_phase(co,ce))) then
if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill'
terminallyIll = .true.
if (.not. broken) print*, ' Cell ', ce, ' failed (thermal)'
broken = .true.
end if
end do
end do
!$OMP END PARALLEL DO
broken = broken
end subroutine homogenization_thermal_response

5
src/mesh/FEM_utilities.f90
View File

@ -120,8 +120,9 @@ end subroutine FEM_utilities_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates constitutive response
!--------------------------------------------------------------------------------------------------
subroutine utilities_constitutiveResponse(Delta_t,P_av,forwardData)
subroutine utilities_constitutiveResponse(broken, Delta_t,P_av,forwardData)
logical, intent(out) :: broken
real(pREAL), intent(in) :: Delta_t !< loading time
logical, intent(in) :: forwardData !< age results
real(pREAL),intent(out), dimension(3,3) :: P_av !< average PK stress
@ -131,7 +132,7 @@ subroutine utilities_constitutiveResponse(Delta_t,P_av,forwardData)
print'(/,1x,a)', '... evaluating constitutive response ......................................'
call homogenization_mechanical_response(Delta_t,1,mesh_maxNips*mesh_NcpElems) ! calculate P field
call homogenization_mechanical_response(broken,Delta_t,1,mesh_maxNips*mesh_NcpElems) ! calculate P field
cutBack = .false.
P_av = sum(homogenization_P,dim=3) * wgt

13
src/mesh/mesh_mech_FEM.f90
View File

@ -68,7 +68,7 @@ module mesh_mechanical_FEM
character(len=pSTRLEN) :: incInfo
real(pREAL), dimension(3,3) :: &
P_av = 0.0_pREAL
logical :: ForwardData
logical :: ForwardData, broken
real(pREAL), parameter :: eps = 1.0e-18_pREAL
external :: & ! ToDo: write interfaces
@ -311,7 +311,7 @@ subroutine FEM_mechanical_init(mechBC,num_mesh)
call DMPlexVecSetClosure(mechanical_mesh,section,solution_local,cell,px_scal,5,err_PETSc)
CHKERRQ(err_PETSc)
end do
call utilities_constitutiveResponse(0.0_pREAL,devNull,.true.)
call utilities_constitutiveResponse(broken,0.0_pREAL,devNull,.true.)
end subroutine FEM_mechanical_init
@ -346,12 +346,11 @@ type(tSolutionState) function FEM_mechanical_solution( &
CHKERRQ(err_PETSc)
call SNESGetConvergedReason(mechanical_snes,reason,err_PETSc) ! solution converged?
CHKERRQ(err_PETSc)
terminallyIll = .false.
if (reason < 1) then ! 0: still iterating (will not occur), negative -> convergence error
FEM_mechanical_solution%converged = .false.
FEM_mechanical_solution%iterationsNeeded = num%itmax
else ! >= 1 proper convergence (or terminally ill)
else ! >= 1 proper convergence (or broken)
FEM_mechanical_solution%converged = .true.
call SNESGetIterationNumber(mechanical_snes,FEM_mechanical_solution%iterationsNeeded,err_PETSc)
CHKERRQ(err_PETSc)
@ -459,8 +458,8 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
call utilities_constitutiveResponse(params%Delta_t,P_av,ForwardData)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
call utilities_constitutiveResponse(broken,params%Delta_t,P_av,ForwardData)
call MPI_Allreduce(MPI_IN_PLACE,broken,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
call parallelization_chkerr(err_MPI)
ForwardData = .false.
@ -748,7 +747,7 @@ subroutine FEM_mechanical_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reaso
divTol = max(maxval(abs(P_av(1:dimPlex,1:dimPlex)))*num%eps_struct_rtol,num%eps_struct_atol)
call SNESConvergedDefault(snes_local,PETScIter,xnorm,snorm,fnorm/divTol,reason,dummy,err_PETSc)
CHKERRQ(err_PETSc)
if (terminallyIll) reason = SNES_DIVERGED_FUNCTION_DOMAIN
if (broken) reason = SNES_DIVERGED_FUNCTION_DOMAIN
print'(/,1x,a,a,i0,a,f0.3)', trim(incInfo), &
' @ Iteration ',PETScIter,' mechanical residual norm = ',fnorm/divTol
print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &