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

removed terminallyIll as global variable

This commit is contained in:
Martin Diehl 2024-01-08 15:22:21 +00:00 committed by Philip Eisenlohr
parent 303061fb96
commit 1f9a372aa7
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 elseif (lovl == 6) then ! stress requested by marc
computationMode = materialpoint_CALCRESULTS computationMode = materialpoint_CALCRESULTS
if (cptim > theTime .or. inc /= theInc) then ! reached "convergence" if (cptim > theTime .or. inc /= theInc) then ! reached "convergence"
terminallyIll = .false.
cycleCounter = -1 ! first calc step increments this to cycle = 0 cycleCounter = -1 ! first calc step increments this to cycle = 0
if (inc == 0) then ! >> start of analysis << if (inc == 0) then ! >> start of analysis <<
lastIncConverged = .false. 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 << else if ( timinc < theDelta ) then ! >> cutBack <<
lastIncConverged = .false. lastIncConverged = .false.
outdatedByNewInc = .false. outdatedByNewInc = .false.
terminallyIll = .false.
cycleCounter = -1 ! first calc step increments this to cycle = 0 cycleCounter = -1 ! first calc step increments this to cycle = 0
print'(a,i6,1x,i2)', '<< HYPELA2 >> cutback detected..! ',m(1),nn print'(a,i6,1x,i2)', '<< HYPELA2 >> cutback detected..! ',m(1),nn
end if ! convergence treatment end end if ! convergence treatment end

15
src/Marc/materialpoint_Marc.f90
View File

@ -37,6 +37,8 @@ module materialpoint_Marc
integer, public :: & integer, public :: &
cycleCounter = 0 !< needs description cycleCounter = 0 !< needs description
logical, public :: &
broken = .false. !< needs description
integer, parameter, public :: & integer, parameter, public :: &
materialpoint_CALCRESULTS = 2**0, & 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 integer elCP, & ! crystal plasticity element number
i, j, k, l, m, n, ph, homog, mySource,ce i, j, k, l, m, n, ph, homog, mySource,ce
real(pREAL), parameter :: ODD_STRESS = 1e15_pREAL, & !< return value for stress if terminallyIll real(pREAL), parameter :: ODD_STRESS = 1e15_pREAL, & !< return value for stress if broken
ODD_JACOBIAN = 1e50_pREAL !< return value for jacobian if terminallyIll ODD_JACOBIAN = 1e50_pREAL !< return value for jacobian if broken
elCP = discretization_Marc_FEM2DAMASK_elem(elFE) elCP = discretization_Marc_FEM2DAMASK_elem(elFE)
ce = discretization_Marc_FEM2DAMASK_cell(ip,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 if (iand(mode, materialpoint_CALCRESULTS) /= 0) then
validCalculation: if (terminallyIll) then validCalculation: if (broken) then
call random_number(rnd) call random_number(rnd)
if (rnd < 0.5_pREAL) rnd = rnd - 1.0_pREAL if (rnd < 0.5_pREAL) rnd = rnd - 1.0_pREAL
materialpoint_cs(1:6,ip,elCP) = ODD_STRESS * rnd materialpoint_cs(1:6,ip,elCP) = ODD_STRESS * rnd
materialpoint_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_eye(6) materialpoint_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_eye(6)
else validCalculation else validCalculation
call homogenization_mechanical_response(dt,(elCP-1)*discretization_nIPs + ip, & call homogenization_mechanical_response(broken, &
(elCP-1)*discretization_nIPs + ip) dt,(elCP-1)*discretization_nIPs + ip, (elCP-1)*discretization_nIPs + ip)
terminalIllness: if (terminallyIll) then terminalIllness: if (broken) then
call random_number(rnd) call random_number(rnd)
if (rnd < 0.5_pREAL) rnd = rnd - 1.0_pREAL 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) :: & real(pREAL) :: &
err_BC !< deviation from stress BC err_BC !< deviation from stress BC
integer :: & integer :: totalIter = 0 !< total iteration in current increment
totalIter = 0 !< total iteration in current increment logical :: broken
public :: & public :: &
grid_mechanical_FEM_init, & grid_mechanical_FEM_init, &
@ -271,7 +271,7 @@ subroutine grid_mechanical_FEM_init(num_grid)
end if restartRead end if restartRead
call utilities_updateCoords(F) 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 F, & ! target F
0.0_pREAL) ! time increment 0.0_pREAL) ! time increment
call DMDAVecRestoreArrayF90(DM_mech,u_PETSc,u,err_PETSc) 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%converged = reason > 0
solution%iterationsNeeded = totalIter solution%iterationsNeeded = totalIter
solution%termIll = terminallyIll solution%termIll = broken
terminallyIll = .false.
P_aim = merge(P_av,P_aim,params%stress_mask) P_aim = merge(P_av,P_aim,params%stress_mask)
end function grid_mechanical_FEM_solution end function grid_mechanical_FEM_solution
@ -336,7 +335,7 @@ end function grid_mechanical_FEM_solution
!> @brief forwarding routine !> @brief forwarding routine
!> @details find new boundary conditions and best F estimate for end of current timestep !> @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) deformation_BC,stress_BC,rotation_BC)
logical, intent(in) :: & 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) 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)) & if ((totalIter >= num%itmin .and. all([err_div/divTol, err_BC/BCTol] < 1.0_pREAL)) &
.or. terminallyIll) then .or. broken) then
reason = 1 reason = 1
elseif (totalIter >= num%itmax) then elseif (totalIter >= num%itmax) then
reason = -1 reason = -1
@ -571,10 +570,10 @@ subroutine formResidual(da_local,x_local, &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate constitutive response ! evaluate constitutive response
call utilities_constitutiveResponse(P_current,& call utilities_constitutiveResponse(broken,P_current,&
P_av,C_volAvg,devNull, & P_av,C_volAvg,devNull, &
F,params%Delta_t,params%rotation_BC) 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) 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_BC, & !< deviation from stress BC
err_div !< RMS of div of P err_div !< RMS of div of P
integer :: & integer :: totalIter = 0 !< total iteration in current increment
totalIter = 0 !< total iteration in current increment logical :: broken
public :: & public :: &
grid_mechanical_spectral_basic_init, & grid_mechanical_spectral_basic_init, &
@ -228,7 +228,7 @@ subroutine grid_mechanical_spectral_basic_init(num_grid)
end if restartRead end if restartRead
call utilities_updateCoords(reshape(F,shape(F_lastInc))) 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 reshape(F,shape(F_lastInc)), & ! target F
0.0_pREAL) ! time increment 0.0_pREAL) ! time increment
call DMDAVecRestoreArrayF90(DM_mech,F_PETSc,F,err_PETSc) ! deassociate pointer 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%converged = reason > 0
solution%iterationsNeeded = totalIter solution%iterationsNeeded = totalIter
solution%termIll = terminallyIll solution%termIll = broken
terminallyIll = .false.
P_aim = merge(P_av,P_aim,params%stress_mask) P_aim = merge(P_av,P_aim,params%stress_mask)
end function grid_mechanical_spectral_basic_solution 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) 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)) & if ((totalIter >= num%itmin .and. all([err_div/divTol, err_BC/BCTol] < 1.0_pREAL)) &
.or. terminallyIll) then .or. broken) then
reason = 1 reason = 1
elseif (totalIter >= num%itmax) then elseif (totalIter >= num%itmax) then
reason = -1 reason = -1
@ -515,10 +514,10 @@ subroutine formResidual(residual_subdomain, F, &
end if newIteration end if newIteration
associate (P => r) associate (P => r)
call utilities_constitutiveResponse(P, & call utilities_constitutiveResponse(broken,P, &
P_av,C_volAvg,C_minMaxAvg, & P_av,C_volAvg,C_minMaxAvg, &
F,params%Delta_t,params%rotation_BC) 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) call parallelization_chkerr(err_MPI)
err_div = utilities_divergenceRMS(P) err_div = utilities_divergenceRMS(P)
end associate 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_curl, & !< RMS of curl of F
err_div !< RMS of div of P err_div !< RMS of div of P
integer :: & integer :: totalIter = 0 !< total iteration in current increment
totalIter = 0 !< total iteration in current increment logical :: broken
public :: & public :: &
grid_mechanical_spectral_polarization_init, & grid_mechanical_spectral_polarization_init, &
@ -257,7 +257,7 @@ subroutine grid_mechanical_spectral_polarization_init(num_grid)
end if restartRead end if restartRead
call utilities_updateCoords(reshape(F,shape(F_lastInc))) 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 reshape(F,shape(F_lastInc)), & ! target F
0.0_pREAL) ! time increment 0.0_pREAL) ! time increment
call DMDAVecRestoreArrayF90(DM_mech,FandF_tau_PETSc,FandF_tau,err_PETSc) ! deassociate pointer 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%converged = reason > 0
solution%iterationsNeeded = totalIter solution%iterationsNeeded = totalIter
solution%termIll = terminallyIll solution%termIll = broken
terminallyIll = .false.
P_aim = merge(P_av,P_aim,params%stress_mask) P_aim = merge(P_av,P_aim,params%stress_mask)
end function grid_mechanical_spectral_polarization_solution 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) 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)) & 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 reason = 1
elseif (totalIter >= num%itmax) then elseif (totalIter >= num%itmax) then
reason = -1 reason = -1
@ -605,14 +604,14 @@ subroutine formResidual(residual_subdomain, FandF_tau, &
err_curl = utilities_curlRMS(F) err_curl = utilities_curlRMS(F)
#ifdef __GFORTRAN__ #ifdef __GFORTRAN__
call utilities_constitutiveResponse(r_F, & call utilities_constitutiveResponse(broken,r_F, &
#else #else
associate (P => r_F) associate (P => r_F)
call utilities_constitutiveResponse(P, & call utilities_constitutiveResponse(broken, P, &
#endif #endif
P_av,C_volAvg,C_minMaxAvg, & P_av,C_volAvg,C_minMaxAvg, &
F - r_F_tau/num%beta,params%Delta_t,params%rotation_BC) 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__ #ifdef __GFORTRAN__
err_div = utilities_divergenceRMS(r_F) err_div = utilities_divergenceRMS(r_F)
#else #else

7
src/grid/grid_mech_utilities.f90
View File

@ -113,9 +113,10 @@ end function utilities_maskedCompliance
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Calculate constitutive response. !> @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) 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,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) :: P_av !< average PK stress
real(pREAL), intent(out), dimension(3,3,cells(1),cells(2),cells3) :: P !< 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 real(pREAL), intent(in) :: Delta_t !< loading time
type(tRotation), intent(in), optional :: rotation_BC !< rotation of load frame type(tRotation), intent(in), optional :: rotation_BC !< rotation of load frame
integer :: i integer :: i
integer(MPI_INTEGER_KIND) :: err_MPI integer(MPI_INTEGER_KIND) :: err_MPI
real(pREAL), dimension(3,3,3,3) :: dPdF_max, dPdF_min real(pREAL), dimension(3,3,3,3) :: dPdF_max, dPdF_min
real(pREAL) :: dPdF_norm_max, dPdF_norm_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 real(pREAL), dimension(2) :: valueAndRank !< pair of min/max norm of dPdF to synchronize min/max of dPdF
print'(/,1x,a)', '... evaluating constitutive response ......................................' print'(/,1x,a)', '... evaluating constitutive response ......................................'
flush(IO_STDOUT) flush(IO_STDOUT)
homogenization_F = reshape(F,[3,3,product(cells(1:2))*cells3]) ! set materialpoint target F to estimated field 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 = reshape(homogenization_P, [3,3,cells(1),cells(2),cells3])
P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt 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 integer :: totalIter = 0 !< total iteration in current increment
real(pREAL), dimension(3,3) :: K_ref real(pREAL), dimension(3,3) :: K_ref
real(pREAL) :: mu_ref, Delta_t_ real(pREAL) :: mu_ref, Delta_t_
logical :: broken
public :: & public :: &
grid_thermal_spectral_init, & grid_thermal_spectral_init, &
@ -206,7 +207,7 @@ function grid_thermal_spectral_solution(Delta_t) result(solution)
call SNESGetConvergedReason(SNES_thermal,reason,err_PETSc) call SNESGetConvergedReason(SNES_thermal,reason,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
solution%converged = reason > 0 solution%converged = reason > 0 .and. .not. broken
solution%iterationsNeeded = merge(totalIter,num%itmax,solution%converged) solution%iterationsNeeded = merge(totalIter,num%itmax,solution%converged)
call SNESGetDM(SNES_thermal,DM_thermal,err_PETSc) 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 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) associate(T => x_scal)
vectorField = utilities_ScalarGradient(T) vectorField = utilities_ScalarGradient(T)

42
src/homogenization.f90
View File

@ -46,9 +46,6 @@ module homogenization
thermal_active, & thermal_active, &
damage_active damage_active
logical, public :: &
terminallyIll = .false. !< at least one material point is terminally ill
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! General variables for the homogenization at a material point ! General variables for the homogenization at a material point
real(pREAL), dimension(:,:,:), allocatable, public :: & real(pREAL), dimension(:,:,:), allocatable, public :: &
@ -213,8 +210,9 @@ end subroutine homogenization_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief !> @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 real(pREAL), intent(in) :: Delta_t !< time increment
integer, intent(in) :: & integer, intent(in) :: &
cell_start, cell_end cell_start, cell_end
@ -226,8 +224,8 @@ subroutine homogenization_mechanical_response(Delta_t,cell_start,cell_end)
doneAndHappy doneAndHappy
!$OMP PARALLEL broken = .false.
!$OMP DO PRIVATE(en,ho,co,converged,doneAndHappy) !$OMP PARALLEL DO PRIVATE(en,ho,co,converged,doneAndHappy)
do ce = cell_start, cell_end do ce = cell_start, cell_end
en = material_entry_homogenization(ce) en = material_entry_homogenization(ce)
@ -241,7 +239,7 @@ subroutine homogenization_mechanical_response(Delta_t,cell_start,cell_end)
doneAndHappy = [.false.,.true.] 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) 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))]) 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.] doneAndHappy = [.true.,.false.]
end if end if
end do convergenceLooping 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))]) converged = converged .and. all([(phase_damage_constitutive(Delta_t,co,ce),co=1,homogenization_Nconstituents(ho))])
if (.not. converged) then if (.not. converged) then
if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill' if (.not. broken) print*, ' Cell ', ce, ' failed (damage)'
terminallyIll = .true. broken = .true.
end if end if
end do 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 do ce = cell_start, cell_end
ho = material_ID_homogenization(ce) ho = material_ID_homogenization(ce)
do co = 1, homogenization_Nconstituents(ho) do co = 1, homogenization_Nconstituents(ho)
@ -270,8 +273,7 @@ subroutine homogenization_mechanical_response(Delta_t,cell_start,cell_end)
end do end do
call mechanical_homogenize(Delta_t,ce) call mechanical_homogenize(Delta_t,ce)
end do end do
!$OMP END DO !$OMP END PARALLEL DO
!$OMP END PARALLEL
end subroutine homogenization_mechanical_response end subroutine homogenization_mechanical_response
@ -279,8 +281,10 @@ end subroutine homogenization_mechanical_response
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief !> @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 real(pREAL), intent(in) :: Delta_t !< time increment
integer, intent(in) :: & integer, intent(in) :: &
cell_start, cell_end cell_start, cell_end
@ -289,18 +293,20 @@ subroutine homogenization_thermal_response(Delta_t,cell_start,cell_end)
co, ce, ho co, ce, ho
broken = .false.
!$OMP PARALLEL DO PRIVATE(ho) !$OMP PARALLEL DO PRIVATE(ho)
do ce = cell_start, cell_end do ce = cell_start, cell_end
if (terminallyIll) continue if (broken) continue
ho = material_ID_homogenization(ce) ho = material_ID_homogenization(ce)
do co = 1, homogenization_Nconstituents(ho) 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. phase_thermal_constitutive(Delta_t,material_ID_phase(co,ce),material_entry_phase(co,ce))) then
if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill' if (.not. broken) print*, ' Cell ', ce, ' failed (thermal)'
terminallyIll = .true. broken = .true.
end if end if
end do end do
end do end do
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
broken = broken
end subroutine homogenization_thermal_response 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 !> @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 real(pREAL), intent(in) :: Delta_t !< loading time
logical, intent(in) :: forwardData !< age results logical, intent(in) :: forwardData !< age results
real(pREAL),intent(out), dimension(3,3) :: P_av !< average PK stress 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 ......................................' 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. cutBack = .false.
P_av = sum(homogenization_P,dim=3) * wgt 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 character(len=pSTRLEN) :: incInfo
real(pREAL), dimension(3,3) :: & real(pREAL), dimension(3,3) :: &
P_av = 0.0_pREAL P_av = 0.0_pREAL
logical :: ForwardData logical :: ForwardData, broken
real(pREAL), parameter :: eps = 1.0e-18_pREAL real(pREAL), parameter :: eps = 1.0e-18_pREAL
external :: & ! ToDo: write interfaces 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) call DMPlexVecSetClosure(mechanical_mesh,section,solution_local,cell,px_scal,5,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
end do end do
call utilities_constitutiveResponse(0.0_pREAL,devNull,.true.) call utilities_constitutiveResponse(broken,0.0_pREAL,devNull,.true.)
end subroutine FEM_mechanical_init end subroutine FEM_mechanical_init
@ -346,12 +346,11 @@ type(tSolutionState) function FEM_mechanical_solution( &
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call SNESGetConvergedReason(mechanical_snes,reason,err_PETSc) ! solution converged? call SNESGetConvergedReason(mechanical_snes,reason,err_PETSc) ! solution converged?
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
terminallyIll = .false.
if (reason < 1) then ! 0: still iterating (will not occur), negative -> convergence error if (reason < 1) then ! 0: still iterating (will not occur), negative -> convergence error
FEM_mechanical_solution%converged = .false. FEM_mechanical_solution%converged = .false.
FEM_mechanical_solution%iterationsNeeded = num%itmax FEM_mechanical_solution%iterationsNeeded = num%itmax
else ! >= 1 proper convergence (or terminally ill) else ! >= 1 proper convergence (or broken)
FEM_mechanical_solution%converged = .true. FEM_mechanical_solution%converged = .true.
call SNESGetIterationNumber(mechanical_snes,FEM_mechanical_solution%iterationsNeeded,err_PETSc) call SNESGetIterationNumber(mechanical_snes,FEM_mechanical_solution%iterationsNeeded,err_PETSc)
CHKERRQ(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 ! evaluate constitutive response
call utilities_constitutiveResponse(params%Delta_t,P_av,ForwardData) call utilities_constitutiveResponse(broken,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 MPI_Allreduce(MPI_IN_PLACE,broken,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
call parallelization_chkerr(err_MPI) call parallelization_chkerr(err_MPI)
ForwardData = .false. 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) 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) call SNESConvergedDefault(snes_local,PETScIter,xnorm,snorm,fnorm/divTol,reason,dummy,err_PETSc)
CHKERRQ(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), & print'(/,1x,a,a,i0,a,f0.3)', trim(incInfo), &
' @ Iteration ',PETScIter,' mechanical residual norm = ',fnorm/divTol ' @ Iteration ',PETScIter,' mechanical residual norm = ',fnorm/divTol
print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', & print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &