From 1f9a372aa751a22d0c4a023489658404c1491d00 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 8 Jan 2024 15:22:21 +0000 Subject: [PATCH] removed terminallyIll as global variable --- src/Marc/DAMASK_Marc.f90 | 2 - src/Marc/materialpoint_Marc.f90 | 15 +++---- src/grid/grid_mech_FEM.f90 | 19 +++++---- src/grid/grid_mech_spectral_basic.f90 | 15 ++++--- src/grid/grid_mech_spectral_polarization.f90 | 17 ++++---- src/grid/grid_mech_utilities.f90 | 7 ++-- src/grid/grid_thermal_spectral.f90 | 5 ++- src/homogenization.f90 | 42 +++++++++++--------- src/mesh/FEM_utilities.f90 | 5 ++- src/mesh/mesh_mech_FEM.f90 | 13 +++--- 10 files changed, 72 insertions(+), 68 deletions(-) diff --git a/src/Marc/DAMASK_Marc.f90 b/src/Marc/DAMASK_Marc.f90 index ec6f34923..539bc49d2 100644 --- a/src/Marc/DAMASK_Marc.f90 +++ b/src/Marc/DAMASK_Marc.f90 @@ -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 diff --git a/src/Marc/materialpoint_Marc.f90 b/src/Marc/materialpoint_Marc.f90 index 4d8f333ec..18ba90a61 100644 --- a/src/Marc/materialpoint_Marc.f90 +++ b/src/Marc/materialpoint_Marc.f90 @@ -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 diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 2b856ff1c..140b22817 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -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,8 +335,8 @@ 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,& - deformation_BC,stress_BC,rotation_BC) +subroutine grid_mechanical_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining, & + deformation_BC,stress_BC,rotation_BC) logical, intent(in) :: & cutBack, & @@ -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) !-------------------------------------------------------------------------------------------------- diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 22113a6fd..cc92d4e13 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -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 diff --git a/src/grid/grid_mech_spectral_polarization.f90 b/src/grid/grid_mech_spectral_polarization.f90 index a4da7452d..4328e591f 100644 --- a/src/grid/grid_mech_spectral_polarization.f90 +++ b/src/grid/grid_mech_spectral_polarization.f90 @@ -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 diff --git a/src/grid/grid_mech_utilities.f90 b/src/grid/grid_mech_utilities.f90 index bae5c309b..a767c8623 100644 --- a/src/grid/grid_mech_utilities.f90 +++ b/src/grid/grid_mech_utilities.f90 @@ -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 diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 1c46050ff..35f2b47b4 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -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) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 2da1d73b2..ab2b8a8b1 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -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 diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 220852b89..9f1584795 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -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 diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 13ae970ae..0baebe6c3 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -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))', &