diff --git a/src/Marc/DAMASK_Marc.f90 b/src/Marc/DAMASK_Marc.f90 index 539bc49d2..6f34f5354 100644 --- a/src/Marc/DAMASK_Marc.f90 +++ b/src/Marc/DAMASK_Marc.f90 @@ -16,8 +16,8 @@ #endif #include "../prec.f90" -#include "../parallelization.f90" #include "../constants.f90" +#include "../parallelization.f90" #include "../misc.f90" #include "../IO.f90" #include "../YAML_types.f90" diff --git a/src/Marc/materialpoint_Marc.f90 b/src/Marc/materialpoint_Marc.f90 index 18ba90a61..f0040cba2 100644 --- a/src/Marc/materialpoint_Marc.f90 +++ b/src/Marc/materialpoint_Marc.f90 @@ -20,6 +20,7 @@ module materialpoint_Marc use material use phase use homogenization + use constants use discretization use discretization_Marc @@ -27,18 +28,18 @@ module materialpoint_Marc implicit none(type,external) private - real(pREAL), dimension (:,:,:), allocatable, private :: & + real(pREAL), dimension (:,:,:), allocatable :: & materialpoint_cs !< Cauchy stress - real(pREAL), dimension (:,:,:,:), allocatable, private :: & + real(pREAL), dimension (:,:,:,:), allocatable :: & materialpoint_dcsdE, & !< Cauchy stress tangent materialpoint_F !< deformation gradient - real(pREAL), dimension (:,:,:,:), allocatable, private :: & + real(pREAL), dimension (:,:,:,:), allocatable :: & materialpoint_dcsdE_knownGood !< known good tangent integer, public :: & cycleCounter = 0 !< needs description - logical, public :: & - broken = .false. !< needs description + integer(kind(STATUS_OK)) :: & + status integer, parameter, public :: & materialpoint_CALCRESULTS = 2**0, & @@ -150,17 +151,17 @@ subroutine materialpoint_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, if (iand(mode, materialpoint_CALCRESULTS) /= 0) then - validCalculation: if (broken) then + validCalculation: if (status /= STATUS_OK) 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(broken, & + call homogenization_mechanical_response(status, & dt,(elCP-1)*discretization_nIPs + ip, (elCP-1)*discretization_nIPs + ip) - terminalIllness: if (broken) then + terminalIllness: if (status /= STATUS_OK) then call random_number(rnd) if (rnd < 0.5_pREAL) rnd = rnd - 1.0_pREAL diff --git a/src/constants.f90 b/src/constants.f90 index 1402154c7..38cdcb195 100644 --- a/src/constants.f90 +++ b/src/constants.f90 @@ -1,6 +1,6 @@ !-------------------------------------------------------------------------------------------------- !> @author Martin Diehl, KU Leuven -!> @brief physical constants +!> @brief Constants. !-------------------------------------------------------------------------------------------------- module constants use prec @@ -20,4 +20,18 @@ module constants character(len=*), parameter :: LOWER = 'abcdefghijklmnopqrstuvwxyz' character(len=len(LOWER)), parameter :: UPPER = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' + enum, bind(c); enumerator :: & + STATUS_OK, & + STATUS_ITERATING, & + STATUS_FAIL_PHASE_MECHANICAL, & + STATUS_FAIL_PHASE_MECHANICAL_STATE, & + STATUS_FAIL_PHASE_MECHANICAL_DELTASTATE, & + STATUS_FAIL_PHASE_MECHANICAL_STRESS, & + STATUS_FAIL_PHASE_DAMAGE, & + STATUS_FAIL_PHASE_DAMAGE_STATE, & + STATUS_FAIL_PHASE_DAMAGE_DELTASTATE, & + STATUS_FAIL_PHASE_THERMAL, & + STATUS_FAIL_PHASE_THERMAL_DOTSTATE + end enum + end module constants diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 140b22817..00653a508 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -28,7 +28,7 @@ module grid_mechanical_FEM use homogenization use discretization use discretization_grid - + use constants #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) implicit none(type,external) @@ -84,7 +84,7 @@ module grid_mechanical_FEM err_BC !< deviation from stress BC integer :: totalIter = 0 !< total iteration in current increment - logical :: broken + integer(kind(STATUS_OK)) :: status 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(broken,P_current,P_av,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2 + call utilities_constitutiveResponse(status,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,7 +325,7 @@ function grid_mechanical_FEM_solution(incInfoIn) result(solution) solution%converged = reason > 0 solution%iterationsNeeded = totalIter - solution%termIll = broken + solution%termIll = status /= STATUS_OK P_aim = merge(P_av,P_aim,params%stress_mask) end function grid_mechanical_FEM_solution @@ -492,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. broken) then + .or. status /= STATUS_OK) then reason = 1 elseif (totalIter >= num%itmax) then reason = -1 @@ -525,7 +525,7 @@ subroutine formResidual(da_local,x_local, & real(pREAL), pointer,dimension(:,:,:,:) :: x_scal, r real(pREAL), dimension(8,3) :: x_elem, f_elem - PetscInt :: i, ii, j, jj, k, kk, ctr, ele + PetscInt :: i, ii, j, jj, k, kk, ctr, ce PetscInt :: & PETScIter, & nfuncs @@ -570,10 +570,10 @@ subroutine formResidual(da_local,x_local, & !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - call utilities_constitutiveResponse(broken,P_current,& + call utilities_constitutiveResponse(status,P_current,& P_av,C_volAvg,devNull, & F,params%Delta_t,params%rotation_BC) - call MPI_Allreduce(MPI_IN_PLACE,broken,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) + call MPI_Allreduce(MPI_IN_PLACE,status,1_MPI_INTEGER_KIND,MPI_INTEGER,MPI_MAX,MPI_COMM_WORLD,err_MPI) call parallelization_chkerr(err_MPI) !-------------------------------------------------------------------------------------------------- @@ -587,7 +587,7 @@ subroutine formResidual(da_local,x_local, & CHKERRQ(err_PETSc) call DMDAVecGetArrayReadF90(da_local,x_local,x_scal,err_PETSc) CHKERRQ(err_PETSc) - ele = 0 + ce = 0 r = 0.0_pREAL do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, cells(1) ctr = 0 @@ -595,11 +595,11 @@ subroutine formResidual(da_local,x_local, & ctr = ctr + 1 x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk) end do; end do; end do - ele = ele + 1 + ce = ce + 1 f_elem = matmul(transpose(BMat),transpose(P_current(1:3,1:3,i,j,k-cells3Offset)))*detJ + & - matmul(HGMat,x_elem)*(homogenization_dPdF(1,1,1,1,ele) + & - homogenization_dPdF(2,2,2,2,ele) + & - homogenization_dPdF(3,3,3,3,ele))/3.0_pREAL + matmul(HGMat,x_elem)*(homogenization_dPdF(1,1,1,1,ce) + & + homogenization_dPdF(2,2,2,2,ce) + & + homogenization_dPdF(3,3,3,3,ce))/3.0_pREAL ctr = 0 do kk = -1, 0; do jj = -1, 0; do ii = -1, 0 ctr = ctr + 1 diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index cc92d4e13..d71181479 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -26,6 +26,7 @@ module grid_mechanical_spectral_basic use grid_mech_utilities use homogenization use discretization_grid + use constants #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) implicit none(type,external) @@ -84,7 +85,7 @@ module grid_mechanical_spectral_basic err_div !< RMS of div of P integer :: totalIter = 0 !< total iteration in current increment - logical :: broken + integer(kind(STATUS_OK)) :: status public :: & grid_mechanical_spectral_basic_init, & @@ -228,7 +229,7 @@ subroutine grid_mechanical_spectral_basic_init(num_grid) end if restartRead call utilities_updateCoords(reshape(F,shape(F_lastInc))) - call utilities_constitutiveResponse(broken,P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 + call utilities_constitutiveResponse(status,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,7 +288,7 @@ function grid_mechanical_spectral_basic_solution(incInfoIn) result(solution) solution%converged = reason > 0 solution%iterationsNeeded = totalIter - solution%termIll = broken + solution%termIll = status /= STATUS_OK P_aim = merge(P_av,P_aim,params%stress_mask) end function grid_mechanical_spectral_basic_solution @@ -452,7 +453,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. broken) then + .or. status /= STATUS_OK) then reason = 1 elseif (totalIter >= num%itmax) then reason = -1 @@ -514,10 +515,10 @@ subroutine formResidual(residual_subdomain, F, & end if newIteration associate (P => r) - call utilities_constitutiveResponse(broken,P, & + call utilities_constitutiveResponse(status,P, & P_av,C_volAvg,C_minMaxAvg, & F,params%Delta_t,params%rotation_BC) - call MPI_Allreduce(MPI_IN_PLACE,broken,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) + call MPI_Allreduce(MPI_IN_PLACE,status,1_MPI_INTEGER_KIND,MPI_INTEGER,MPI_MAX,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 4328e591f..7304d3640 100644 --- a/src/grid/grid_mech_spectral_polarization.f90 +++ b/src/grid/grid_mech_spectral_polarization.f90 @@ -27,6 +27,7 @@ module grid_mechanical_spectral_polarization use config use homogenization use discretization_grid + use constants #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) implicit none(type,external) @@ -96,7 +97,7 @@ module grid_mechanical_spectral_polarization err_div !< RMS of div of P integer :: totalIter = 0 !< total iteration in current increment - logical :: broken + integer(kind(STATUS_OK)) :: status public :: & grid_mechanical_spectral_polarization_init, & @@ -257,7 +258,7 @@ subroutine grid_mechanical_spectral_polarization_init(num_grid) end if restartRead call utilities_updateCoords(reshape(F,shape(F_lastInc))) - call utilities_constitutiveResponse(broken,P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 + call utilities_constitutiveResponse(status,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,7 +323,7 @@ function grid_mechanical_spectral_polarization_solution(incInfoIn) result(soluti solution%converged = reason > 0 solution%iterationsNeeded = totalIter - solution%termIll = broken + solution%termIll = status /= STATUS_OK P_aim = merge(P_av,P_aim,params%stress_mask) end function grid_mechanical_spectral_polarization_solution @@ -516,7 +517,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. broken) then + .or. status /= STATUS_OK) then reason = 1 elseif (totalIter >= num%itmax) then reason = -1 @@ -562,7 +563,7 @@ subroutine formResidual(residual_subdomain, FandF_tau, & nfuncs integer(MPI_INTEGER_KIND) :: err_MPI integer :: & - i, j, k, e + i, j, k, ce F => FandF_tau(1:3,1:3,1,1:cells(1),1:cells(2),1:cells3) @@ -604,24 +605,24 @@ subroutine formResidual(residual_subdomain, FandF_tau, & err_curl = utilities_curlRMS(F) #ifdef __GFORTRAN__ - call utilities_constitutiveResponse(broken,r_F, & + call utilities_constitutiveResponse(status,r_F, & #else associate (P => r_F) - call utilities_constitutiveResponse(broken, P, & + call utilities_constitutiveResponse(status, 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,broken,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) + call MPI_Allreduce(MPI_IN_PLACE,status,1_MPI_INTEGER_KIND,MPI_INTEGER,MPI_MAX,MPI_COMM_WORLD,err_MPI) #ifdef __GFORTRAN__ err_div = utilities_divergenceRMS(r_F) #else err_div = utilities_divergenceRMS(P) #endif - e = 0 + ce = 0 do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1) - e = e + 1 + ce = ce + 1 r_F(1:3,1:3,i,j,k) = & - math_mul3333xx33(math_invSym3333(homogenization_dPdF(1:3,1:3,1:3,1:3,e) + C_scale), & + math_mul3333xx33(math_invSym3333(homogenization_dPdF(1:3,1:3,1:3,1:3,ce) + C_scale), & #ifdef __GFORTRAN__ r_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), & #else diff --git a/src/grid/grid_mech_utilities.f90 b/src/grid/grid_mech_utilities.f90 index a767c8623..be426497e 100644 --- a/src/grid/grid_mech_utilities.f90 +++ b/src/grid/grid_mech_utilities.f90 @@ -18,7 +18,7 @@ module grid_mech_utilities use discretization use spectral_utilities use homogenization - + use constants #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) implicit none(type,external) @@ -113,10 +113,10 @@ end function utilities_maskedCompliance !-------------------------------------------------------------------------------------------------- !> @brief Calculate constitutive response. !-------------------------------------------------------------------------------------------------- -subroutine utilities_constitutiveResponse(broken, P,P_av,C_volAvg,C_minmaxAvg,& +subroutine utilities_constitutiveResponse(status, P,P_av,C_volAvg,C_minmaxAvg,& F,Delta_t,rotation_BC) - logical, intent(out) :: broken + integer(kind(STATUS_OK)), intent(out) :: status 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 @@ -124,7 +124,7 @@ subroutine utilities_constitutiveResponse(broken, 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 :: ce 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 @@ -136,7 +136,7 @@ subroutine utilities_constitutiveResponse(broken, P,P_av,C_volAvg,C_minmaxAvg,& homogenization_F = reshape(F,[3,3,product(cells(1:2))*cells3]) ! set materialpoint target F to estimated field - call homogenization_mechanical_response(broken,Delta_t,1,product(cells(1:2))*cells3) ! calculate P field + call homogenization_mechanical_response(status,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 @@ -156,14 +156,14 @@ subroutine utilities_constitutiveResponse(broken, P,P_av,C_volAvg,C_minmaxAvg,& dPdF_norm_max = 0.0_pREAL dPdF_min = huge(1.0_pREAL) dPdF_norm_min = huge(1.0_pREAL) - do i = 1, product(cells(1:2))*cells3 - if (dPdF_norm_max < sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)) then - dPdF_max = homogenization_dPdF(1:3,1:3,1:3,1:3,i) - dPdF_norm_max = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2) + do ce = 1, product(cells(1:2))*cells3 + if (dPdF_norm_max < sum(homogenization_dPdF(1:3,1:3,1:3,1:3,ce)**2)) then + dPdF_max = homogenization_dPdF(1:3,1:3,1:3,1:3,ce) + dPdF_norm_max = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,ce)**2) end if - if (dPdF_norm_min > sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)) then - dPdF_min = homogenization_dPdF(1:3,1:3,1:3,1:3,i) - dPdF_norm_min = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2) + if (dPdF_norm_min > sum(homogenization_dPdF(1:3,1:3,1:3,1:3,ce)**2)) then + dPdF_min = homogenization_dPdF(1:3,1:3,1:3,1:3,ce) + dPdF_norm_min = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,ce)**2) end if end do diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 35f2b47b4..cc4b89afb 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -25,6 +25,7 @@ module grid_thermal_spectral use homogenization use YAML_types use config + use constants #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) implicit none(type,external) @@ -56,7 +57,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 + integer(kind(STATUS_OK)) :: status public :: & grid_thermal_spectral_init, & @@ -207,7 +208,7 @@ function grid_thermal_spectral_solution(Delta_t) result(solution) call SNESGetConvergedReason(SNES_thermal,reason,err_PETSc) CHKERRQ(err_PETSc) - solution%converged = reason > 0 .and. .not. broken + solution%converged = reason > 0 .and. status == STATUS_OK solution%iterationsNeeded = merge(totalIter,num%itmax,solution%converged) call SNESGetDM(SNES_thermal,DM_thermal,err_PETSc) @@ -324,7 +325,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(broken,Delta_t_,1,product(cells(1:2))*cells3) + call homogenization_thermal_response(status,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 ab2b8a8b1..13e6910ab 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -50,9 +50,9 @@ module homogenization ! General variables for the homogenization at a material point real(pREAL), dimension(:,:,:), allocatable, public :: & homogenization_F !< def grad of IP to be reached at end of FE increment - real(pREAL), dimension(:,:,:), allocatable, public :: & !, protected :: & Issue with ifort + real(pREAL), dimension(:,:,:), allocatable, public, protected :: & homogenization_P !< first P--K stress of IP - real(pREAL), dimension(:,:,:,:,:), allocatable, public :: & !, protected :: & + real(pREAL), dimension(:,:,:,:,:), allocatable, public, protected :: & homogenization_dPdF !< tangent of first P--K stress at IP !-------------------------------------------------------------------------------------------------- @@ -210,9 +210,9 @@ end subroutine homogenization_init !-------------------------------------------------------------------------------------------------- !> @brief !-------------------------------------------------------------------------------------------------- -subroutine homogenization_mechanical_response(broken,Delta_t,cell_start,cell_end) +subroutine homogenization_mechanical_response(status,Delta_t,cell_start,cell_end) - logical, intent(out) :: broken + integer(kind(STATUS_OK)), intent(out) :: status real(pREAL), intent(in) :: Delta_t !< time increment integer, intent(in) :: & cell_start, cell_end @@ -224,7 +224,7 @@ subroutine homogenization_mechanical_response(broken,Delta_t,cell_start,cell_end doneAndHappy - broken = .false. + status = STATUS_OK !$OMP PARALLEL DO PRIVATE(en,ho,co,converged,doneAndHappy) do ce = cell_start, cell_end @@ -239,10 +239,10 @@ subroutine homogenization_mechanical_response(broken,Delta_t,cell_start,cell_end doneAndHappy = [.false.,.true.] - convergenceLooping: do while (.not. (broken .or. doneAndHappy(1))) + convergenceLooping: do while (status == STATUS_OK .and. .not. 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))]) + converged = all([(phase_mechanical_constitutive(Delta_t,co,ce) == STATUS_OK,co=1,homogenization_Nconstituents(ho))]) if (converged) then doneAndHappy = mechanical_updateState(Delta_t,homogenization_F(1:3,1:3,ce),ce) converged = all(doneAndHappy) @@ -251,19 +251,19 @@ subroutine homogenization_mechanical_response(broken,Delta_t,cell_start,cell_end end if end do convergenceLooping if (.not. converged) then - if (.not. broken) print*, ' Cell ', ce, ' failed (mechanics)' - broken = .true. + if (status == STATUS_OK) print*, ' Cell ', ce, ' failed (mechanics)' + status = STATUS_FAIL_PHASE_MECHANICAL 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)==STATUS_OK,co=1,homogenization_Nconstituents(ho))]) if (.not. converged) then - if (.not. broken) print*, ' Cell ', ce, ' failed (damage)' - broken = .true. + if (status == STATUS_OK) print*, ' Cell ', ce, ' failed (damage)' + status = STATUS_FAIL_PHASE_DAMAGE end if end do !$OMP END PARALLEL DO - if (broken) return + if (status /= STATUS_OK) return !$OMP PARALLEL DO PRIVATE(ho) do ce = cell_start, cell_end @@ -281,10 +281,10 @@ end subroutine homogenization_mechanical_response !-------------------------------------------------------------------------------------------------- !> @brief !-------------------------------------------------------------------------------------------------- -subroutine homogenization_thermal_response(broken, & +subroutine homogenization_thermal_response(status, & Delta_t,cell_start,cell_end) - logical, intent(out) :: broken + integer(kind(STATUS_OK)), intent(out) :: status real(pREAL), intent(in) :: Delta_t !< time increment integer, intent(in) :: & cell_start, cell_end @@ -293,20 +293,19 @@ subroutine homogenization_thermal_response(broken, & co, ce, ho - broken = .false. + status = STATUS_OK !$OMP PARALLEL DO PRIVATE(ho) do ce = cell_start, cell_end - if (broken) continue + if (status /= STATUS_OK) 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. broken) print*, ' Cell ', ce, ' failed (thermal)' - broken = .true. + if (phase_thermal_constitutive(Delta_t,material_ID_phase(co,ce),material_entry_phase(co,ce)) /= STATUS_OK) then + if (status == STATUS_OK) print*, ' Cell ', ce, ' failed (thermal)' + status = STATUS_FAIL_PHASE_THERMAL 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 9f1584795..ea4e4b87f 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -22,6 +22,7 @@ module FEM_utilities use discretization_mesh use homogenization use FEM_quadrature + use constants #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) implicit none(type,external) @@ -120,9 +121,9 @@ end subroutine FEM_utilities_init !-------------------------------------------------------------------------------------------------- !> @brief calculates constitutive response !-------------------------------------------------------------------------------------------------- -subroutine utilities_constitutiveResponse(broken, Delta_t,P_av,forwardData) +subroutine utilities_constitutiveResponse(status, Delta_t,P_av,forwardData) - logical, intent(out) :: broken + integer(kind(STATUS_OK)), intent(out) :: status 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 @@ -132,7 +133,7 @@ subroutine utilities_constitutiveResponse(broken, Delta_t,P_av,forwardData) print'(/,1x,a)', '... evaluating constitutive response ......................................' - call homogenization_mechanical_response(broken,Delta_t,1,mesh_maxNips*mesh_NcpElems) ! calculate P field + call homogenization_mechanical_response(status,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 0baebe6c3..17475ba48 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -25,6 +25,7 @@ module mesh_mechanical_FEM use FEM_quadrature use homogenization use math + use constants #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) implicit none(type,external) @@ -68,7 +69,8 @@ module mesh_mechanical_FEM character(len=pSTRLEN) :: incInfo real(pREAL), dimension(3,3) :: & P_av = 0.0_pREAL - logical :: ForwardData, broken + logical :: ForwardData + integer(kind(STATUS_OK)) :: status real(pREAL), parameter :: eps = 1.0e-18_pREAL external :: & ! ToDo: write interfaces @@ -311,7 +313,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(broken,0.0_pREAL,devNull,.true.) + call utilities_constitutiveResponse(status,0.0_pREAL,devNull,.true.) end subroutine FEM_mechanical_init @@ -458,8 +460,8 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - 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 utilities_constitutiveResponse(status,params%Delta_t,P_av,ForwardData) + call MPI_Allreduce(MPI_IN_PLACE,status,1_MPI_INTEGER_KIND,MPI_INTEGER,MPI_MAX,MPI_COMM_WORLD,err_MPI) call parallelization_chkerr(err_MPI) ForwardData = .false. @@ -529,7 +531,7 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P real(pREAL),dimension(cellDOF,cellDOF) :: K_eA, K_eB PetscInt :: cellStart, cellEnd, cell, component, face, & - qPt, basis, comp, cidx,bcSize, m, i + qPt, basis, comp, cidx,bcSize, ce, i IS :: bcPoints @@ -581,7 +583,7 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P FAvg = 0.0_pREAL BMatAvg = 0.0_pREAL do qPt = 0_pPETSCINT, nQuadrature-1_pPETSCINT - m = cell*nQuadrature + qPt + 1_pPETSCINT + ce = cell*nQuadrature + qPt + 1_pPETSCINT BMat = 0.0_pREAL do basis = 0_pPETSCINT, nBasis-1_pPETSCINT do comp = 0_pPETSCINT, dimPlex-1_pPETSCINT @@ -591,7 +593,7 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P matmul(reshape(pInvcellJ,[dimPlex,dimPlex]),basisFieldDer(i*dimPlex+1_pPETSCINT:(i+1_pPETSCINT)*dimPlex)) end do end do - MatA = matmul(reshape(reshape(homogenization_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,m), & + MatA = matmul(reshape(reshape(homogenization_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,ce), & shape=[dimPlex,dimPlex,dimPlex,dimPlex], order=[2,1,4,3]), & shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1_pPETSCINT) if (num%BBarStabilization) then @@ -599,11 +601,11 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P FInv = math_inv33(F) K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0_pREAL/real(dimPlex,pREAL)) K_eB = K_eB - & - matmul(transpose(matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,m),shape=[dimPlex**2,1_pPETSCINT]), & + matmul(transpose(matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,ce),shape=[dimPlex**2,1_pPETSCINT]), & matmul(reshape(FInv(1:dimPlex,1:dimPlex), & shape=[1_pPETSCINT,dimPlex**2],order=[2,1]),BMat))),MatA) MatB = MatB & - + matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,m),shape=[1_pPETSCINT,dimPlex**2]),MatA) + + matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,ce),shape=[1_pPETSCINT,dimPlex**2]),MatA) FAvg = FAvg + F BMatAvg = BMatAvg + BMat else @@ -747,7 +749,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 (broken) reason = SNES_DIVERGED_FUNCTION_DOMAIN + if (status /= STATUS_OK) 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))', & diff --git a/src/parallelization.f90 b/src/parallelization.f90 index 20ab75336..775a68f6c 100644 --- a/src/parallelization.f90 +++ b/src/parallelization.f90 @@ -7,6 +7,8 @@ module parallelization OUTPUT_UNIT, & ERROR_UNIT + use constants + #ifdef PETSC #include use PETScSys @@ -63,6 +65,7 @@ subroutine parallelization_init() !$ integer(pI32) :: OMP_NUM_THREADS !$ character(len=6) NumThreadsString PetscErrorCode :: err_PETSc + integer(kind(STATUS_OK)) :: status #ifdef _OPENMP @@ -116,28 +119,30 @@ subroutine parallelization_init() #endif call MPI_Comm_size(MPI_COMM_WORLD,worldsize,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) & - error stop 'Could not determine worldsize' + call parallelization_chkerr(err_MPI) if (worldrank == 0) print'(/,1x,a,i0)', 'MPI processes: ',worldsize call MPI_Type_size(MPI_INTEGER,typeSize,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) & - error stop 'Could not determine size of MPI_INTEGER' + call parallelization_chkerr(err_MPI) if (typeSize*8_MPI_INTEGER_KIND /= int(bit_size(0),MPI_INTEGER_KIND)) & error stop 'Mismatch between MPI_INTEGER and DAMASK default integer' call MPI_Type_size(MPI_INTEGER8,typeSize,err_MPI) - if (err_MPI /= 0) & - error stop 'Could not determine size of MPI_INTEGER8' + call parallelization_chkerr(err_MPI) if (typeSize*8_MPI_INTEGER_KIND /= int(bit_size(0_pI64),MPI_INTEGER_KIND)) & error stop 'Mismatch between MPI_INTEGER8 and DAMASK pI64' call MPI_Type_size(MPI_DOUBLE,typeSize,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) & - error stop 'Could not determine size of MPI_DOUBLE' + call parallelization_chkerr(err_MPI) if (typeSize*8_MPI_INTEGER_KIND /= int(storage_size(0.0_pREAL),MPI_INTEGER_KIND)) & error stop 'Mismatch between MPI_DOUBLE and DAMASK pREAL' + call MPI_Type_size(MPI_INTEGER,typeSize,err_MPI) + call parallelization_chkerr(err_MPI) + if (typeSize*8_MPI_INTEGER_KIND /= int(bit_size(status),MPI_INTEGER_KIND)) & + error stop 'Mismatch between MPI_INTEGER and DAMASK status' + + !$ call get_environment_variable(name='OMP_NUM_THREADS',value=NumThreadsString,STATUS=got_env) !$ if (got_env /= 0) then !$ print'(1x,a)', 'Could not get $OMP_NUM_THREADS, using default' diff --git a/src/phase.f90 b/src/phase.f90 index 39a8679be..ed2f6fa02 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -282,24 +282,22 @@ module phase ! == cleaned:end =================================================================================== - module function phase_thermal_constitutive(Delta_t,ph,en) result(converged_) - + module function phase_thermal_constitutive(Delta_t,ph,en) result(status) real(pREAL), intent(in) :: Delta_t integer, intent(in) :: ph, en - logical :: converged_ - + integer(kind(STATUS_OK)) :: status end function phase_thermal_constitutive - module function phase_damage_constitutive(Delta_t,co,ce) result(converged_) + module function phase_damage_constitutive(Delta_t,co,ce) result(status) real(pREAL), intent(in) :: Delta_t integer, intent(in) :: co, ce - logical :: converged_ + integer(kind(STATUS_OK)) :: status end function phase_damage_constitutive - module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_) + module function phase_mechanical_constitutive(Delta_t,co,ce) result(status) real(pREAL), intent(in) :: Delta_t integer, intent(in) :: co, ce - logical :: converged_ + integer(kind(STATUS_OK)) :: status end function phase_mechanical_constitutive !ToDo: Merge all the stiffness functions diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 2c605c963..ade9e9930 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -120,13 +120,13 @@ end subroutine damage_init !-------------------------------------------------------------------------------------------------- !> @brief calculate stress (P) !-------------------------------------------------------------------------------------------------- -module function phase_damage_constitutive(Delta_t,co,ce) result(converged_) +module function phase_damage_constitutive(Delta_t,co,ce) result(status) real(pREAL), intent(in) :: Delta_t integer, intent(in) :: & co, & ce - logical :: converged_ + integer(kind(STATUS_OK)) :: status integer :: & ph, en @@ -135,7 +135,7 @@ module function phase_damage_constitutive(Delta_t,co,ce) result(converged_) ph = material_ID_phase(co,ce) en = material_entry_phase(co,ce) - converged_ = .not. integrateDamageState(Delta_t,ph,en) + status = integrateDamageState(Delta_t,ph,en) end function phase_damage_constitutive @@ -214,13 +214,13 @@ end function phase_f_phi !> @brief integrate stress, state with adaptive 1st order explicit Euler method !> using Fixed Point Iteration to adapt the stepsize !-------------------------------------------------------------------------------------------------- -function integrateDamageState(Delta_t,ph,en) result(broken) +function integrateDamageState(Delta_t,ph,en) result(status) real(pREAL), intent(in) :: Delta_t integer, intent(in) :: & ph, & en - logical :: broken + integer(kind(STATUS_OK)) :: status integer :: & NiterationState, & !< number of iterations in state loop @@ -235,13 +235,13 @@ function integrateDamageState(Delta_t,ph,en) result(broken) if (damageState(ph)%sizeState == 0) then - broken = .false. + status = STATUS_OK return end if converged_ = .true. - broken = phase_damage_collectDotState(ph,en) - if (broken) return + status = phase_damage_collectDotState(ph,en) + if (status /= STATUS_OK) return size_so = damageState(ph)%sizeDotState damageState(ph)%state(1:size_so,en) = damageState(ph)%state0 (1:size_so,en) & @@ -253,8 +253,8 @@ function integrateDamageState(Delta_t,ph,en) result(broken) if (nIterationState > 1) source_dotState(1:size_so,2) = source_dotState(1:size_so,1) source_dotState(1:size_so,1) = damageState(ph)%dotState(:,en) - broken = phase_damage_collectDotState(ph,en) - if (broken) exit iteration + status = phase_damage_collectDotState(ph,en) + if (status /= STATUS_OK) exit iteration zeta = damper(damageState(ph)%dotState(:,en),source_dotState(1:size_so,1),source_dotState(1:size_so,2)) @@ -270,14 +270,13 @@ function integrateDamageState(Delta_t,ph,en) result(broken) if (converged_) then - broken = phase_damage_deltaState(mechanical_F_e(ph,en),ph,en) + status = phase_damage_deltaState(mechanical_F_e(ph,en),ph,en) exit iteration end if end do iteration - broken = broken .or. .not. converged_ - + if (.not. converged_) status = STATUS_FAIL_PHASE_DAMAGE_STATE contains !-------------------------------------------------------------------------------------------------- @@ -361,15 +360,15 @@ end subroutine damage_result !-------------------------------------------------------------------------------------------------- !> @brief Constitutive equation for calculating the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- -function phase_damage_collectDotState(ph,en) result(broken) +function phase_damage_collectDotState(ph,en) result(status) integer, intent(in) :: & ph, & en !< counter in source loop - logical :: broken + integer(kind(STATUS_OK)) :: status - broken = .false. + status = STATUS_OK if (damageState(ph)%sizeState > 0) then @@ -380,7 +379,7 @@ function phase_damage_collectDotState(ph,en) result(broken) end select sourceType - broken = broken .or. any(IEEE_is_NaN(damageState(ph)%dotState(:,en))) + if (any(IEEE_is_NaN(damageState(ph)%dotState(:,en)))) status = STATUS_FAIL_PHASE_DAMAGE_STATE end if @@ -419,7 +418,7 @@ end function phase_K_phi !> @brief for constitutive models having an instantaneous change of state !> will return false if delta state is not needed/supported by the constitutive model !-------------------------------------------------------------------------------------------------- -function phase_damage_deltaState(Fe, ph, en) result(broken) +function phase_damage_deltaState(Fe, ph, en) result(status) integer, intent(in) :: & ph, & @@ -430,11 +429,10 @@ function phase_damage_deltaState(Fe, ph, en) result(broken) integer :: & myOffset, & mySize - logical :: & - broken + integer(kind(STATUS_OK)) :: status - broken = .false. + status = STATUS_OK if (damageState(ph)%sizeState == 0) return @@ -442,8 +440,8 @@ function phase_damage_deltaState(Fe, ph, en) result(broken) case (DAMAGE_ISOBRITTLE) sourceType call isobrittle_deltaState(phase_homogenizedC66(ph,en), Fe, ph,en) - broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,en))) - if (.not. broken) then + if (any(IEEE_is_NaN(damageState(ph)%deltaState(:,en)))) status = STATUS_FAIL_PHASE_DAMAGE_DELTASTATE + if (status == STATUS_OK) then myOffset = damageState(ph)%offsetDeltaState mySize = damageState(ph)%sizeDeltaState damageState(ph)%state(myOffset + 1: myOffset + mySize,en) = & diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index beb968875..ef4914db4 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -71,12 +71,12 @@ submodule(phase) mechanical dotState end function plastic_dotState - module function plastic_deltaState(ph, en) result(broken) + module function plastic_deltaState(ph, en) result(status) integer, intent(in) :: & ph, & en - logical :: & - broken + integer(kind(STATUS_OK)) :: & + status end function plastic_deltaState module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & @@ -377,11 +377,12 @@ end subroutine mechanical_result !> @brief calculation of stress (P) with time integration based on a residuum in Lp and !> intermediate acceleration of the Newton-Raphson correction !-------------------------------------------------------------------------------------------------- -function integrateStress(F,Fp0,Fi0,Delta_t,ph,en) result(broken) +function integrateStress(F,Fp0,Fi0,Delta_t,ph,en) result(status) real(pREAL), dimension(3,3), intent(in) :: F,Fp0,Fi0 real(pREAL), intent(in) :: Delta_t integer, intent(in) :: ph, en + integer(kind(STATUS_OK)) :: status real(pREAL), dimension(3,3):: Fp_new, & ! plastic deformation gradient at end of timestep invFp_new, & ! inverse of Fp_new @@ -430,10 +431,10 @@ function integrateStress(F,Fp0,Fi0,Delta_t,ph,en) result(broken) p, & jacoCounterLp, & jacoCounterLi ! counters to check for Jacobian update - logical :: error,broken + logical :: error - broken = .true. + status = STATUS_FAIL_PHASE_MECHANICAL_STRESS call plastic_dependentState(ph,en) Lpguess = phase_mechanical_Lp(ph)%data(1:3,1:3,en) ! take as first guess @@ -573,7 +574,7 @@ function integrateStress(F,Fp0,Fi0,Delta_t,ph,en) result(broken) phase_mechanical_Fp(ph)%data(1:3,1:3,en) = Fp_new / math_det33(Fp_new)**(1.0_pREAL/3.0_pREAL) ! regularize phase_mechanical_Fi(ph)%data(1:3,1:3,en) = Fi_new phase_mechanical_Fe(ph)%data(1:3,1:3,en) = matmul(matmul(F,invFp_new),invFi_new) - broken = .false. + status = STATUS_OK end function integrateStress @@ -582,7 +583,7 @@ end function integrateStress !> @brief integrate stress, state with adaptive 1st order explicit Euler method !> using Fixed Point Iteration to adapt the stepsize !-------------------------------------------------------------------------------------------------- -function integrateStateFPI(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken) +function integrateStateFPI(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(status) real(pREAL), intent(in),dimension(3,3) :: F_0,F,Fp0,Fi0 real(pREAL), intent(in),dimension(:) :: state0 @@ -590,8 +591,7 @@ function integrateStateFPI(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken) integer, intent(in) :: & ph, & en - logical :: & - broken + integer(kind(STATUS_OK)) :: status integer :: & NiterationState, & !< number of iterations in state loop @@ -605,7 +605,7 @@ function integrateStateFPI(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken) dotState_last - broken = .true. + status = STATUS_FAIL_PHASE_MECHANICAL_STATE dotState = plastic_dotState(Delta_t,ph,en) if (any(IEEE_is_NaN(dotState))) return @@ -618,8 +618,8 @@ function integrateStateFPI(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken) dotState_last(1:sizeDotState,2) = merge(dotState_last(1:sizeDotState,1),0.0_pREAL, nIterationState > 1) dotState_last(1:sizeDotState,1) = dotState - broken = integrateStress(F,Fp0,Fi0,Delta_t,ph,en) - if (broken) exit iteration + status = integrateStress(F,Fp0,Fi0,Delta_t,ph,en) + if (status /= STATUS_OK) exit iteration dotState = plastic_dotState(Delta_t,ph,en) if (any(IEEE_is_NaN(dotState))) exit iteration @@ -633,7 +633,7 @@ function integrateStateFPI(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken) plasticState(ph)%state(1:sizeDotState,en) = plasticState(ph)%state(1:sizeDotState,en) - r if (converged(r,plasticState(ph)%state(1:sizeDotState,en),plasticState(ph)%atol(1:sizeDotState))) then - broken = plastic_deltaState(ph,en) + status = plastic_deltaState(ph,en) exit iteration end if @@ -670,7 +670,7 @@ end function integrateStateFPI !-------------------------------------------------------------------------------------------------- !> @brief integrate state with 1st order explicit Euler method !-------------------------------------------------------------------------------------------------- -function integrateStateEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken) +function integrateStateEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(status) real(pREAL), intent(in),dimension(3,3) :: F_0,F,Fp0,Fi0 real(pREAL), intent(in),dimension(:) :: state0 @@ -678,8 +678,8 @@ function integrateStateEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken) integer, intent(in) :: & ph, & en !< grain index in grain loop - logical :: & - broken + integer(kind(STATUS_OK)) :: & + status real(pREAL), dimension(plasticState(ph)%sizeDotState) :: & dotState @@ -687,7 +687,7 @@ function integrateStateEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken) sizeDotState - broken = .true. + status = STATUS_FAIL_PHASE_MECHANICAL_STATE dotState = plastic_dotState(Delta_t,ph,en) if (any(IEEE_is_NaN(dotState))) return @@ -695,10 +695,10 @@ function integrateStateEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken) sizeDotState = plasticState(ph)%sizeDotState plasticState(ph)%state(1:sizeDotState,en) = state0 + dotState*Delta_t - broken = plastic_deltaState(ph,en) - if (broken) return + status = plastic_deltaState(ph,en) + if (status /= STATUS_OK) return - broken = integrateStress(F,Fp0,Fi0,Delta_t,ph,en) + status = integrateStress(F,Fp0,Fi0,Delta_t,ph,en) end function integrateStateEuler @@ -706,7 +706,7 @@ end function integrateStateEuler !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, state with 1st order Euler method with adaptive step size !-------------------------------------------------------------------------------------------------- -function integrateStateAdaptiveEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken) +function integrateStateAdaptiveEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(status) real(pREAL), intent(in),dimension(3,3) :: F_0,F,Fp0,Fi0 real(pREAL), intent(in),dimension(:) :: state0 @@ -714,8 +714,8 @@ function integrateStateAdaptiveEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result( integer, intent(in) :: & ph, & en - logical :: & - broken + integer(kind(STATUS_OK)) :: & + status integer :: & sizeDotState @@ -724,7 +724,7 @@ function integrateStateAdaptiveEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result( dotState - broken = .true. + status = STATUS_FAIL_PHASE_MECHANICAL_STATE dotState = plastic_dotState(Delta_t,ph,en) if (any(IEEE_is_NaN(dotState))) return @@ -734,18 +734,21 @@ function integrateStateAdaptiveEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result( r = - dotState * 0.5_pREAL * Delta_t plasticState(ph)%state(1:sizeDotState,en) = state0 + dotState*Delta_t - broken = plastic_deltaState(ph,en) - if (broken) return + status = plastic_deltaState(ph,en) + if (status /= STATUS_OK) return - broken = integrateStress(F,Fp0,Fi0,Delta_t,ph,en) - if (broken) return + status = integrateStress(F,Fp0,Fi0,Delta_t,ph,en) + if (status /= STATUS_OK) return dotState = plastic_dotState(Delta_t,ph,en) + ! ToDo: MD: need to set status to failed if (any(IEEE_is_NaN(dotState))) return - broken = .not. converged(r + 0.5_pREAL * dotState * Delta_t, & + status = merge(STATUS_OK, & + STATUS_FAIL_PHASE_MECHANICAL_STATE, & + converged(r + 0.5_pREAL * dotState * Delta_t, & plasticState(ph)%state(1:sizeDotState,en), & - plasticState(ph)%atol(1:sizeDotState)) + plasticState(ph)%atol(1:sizeDotState))) end function integrateStateAdaptiveEuler @@ -753,13 +756,13 @@ end function integrateStateAdaptiveEuler !--------------------------------------------------------------------------------------------------- !> @brief Integrate state (including stress integration) with the classic Runge Kutta method !--------------------------------------------------------------------------------------------------- -function integrateStateRK4(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken) +function integrateStateRK4(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(status) real(pREAL), intent(in),dimension(3,3) :: F_0,F,Fp0,Fi0 real(pREAL), intent(in),dimension(:) :: state0 real(pREAL), intent(in) :: Delta_t integer, intent(in) :: ph, en - logical :: broken + integer(kind(STATUS_OK)) :: status real(pREAL), dimension(3,3), parameter :: & A = reshape([& @@ -773,7 +776,7 @@ function integrateStateRK4(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken) B = [6.0_pREAL, 3.0_pREAL, 3.0_pREAL, 6.0_pREAL]**(-1) - broken = integrateStateRK(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en,A,B,C) + status = integrateStateRK(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en,A,B,C) end function integrateStateRK4 @@ -781,13 +784,13 @@ end function integrateStateRK4 !--------------------------------------------------------------------------------------------------- !> @brief Integrate state (including stress integration) with the Cash-Carp method !--------------------------------------------------------------------------------------------------- -function integrateStateRKCK45(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken) +function integrateStateRKCK45(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(status) real(pREAL), intent(in),dimension(3,3) :: F_0,F,Fp0,Fi0 real(pREAL), intent(in),dimension(:) :: state0 real(pREAL), intent(in) :: Delta_t integer, intent(in) :: ph, en - logical :: broken + integer(kind(STATUS_OK)) :: status real(pREAL), dimension(5,5), parameter :: & A = reshape([& @@ -808,7 +811,7 @@ function integrateStateRKCK45(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken) 13525.0_pREAL/55296.0_pREAL, 277.0_pREAL/14336.0_pREAL, 1._pREAL/4._pREAL] - broken = integrateStateRK(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en,A,B,C,DB) + status = integrateStateRK(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en,A,B,C,DB) end function integrateStateRKCK45 @@ -817,7 +820,7 @@ end function integrateStateRKCK45 !> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an !! embedded explicit Runge-Kutta method !-------------------------------------------------------------------------------------------------- -function integrateStateRK(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en,A,B,C,DB) result(broken) +function integrateStateRK(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en,A,B,C,DB) result(status) real(pREAL), intent(in),dimension(3,3) :: F_0,F,Fp0,Fi0 real(pREAL), intent(in),dimension(:) :: state0 @@ -828,7 +831,7 @@ function integrateStateRK(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en,A,B,C,DB) result(br integer, intent(in) :: & ph, & en - logical :: broken + integer(kind(STATUS_OK)) :: status integer :: & stage, & ! stage index in integration stage loop @@ -840,7 +843,7 @@ function integrateStateRK(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en,A,B,C,DB) result(br plastic_RKdotState - broken = .true. + status = STATUS_FAIL_PHASE_MECHANICAL_STATE dotState = plastic_dotState(Delta_t,ph,en) if (any(IEEE_is_NaN(dotState))) return @@ -858,14 +861,15 @@ function integrateStateRK(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en,A,B,C,DB) result(br plasticState(ph)%state(1:sizeDotState,en) = state0 + dotState*Delta_t - broken = integrateStress(F_0+(F-F_0)*Delta_t*C(stage),Fp0,Fi0,Delta_t*C(stage), ph,en) - if (broken) exit + status = integrateStress(F_0+(F-F_0)*Delta_t*C(stage),Fp0,Fi0,Delta_t*C(stage), ph,en) + if (status /= STATUS_OK) return dotState = plastic_dotState(Delta_t*C(stage), ph,en) + ! ToDo: MD: need to set status to failed if (any(IEEE_is_NaN(dotState))) exit end do - if (broken) return + if (status /= STATUS_OK) return plastic_RKdotState(1:sizeDotState,size(B)) = dotState @@ -873,16 +877,18 @@ function integrateStateRK(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en,A,B,C,DB) result(br plasticState(ph)%state(1:sizeDotState,en) = state0 + dotState*Delta_t if (present(DB)) & - broken = .not. converged(matmul(plastic_RKdotState(1:sizeDotState,1:size(DB)),DB) * Delta_t, & + status = merge(STATUS_OK, & + STATUS_FAIL_PHASE_MECHANICAL_STATE, & + converged(matmul(plastic_RKdotState(1:sizeDotState,1:size(DB)),DB) * Delta_t, & plasticState(ph)%state(1:sizeDotState,en), & - plasticState(ph)%atol(1:sizeDotState)) + plasticState(ph)%atol(1:sizeDotState))) - if (broken) return + if (status /= STATUS_OK) return - broken = plastic_deltaState(ph,en) - if (broken) return + status = plastic_deltaState(ph,en) + if (status /= STATUS_OK) return - broken = integrateStress(F,Fp0,Fi0,Delta_t,ph,en) + status = integrateStress(F,Fp0,Fi0,Delta_t,ph,en) end function integrateStateRK @@ -984,13 +990,13 @@ end subroutine mechanical_forward !-------------------------------------------------------------------------------------------------- !> @brief calculate stress (P) !-------------------------------------------------------------------------------------------------- -module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_) +module function phase_mechanical_constitutive(Delta_t,co,ce) result(status) real(pREAL), intent(in) :: Delta_t integer, intent(in) :: & co, & ce - logical :: converged_ + integer(kind(STATUS_OK)) :: status real(pREAL) :: & formerStep @@ -1019,13 +1025,13 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_) F0 = phase_mechanical_F0(ph)%data(1:3,1:3,en) stepFrac = 0.0_pREAL todo = .true. - step = 1.0_pREAL/num%stepSizeCryst - converged_ = .false. ! pretend failed step of 1/stepSizeCryst + step = 1.0_pREAL/num%stepSizeCryst ! pretend failed step of 1/stepSizeCryst + status = STATUS_ITERATING todo = .true. cutbackLooping: do while (todo) - if (converged_) then + if (status == STATUS_OK) then formerStep = step stepFrac = stepFrac + step step = min(1.0_pREAL - stepFrac, num%stepIncreaseCryst * step) @@ -1061,7 +1067,7 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_) sizeDotState = plasticState(ph)%sizeDotState F = F0 & + step * (phase_mechanical_F(ph)%data(1:3,1:3,en) - phase_mechanical_F0(ph)%data(1:3,1:3,en)) - converged_ = .not. integrateState(F0,F,Fp0,Fi0,state0(1:sizeDotState),step * Delta_t,ph,en) + status = integrateState(F0,F,Fp0,Fi0,state0(1:sizeDotState),step * Delta_t,ph,en) end if end do cutbackLooping diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index 2d0324742..548957764 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -375,12 +375,12 @@ end subroutine plastic_dependentState !> @brief for constitutive models that have an instantaneous change of state !> will return false if delta state is not needed/supported by the constitutive model !-------------------------------------------------------------------------------------------------- -module function plastic_deltaState(ph, en) result(broken) +module function plastic_deltaState(ph, en) result(status) integer, intent(in) :: & ph, & en - logical :: broken + integer(kind(STATUS_OK)) :: status real(pREAL), dimension(3,3) :: & Mp @@ -388,7 +388,7 @@ module function plastic_deltaState(ph, en) result(broken) mySize - broken = .false. + status = STATUS_OK select case (mechanical_plasticity_type(ph)) case (MECHANICAL_PLASTICITY_NONLOCAL,MECHANICAL_PLASTICITY_KINEHARDENING) @@ -407,8 +407,8 @@ module function plastic_deltaState(ph, en) result(broken) end select plasticType - broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,en))) - if (.not. broken) then + if (any(IEEE_is_NaN(plasticState(ph)%deltaState(:,en)))) status = STATUS_FAIL_PHASE_MECHANICAL_DELTASTATE + if (status == STATUS_OK) then mySize = plasticState(ph)%sizeDeltaState plasticState(ph)%deltaState2(1:mySize,en) = plasticState(ph)%deltaState2(1:mySize,en) & + plasticState(ph)%deltaState(1:mySize,en) diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90 index f63eb5e19..aad7ef470 100644 --- a/src/phase_thermal.f90 +++ b/src/phase_thermal.f90 @@ -187,22 +187,22 @@ end function phase_f_T !-------------------------------------------------------------------------------------------------- !> @brief tbd. !-------------------------------------------------------------------------------------------------- -function phase_thermal_collectDotState(ph,en) result(ok) +function phase_thermal_collectDotState(ph,en) result(status) integer, intent(in) :: ph, en - logical :: ok + integer(kind(STATUS_OK)) :: status integer :: i - ok = .true. + status = STATUS_OK SourceLoop: do i = 1, thermal_Nsources(ph) if (thermal_source_type(i,ph) == THERMAL_SOURCE_EXTERNALHEAT) & call source_externalheat_dotState(ph,en) - ok = ok .and. .not. any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,en))) + if (any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,en)))) status = STATUS_FAIL_PHASE_THERMAL_DOTSTATE end do SourceLoop @@ -238,14 +238,14 @@ module function phase_K_T(co,ce) result(K) end function phase_K_T -module function phase_thermal_constitutive(Delta_t,ph,en) result(converged_) +module function phase_thermal_constitutive(Delta_t,ph,en) result(status) real(pREAL), intent(in) :: Delta_t integer, intent(in) :: ph, en - logical :: converged_ + integer(kind(STATUS_OK)) :: status - converged_ = integrateThermalState(Delta_t,ph,en) + status = integrateThermalState(Delta_t,ph,en) end function phase_thermal_constitutive @@ -253,19 +253,19 @@ end function phase_thermal_constitutive !-------------------------------------------------------------------------------------------------- !> @brief Integrate state with 1st order explicit Euler method. !-------------------------------------------------------------------------------------------------- -function integrateThermalState(Delta_t, ph,en) result(converged) +function integrateThermalState(Delta_t, ph,en) result(status) real(pREAL), intent(in) :: Delta_t integer, intent(in) :: ph, en - logical :: converged + integer(kind(STATUS_OK)) :: status integer :: & so, & sizeDotState - converged = phase_thermal_collectDotState(ph,en) - if (converged) then + status = phase_thermal_collectDotState(ph,en) + if (status == STATUS_OK) then do so = 1, thermal_Nsources(ph) sizeDotState = thermalState(ph)%p(so)%sizeDotState