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

Merge branch 'sophisticated-return-type' into 'development' 

Sophisticated return type

See merge request damask/DAMASK!898
This commit is contained in:
Philip Eisenlohr 2024-01-22 13:47:12 +00:00
parent 9d1b96d3bb db263c3af7
commit 7d8a6d054b
17 changed files with 226 additions and 199 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/Marc/DAMASK_Marc.f90
View File

@ -16,8 +16,8 @@
#endif #endif
#include "../prec.f90" #include "../prec.f90"
#include "../parallelization.f90"
#include "../constants.f90" #include "../constants.f90"
#include "../parallelization.f90"
#include "../misc.f90" #include "../misc.f90"
#include "../IO.f90" #include "../IO.f90"
#include "../YAML_types.f90" #include "../YAML_types.f90"

17
src/Marc/materialpoint_Marc.f90
View File

@ -20,6 +20,7 @@ module materialpoint_Marc
use material use material
use phase use phase
use homogenization use homogenization
use constants
use discretization use discretization
use discretization_Marc use discretization_Marc
@ -27,18 +28,18 @@ module materialpoint_Marc
implicit none(type,external) implicit none(type,external)
private private
real(pREAL), dimension (:,:,:), allocatable, private :: & real(pREAL), dimension (:,:,:), allocatable :: &
materialpoint_cs !< Cauchy stress materialpoint_cs !< Cauchy stress
real(pREAL), dimension (:,:,:,:), allocatable, private :: & real(pREAL), dimension (:,:,:,:), allocatable :: &
materialpoint_dcsdE, & !< Cauchy stress tangent materialpoint_dcsdE, & !< Cauchy stress tangent
materialpoint_F !< deformation gradient materialpoint_F !< deformation gradient
real(pREAL), dimension (:,:,:,:), allocatable, private :: & real(pREAL), dimension (:,:,:,:), allocatable :: &
materialpoint_dcsdE_knownGood !< known good tangent materialpoint_dcsdE_knownGood !< known good tangent
integer, public :: & integer, public :: &
cycleCounter = 0 !< needs description cycleCounter = 0 !< needs description
logical, public :: & integer(kind(STATUS_OK)) :: &
broken = .false. !< needs description status
integer, parameter, public :: & integer, parameter, public :: &
materialpoint_CALCRESULTS = 2**0, & 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 if (iand(mode, materialpoint_CALCRESULTS) /= 0) then
validCalculation: if (broken) then validCalculation: if (status /= STATUS_OK) 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(broken, & call homogenization_mechanical_response(status, &
dt,(elCP-1)*discretization_nIPs + ip, (elCP-1)*discretization_nIPs + ip) 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) call random_number(rnd)
if (rnd < 0.5_pREAL) rnd = rnd - 1.0_pREAL if (rnd < 0.5_pREAL) rnd = rnd - 1.0_pREAL

16
src/constants.f90
View File

@ -1,6 +1,6 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, KU Leuven !> @author Martin Diehl, KU Leuven
!> @brief physical constants !> @brief Constants.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module constants module constants
use prec use prec
@ -20,4 +20,18 @@ module constants
character(len=*), parameter :: LOWER = 'abcdefghijklmnopqrstuvwxyz' character(len=*), parameter :: LOWER = 'abcdefghijklmnopqrstuvwxyz'
character(len=len(LOWER)), parameter :: UPPER = '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 end module constants

26
src/grid/grid_mech_FEM.f90
View File

@ -28,7 +28,7 @@ module grid_mechanical_FEM
use homogenization use homogenization
use discretization use discretization
use discretization_grid use discretization_grid
use constants
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
implicit none(type,external) implicit none(type,external)
@ -84,7 +84,7 @@ module grid_mechanical_FEM
err_BC !< deviation from stress BC err_BC !< deviation from stress BC
integer :: totalIter = 0 !< total iteration in current increment integer :: totalIter = 0 !< total iteration in current increment
logical :: broken integer(kind(STATUS_OK)) :: status
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(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 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,7 +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 = broken solution%termIll = status /= STATUS_OK
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
@ -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) 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. broken) then .or. status /= STATUS_OK) then
reason = 1 reason = 1
elseif (totalIter >= num%itmax) then elseif (totalIter >= num%itmax) then
reason = -1 reason = -1
@ -525,7 +525,7 @@ subroutine formResidual(da_local,x_local, &
real(pREAL), pointer,dimension(:,:,:,:) :: x_scal, r real(pREAL), pointer,dimension(:,:,:,:) :: x_scal, r
real(pREAL), dimension(8,3) :: x_elem, f_elem 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 :: & PetscInt :: &
PETScIter, & PETScIter, &
nfuncs nfuncs
@ -570,10 +570,10 @@ subroutine formResidual(da_local,x_local, &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate constitutive response ! evaluate constitutive response
call utilities_constitutiveResponse(broken,P_current,& call utilities_constitutiveResponse(status,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,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) call parallelization_chkerr(err_MPI)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -587,7 +587,7 @@ subroutine formResidual(da_local,x_local, &
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call DMDAVecGetArrayReadF90(da_local,x_local,x_scal,err_PETSc) call DMDAVecGetArrayReadF90(da_local,x_local,x_scal,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
ele = 0 ce = 0
r = 0.0_pREAL r = 0.0_pREAL
do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, cells(1) do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, cells(1)
ctr = 0 ctr = 0
@ -595,11 +595,11 @@ subroutine formResidual(da_local,x_local, &
ctr = ctr + 1 ctr = ctr + 1
x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk) x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk)
end do; end do; end do 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 + & 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) + & matmul(HGMat,x_elem)*(homogenization_dPdF(1,1,1,1,ce) + &
homogenization_dPdF(2,2,2,2,ele) + & homogenization_dPdF(2,2,2,2,ce) + &
homogenization_dPdF(3,3,3,3,ele))/3.0_pREAL homogenization_dPdF(3,3,3,3,ce))/3.0_pREAL
ctr = 0 ctr = 0
do kk = -1, 0; do jj = -1, 0; do ii = -1, 0 do kk = -1, 0; do jj = -1, 0; do ii = -1, 0
ctr = ctr + 1 ctr = ctr + 1

13
src/grid/grid_mech_spectral_basic.f90
View File

@ -26,6 +26,7 @@ module grid_mechanical_spectral_basic
use grid_mech_utilities use grid_mech_utilities
use homogenization use homogenization
use discretization_grid use discretization_grid
use constants
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
implicit none(type,external) implicit none(type,external)
@ -84,7 +85,7 @@ module grid_mechanical_spectral_basic
err_div !< RMS of div of P 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 integer(kind(STATUS_OK)) :: status
public :: & public :: &
grid_mechanical_spectral_basic_init, & grid_mechanical_spectral_basic_init, &
@ -228,7 +229,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(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 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,7 +288,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 = broken solution%termIll = status /= STATUS_OK
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
@ -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) 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. broken) then .or. status /= STATUS_OK) then
reason = 1 reason = 1
elseif (totalIter >= num%itmax) then elseif (totalIter >= num%itmax) then
reason = -1 reason = -1
@ -514,10 +515,10 @@ subroutine formResidual(residual_subdomain, F, &
end if newIteration end if newIteration
associate (P => r) associate (P => r)
call utilities_constitutiveResponse(broken,P, & call utilities_constitutiveResponse(status,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,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) call parallelization_chkerr(err_MPI)
err_div = utilities_divergenceRMS(P) err_div = utilities_divergenceRMS(P)
end associate end associate

23
src/grid/grid_mech_spectral_polarization.f90
View File

@ -27,6 +27,7 @@ module grid_mechanical_spectral_polarization
use config use config
use homogenization use homogenization
use discretization_grid use discretization_grid
use constants
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
implicit none(type,external) implicit none(type,external)
@ -96,7 +97,7 @@ module grid_mechanical_spectral_polarization
err_div !< RMS of div of P 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 integer(kind(STATUS_OK)) :: status
public :: & public :: &
grid_mechanical_spectral_polarization_init, & grid_mechanical_spectral_polarization_init, &
@ -257,7 +258,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(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 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,7 +323,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 = broken solution%termIll = status /= STATUS_OK
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
@ -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) 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. broken) then .or. status /= STATUS_OK) then
reason = 1 reason = 1
elseif (totalIter >= num%itmax) then elseif (totalIter >= num%itmax) then
reason = -1 reason = -1
@ -562,7 +563,7 @@ subroutine formResidual(residual_subdomain, FandF_tau, &
nfuncs nfuncs
integer(MPI_INTEGER_KIND) :: err_MPI integer(MPI_INTEGER_KIND) :: err_MPI
integer :: & 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) 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) err_curl = utilities_curlRMS(F)
#ifdef __GFORTRAN__ #ifdef __GFORTRAN__
call utilities_constitutiveResponse(broken,r_F, & call utilities_constitutiveResponse(status,r_F, &
#else #else
associate (P => r_F) associate (P => r_F)
call utilities_constitutiveResponse(broken, P, & call utilities_constitutiveResponse(status, 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,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__ #ifdef __GFORTRAN__
err_div = utilities_divergenceRMS(r_F) err_div = utilities_divergenceRMS(r_F)
#else #else
err_div = utilities_divergenceRMS(P) err_div = utilities_divergenceRMS(P)
#endif #endif
e = 0 ce = 0
do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1) 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) = & 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__ #ifdef __GFORTRAN__
r_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), & r_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), &
#else #else

24
src/grid/grid_mech_utilities.f90
View File

@ -18,7 +18,7 @@ module grid_mech_utilities
use discretization use discretization
use spectral_utilities use spectral_utilities
use homogenization use homogenization
use constants
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
implicit none(type,external) implicit none(type,external)
@ -113,10 +113,10 @@ end function utilities_maskedCompliance
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Calculate constitutive response. !> @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) 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,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
@ -124,7 +124,7 @@ subroutine utilities_constitutiveResponse(broken, 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 :: ce
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
@ -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 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 = 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
@ -156,14 +156,14 @@ subroutine utilities_constitutiveResponse(broken, P,P_av,C_volAvg,C_minmaxAvg,&
dPdF_norm_max = 0.0_pREAL dPdF_norm_max = 0.0_pREAL
dPdF_min = huge(1.0_pREAL) dPdF_min = huge(1.0_pREAL)
dPdF_norm_min = huge(1.0_pREAL) dPdF_norm_min = huge(1.0_pREAL)
do i = 1, product(cells(1:2))*cells3 do ce = 1, product(cells(1:2))*cells3
if (dPdF_norm_max < sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)) then 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,i) 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,i)**2) dPdF_norm_max = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,ce)**2)
end if end if
if (dPdF_norm_min > sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)) then 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,i) 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,i)**2) dPdF_norm_min = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,ce)**2)
end if end if
end do end do

7
src/grid/grid_thermal_spectral.f90
View File

@ -25,6 +25,7 @@ module grid_thermal_spectral
use homogenization use homogenization
use YAML_types use YAML_types
use config use config
use constants
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
implicit none(type,external) implicit none(type,external)
@ -56,7 +57,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 integer(kind(STATUS_OK)) :: status
public :: & public :: &
grid_thermal_spectral_init, & grid_thermal_spectral_init, &
@ -207,7 +208,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 .and. .not. broken solution%converged = reason > 0 .and. status == STATUS_OK
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)
@ -324,7 +325,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(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) associate(T => x_scal)
vectorField = utilities_ScalarGradient(T) vectorField = utilities_ScalarGradient(T)

41
src/homogenization.f90
View File

@ -50,9 +50,9 @@ module homogenization
! 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 :: &
homogenization_F !< def grad of IP to be reached at end of FE increment 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 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 homogenization_dPdF !< tangent of first P--K stress at IP
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -210,9 +210,9 @@ end subroutine homogenization_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief !> @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 real(pREAL), intent(in) :: Delta_t !< time increment
integer, intent(in) :: & integer, intent(in) :: &
cell_start, cell_end cell_start, cell_end
@ -224,7 +224,7 @@ subroutine homogenization_mechanical_response(broken,Delta_t,cell_start,cell_end
doneAndHappy doneAndHappy
broken = .false. status = STATUS_OK
!$OMP PARALLEL 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
@ -239,10 +239,10 @@ subroutine homogenization_mechanical_response(broken,Delta_t,cell_start,cell_end
doneAndHappy = [.false.,.true.] 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) 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 if (converged) then
doneAndHappy = mechanical_updateState(Delta_t,homogenization_F(1:3,1:3,ce),ce) doneAndHappy = mechanical_updateState(Delta_t,homogenization_F(1:3,1:3,ce),ce)
converged = all(doneAndHappy) converged = all(doneAndHappy)
@ -251,19 +251,19 @@ subroutine homogenization_mechanical_response(broken,Delta_t,cell_start,cell_end
end if end if
end do convergenceLooping end do convergenceLooping
if (.not. converged) then if (.not. converged) then
if (.not. broken) print*, ' Cell ', ce, ' failed (mechanics)' if (status == STATUS_OK) print*, ' Cell ', ce, ' failed (mechanics)'
broken = .true. status = STATUS_FAIL_PHASE_MECHANICAL
end if 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. converged) then
if (.not. broken) print*, ' Cell ', ce, ' failed (damage)' if (status == STATUS_OK) print*, ' Cell ', ce, ' failed (damage)'
broken = .true. status = STATUS_FAIL_PHASE_DAMAGE
end if end if
end do end do
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
if (broken) return if (status /= STATUS_OK) return
!$OMP PARALLEL DO PRIVATE(ho) !$OMP PARALLEL DO PRIVATE(ho)
do ce = cell_start, cell_end do ce = cell_start, cell_end
@ -281,10 +281,10 @@ end subroutine homogenization_mechanical_response
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief !> @brief
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine homogenization_thermal_response(broken, & subroutine homogenization_thermal_response(status, &
Delta_t,cell_start,cell_end) 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 real(pREAL), intent(in) :: Delta_t !< time increment
integer, intent(in) :: & integer, intent(in) :: &
cell_start, cell_end cell_start, cell_end
@ -293,20 +293,19 @@ subroutine homogenization_thermal_response(broken, &
co, ce, ho co, ce, ho
broken = .false. status = STATUS_OK
!$OMP PARALLEL DO PRIVATE(ho) !$OMP PARALLEL DO PRIVATE(ho)
do ce = cell_start, cell_end do ce = cell_start, cell_end
if (broken) continue if (status /= STATUS_OK) 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 (phase_thermal_constitutive(Delta_t,material_ID_phase(co,ce),material_entry_phase(co,ce)) /= STATUS_OK) then
if (.not. broken) print*, ' Cell ', ce, ' failed (thermal)' if (status == STATUS_OK) print*, ' Cell ', ce, ' failed (thermal)'
broken = .true. status = STATUS_FAIL_PHASE_THERMAL
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

7
src/mesh/FEM_utilities.f90
View File

@ -22,6 +22,7 @@ module FEM_utilities
use discretization_mesh use discretization_mesh
use homogenization use homogenization
use FEM_quadrature use FEM_quadrature
use constants
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
implicit none(type,external) implicit none(type,external)
@ -120,9 +121,9 @@ end subroutine FEM_utilities_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates constitutive response !> @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 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
@ -132,7 +133,7 @@ subroutine utilities_constitutiveResponse(broken, Delta_t,P_av,forwardData)
print'(/,1x,a)', '... evaluating constitutive response ......................................' 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. cutBack = .false.
P_av = sum(homogenization_P,dim=3) * wgt P_av = sum(homogenization_P,dim=3) * wgt

22
src/mesh/mesh_mech_FEM.f90
View File

@ -25,6 +25,7 @@ module mesh_mechanical_FEM
use FEM_quadrature use FEM_quadrature
use homogenization use homogenization
use math use math
use constants
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
implicit none(type,external) implicit none(type,external)
@ -68,7 +69,8 @@ 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, broken logical :: ForwardData
integer(kind(STATUS_OK)) :: status
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 +313,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(broken,0.0_pREAL,devNull,.true.) call utilities_constitutiveResponse(status,0.0_pREAL,devNull,.true.)
end subroutine FEM_mechanical_init 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 ! evaluate constitutive response
call utilities_constitutiveResponse(broken,params%Delta_t,P_av,ForwardData) call utilities_constitutiveResponse(status,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 MPI_Allreduce(MPI_IN_PLACE,status,1_MPI_INTEGER_KIND,MPI_INTEGER,MPI_MAX,MPI_COMM_WORLD,err_MPI)
call parallelization_chkerr(err_MPI) call parallelization_chkerr(err_MPI)
ForwardData = .false. 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 real(pREAL),dimension(cellDOF,cellDOF) :: K_eA, K_eB
PetscInt :: cellStart, cellEnd, cell, component, face, & PetscInt :: cellStart, cellEnd, cell, component, face, &
qPt, basis, comp, cidx,bcSize, m, i qPt, basis, comp, cidx,bcSize, ce, i
IS :: bcPoints IS :: bcPoints
@ -581,7 +583,7 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P
FAvg = 0.0_pREAL FAvg = 0.0_pREAL
BMatAvg = 0.0_pREAL BMatAvg = 0.0_pREAL
do qPt = 0_pPETSCINT, nQuadrature-1_pPETSCINT do qPt = 0_pPETSCINT, nQuadrature-1_pPETSCINT
m = cell*nQuadrature + qPt + 1_pPETSCINT ce = cell*nQuadrature + qPt + 1_pPETSCINT
BMat = 0.0_pREAL BMat = 0.0_pREAL
do basis = 0_pPETSCINT, nBasis-1_pPETSCINT do basis = 0_pPETSCINT, nBasis-1_pPETSCINT
do comp = 0_pPETSCINT, dimPlex-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)) matmul(reshape(pInvcellJ,[dimPlex,dimPlex]),basisFieldDer(i*dimPlex+1_pPETSCINT:(i+1_pPETSCINT)*dimPlex))
end do end do
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], order=[2,1,4,3]), &
shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1_pPETSCINT) shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1_pPETSCINT)
if (num%BBarStabilization) then 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) FInv = math_inv33(F)
K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0_pREAL/real(dimPlex,pREAL)) K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0_pREAL/real(dimPlex,pREAL))
K_eB = K_eB - & 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), & matmul(reshape(FInv(1:dimPlex,1:dimPlex), &
shape=[1_pPETSCINT,dimPlex**2],order=[2,1]),BMat))),MatA) shape=[1_pPETSCINT,dimPlex**2],order=[2,1]),BMat))),MatA)
MatB = MatB & 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 FAvg = FAvg + F
BMatAvg = BMatAvg + BMat BMatAvg = BMatAvg + BMat
else 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) 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 (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), & 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))', &

21
src/parallelization.f90
View File

@ -7,6 +7,8 @@ module parallelization
OUTPUT_UNIT, & OUTPUT_UNIT, &
ERROR_UNIT ERROR_UNIT
use constants
#ifdef PETSC #ifdef PETSC
#include <petsc/finclude/petscsys.h> #include <petsc/finclude/petscsys.h>
use PETScSys use PETScSys
@ -63,6 +65,7 @@ subroutine parallelization_init()
!$ integer(pI32) :: OMP_NUM_THREADS !$ integer(pI32) :: OMP_NUM_THREADS
!$ character(len=6) NumThreadsString !$ character(len=6) NumThreadsString
PetscErrorCode :: err_PETSc PetscErrorCode :: err_PETSc
integer(kind(STATUS_OK)) :: status
#ifdef _OPENMP #ifdef _OPENMP
@ -116,28 +119,30 @@ subroutine parallelization_init()
#endif #endif
call MPI_Comm_size(MPI_COMM_WORLD,worldsize,err_MPI) call MPI_Comm_size(MPI_COMM_WORLD,worldsize,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) & call parallelization_chkerr(err_MPI)
error stop 'Could not determine worldsize'
if (worldrank == 0) print'(/,1x,a,i0)', 'MPI processes: ',worldsize if (worldrank == 0) print'(/,1x,a,i0)', 'MPI processes: ',worldsize
call MPI_Type_size(MPI_INTEGER,typeSize,err_MPI) call MPI_Type_size(MPI_INTEGER,typeSize,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) & call parallelization_chkerr(err_MPI)
error stop 'Could not determine size of MPI_INTEGER'
if (typeSize*8_MPI_INTEGER_KIND /= int(bit_size(0),MPI_INTEGER_KIND)) & if (typeSize*8_MPI_INTEGER_KIND /= int(bit_size(0),MPI_INTEGER_KIND)) &
error stop 'Mismatch between MPI_INTEGER and DAMASK default integer' error stop 'Mismatch between MPI_INTEGER and DAMASK default integer'
call MPI_Type_size(MPI_INTEGER8,typeSize,err_MPI) call MPI_Type_size(MPI_INTEGER8,typeSize,err_MPI)
if (err_MPI /= 0) & call parallelization_chkerr(err_MPI)
error stop 'Could not determine size of MPI_INTEGER8'
if (typeSize*8_MPI_INTEGER_KIND /= int(bit_size(0_pI64),MPI_INTEGER_KIND)) & if (typeSize*8_MPI_INTEGER_KIND /= int(bit_size(0_pI64),MPI_INTEGER_KIND)) &
error stop 'Mismatch between MPI_INTEGER8 and DAMASK pI64' error stop 'Mismatch between MPI_INTEGER8 and DAMASK pI64'
call MPI_Type_size(MPI_DOUBLE,typeSize,err_MPI) call MPI_Type_size(MPI_DOUBLE,typeSize,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) & call parallelization_chkerr(err_MPI)
error stop 'Could not determine size of MPI_DOUBLE'
if (typeSize*8_MPI_INTEGER_KIND /= int(storage_size(0.0_pREAL),MPI_INTEGER_KIND)) & if (typeSize*8_MPI_INTEGER_KIND /= int(storage_size(0.0_pREAL),MPI_INTEGER_KIND)) &
error stop 'Mismatch between MPI_DOUBLE and DAMASK pREAL' 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) !$ call get_environment_variable(name='OMP_NUM_THREADS',value=NumThreadsString,STATUS=got_env)
!$ if (got_env /= 0) then !$ if (got_env /= 0) then
!$ print'(1x,a)', 'Could not get $OMP_NUM_THREADS, using default' !$ print'(1x,a)', 'Could not get $OMP_NUM_THREADS, using default'

14
src/phase.f90
View File

@ -282,24 +282,22 @@ module phase
! == cleaned:end =================================================================================== ! == 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 real(pREAL), intent(in) :: Delta_t
integer, intent(in) :: ph, en integer, intent(in) :: ph, en
logical :: converged_ integer(kind(STATUS_OK)) :: status
end function phase_thermal_constitutive 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 real(pREAL), intent(in) :: Delta_t
integer, intent(in) :: co, ce integer, intent(in) :: co, ce
logical :: converged_ integer(kind(STATUS_OK)) :: status
end function phase_damage_constitutive 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 real(pREAL), intent(in) :: Delta_t
integer, intent(in) :: co, ce integer, intent(in) :: co, ce
logical :: converged_ integer(kind(STATUS_OK)) :: status
end function phase_mechanical_constitutive end function phase_mechanical_constitutive
!ToDo: Merge all the stiffness functions !ToDo: Merge all the stiffness functions

44
src/phase_damage.f90
View File

@ -120,13 +120,13 @@ end subroutine damage_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculate stress (P) !> @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 real(pREAL), intent(in) :: Delta_t
integer, intent(in) :: & integer, intent(in) :: &
co, & co, &
ce ce
logical :: converged_ integer(kind(STATUS_OK)) :: status
integer :: & integer :: &
ph, en ph, en
@ -135,7 +135,7 @@ module function phase_damage_constitutive(Delta_t,co,ce) result(converged_)
ph = material_ID_phase(co,ce) ph = material_ID_phase(co,ce)
en = material_entry_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 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 !> @brief integrate stress, state with adaptive 1st order explicit Euler method
!> using Fixed Point Iteration to adapt the stepsize !> 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 real(pREAL), intent(in) :: Delta_t
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
en en
logical :: broken integer(kind(STATUS_OK)) :: status
integer :: & integer :: &
NiterationState, & !< number of iterations in state loop 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 if (damageState(ph)%sizeState == 0) then
broken = .false. status = STATUS_OK
return return
end if end if
converged_ = .true. converged_ = .true.
broken = phase_damage_collectDotState(ph,en) status = phase_damage_collectDotState(ph,en)
if (broken) return if (status /= STATUS_OK) return
size_so = damageState(ph)%sizeDotState size_so = damageState(ph)%sizeDotState
damageState(ph)%state(1:size_so,en) = damageState(ph)%state0 (1:size_so,en) & 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) 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) source_dotState(1:size_so,1) = damageState(ph)%dotState(:,en)
broken = phase_damage_collectDotState(ph,en) status = phase_damage_collectDotState(ph,en)
if (broken) exit iteration if (status /= STATUS_OK) exit iteration
zeta = damper(damageState(ph)%dotState(:,en),source_dotState(1:size_so,1),source_dotState(1:size_so,2)) 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 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 exit iteration
end if end if
end do iteration end do iteration
broken = broken .or. .not. converged_ if (.not. converged_) status = STATUS_FAIL_PHASE_DAMAGE_STATE
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -361,15 +360,15 @@ end subroutine damage_result
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Constitutive equation for calculating the rate of change of microstructure. !> @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) :: & integer, intent(in) :: &
ph, & ph, &
en !< counter in source loop en !< counter in source loop
logical :: broken integer(kind(STATUS_OK)) :: status
broken = .false. status = STATUS_OK
if (damageState(ph)%sizeState > 0) then if (damageState(ph)%sizeState > 0) then
@ -380,7 +379,7 @@ function phase_damage_collectDotState(ph,en) result(broken)
end select sourceType 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 end if
@ -419,7 +418,7 @@ end function phase_K_phi
!> @brief for constitutive models having an instantaneous change of state !> @brief for constitutive models having an instantaneous change of state
!> will return false if delta state is not needed/supported by the constitutive model !> 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) :: & integer, intent(in) :: &
ph, & ph, &
@ -430,11 +429,10 @@ function phase_damage_deltaState(Fe, ph, en) result(broken)
integer :: & integer :: &
myOffset, & myOffset, &
mySize mySize
logical :: & integer(kind(STATUS_OK)) :: status
broken
broken = .false. status = STATUS_OK
if (damageState(ph)%sizeState == 0) return if (damageState(ph)%sizeState == 0) return
@ -442,8 +440,8 @@ function phase_damage_deltaState(Fe, ph, en) result(broken)
case (DAMAGE_ISOBRITTLE) sourceType case (DAMAGE_ISOBRITTLE) sourceType
call isobrittle_deltaState(phase_homogenizedC66(ph,en), Fe, ph,en) call isobrittle_deltaState(phase_homogenizedC66(ph,en), Fe, ph,en)
broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,en))) if (any(IEEE_is_NaN(damageState(ph)%deltaState(:,en)))) status = STATUS_FAIL_PHASE_DAMAGE_DELTASTATE
if (.not. broken) then if (status == STATUS_OK) then
myOffset = damageState(ph)%offsetDeltaState myOffset = damageState(ph)%offsetDeltaState
mySize = damageState(ph)%sizeDeltaState mySize = damageState(ph)%sizeDeltaState
damageState(ph)%state(myOffset + 1: myOffset + mySize,en) = & damageState(ph)%state(myOffset + 1: myOffset + mySize,en) = &

116
src/phase_mechanical.f90
View File

@ -71,12 +71,12 @@ submodule(phase) mechanical
dotState dotState
end function plastic_dotState end function plastic_dotState
module function plastic_deltaState(ph, en) result(broken) module function plastic_deltaState(ph, en) result(status)
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
en en
logical :: & integer(kind(STATUS_OK)) :: &
broken status
end function plastic_deltaState end function plastic_deltaState
module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & 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 !> @brief calculation of stress (P) with time integration based on a residuum in Lp and
!> intermediate acceleration of the Newton-Raphson correction !> 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), dimension(3,3), intent(in) :: F,Fp0,Fi0
real(pREAL), intent(in) :: Delta_t real(pREAL), intent(in) :: Delta_t
integer, intent(in) :: ph, en integer, intent(in) :: ph, en
integer(kind(STATUS_OK)) :: status
real(pREAL), dimension(3,3):: Fp_new, & ! plastic deformation gradient at end of timestep real(pREAL), dimension(3,3):: Fp_new, & ! plastic deformation gradient at end of timestep
invFp_new, & ! inverse of Fp_new invFp_new, & ! inverse of Fp_new
@ -430,10 +431,10 @@ function integrateStress(F,Fp0,Fi0,Delta_t,ph,en) result(broken)
p, & p, &
jacoCounterLp, & jacoCounterLp, &
jacoCounterLi ! counters to check for Jacobian update 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) call plastic_dependentState(ph,en)
Lpguess = phase_mechanical_Lp(ph)%data(1:3,1:3,en) ! take as first guess 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_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_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) 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 end function integrateStress
@ -582,7 +583,7 @@ end function integrateStress
!> @brief integrate stress, state with adaptive 1st order explicit Euler method !> @brief integrate stress, state with adaptive 1st order explicit Euler method
!> using Fixed Point Iteration to adapt the stepsize !> 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(3,3) :: F_0,F,Fp0,Fi0
real(pREAL), intent(in),dimension(:) :: state0 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) :: & integer, intent(in) :: &
ph, & ph, &
en en
logical :: & integer(kind(STATUS_OK)) :: status
broken
integer :: & integer :: &
NiterationState, & !< number of iterations in state loop 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 dotState_last
broken = .true. status = STATUS_FAIL_PHASE_MECHANICAL_STATE
dotState = plastic_dotState(Delta_t,ph,en) dotState = plastic_dotState(Delta_t,ph,en)
if (any(IEEE_is_NaN(dotState))) return 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,2) = merge(dotState_last(1:sizeDotState,1),0.0_pREAL, nIterationState > 1)
dotState_last(1:sizeDotState,1) = dotState dotState_last(1:sizeDotState,1) = dotState
broken = integrateStress(F,Fp0,Fi0,Delta_t,ph,en) status = integrateStress(F,Fp0,Fi0,Delta_t,ph,en)
if (broken) exit iteration if (status /= STATUS_OK) exit iteration
dotState = plastic_dotState(Delta_t,ph,en) dotState = plastic_dotState(Delta_t,ph,en)
if (any(IEEE_is_NaN(dotState))) exit iteration 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 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 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 exit iteration
end if end if
@ -670,7 +670,7 @@ end function integrateStateFPI
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief integrate state with 1st order explicit Euler method !> @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(3,3) :: F_0,F,Fp0,Fi0
real(pREAL), intent(in),dimension(:) :: state0 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) :: & integer, intent(in) :: &
ph, & ph, &
en !< grain index in grain loop en !< grain index in grain loop
logical :: & integer(kind(STATUS_OK)) :: &
broken status
real(pREAL), dimension(plasticState(ph)%sizeDotState) :: & real(pREAL), dimension(plasticState(ph)%sizeDotState) :: &
dotState dotState
@ -687,7 +687,7 @@ function integrateStateEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(broken)
sizeDotState sizeDotState
broken = .true. status = STATUS_FAIL_PHASE_MECHANICAL_STATE
dotState = plastic_dotState(Delta_t,ph,en) dotState = plastic_dotState(Delta_t,ph,en)
if (any(IEEE_is_NaN(dotState))) return 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 sizeDotState = plasticState(ph)%sizeDotState
plasticState(ph)%state(1:sizeDotState,en) = state0 + dotState*Delta_t plasticState(ph)%state(1:sizeDotState,en) = state0 + dotState*Delta_t
broken = plastic_deltaState(ph,en) status = plastic_deltaState(ph,en)
if (broken) return 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 end function integrateStateEuler
@ -706,7 +706,7 @@ end function integrateStateEuler
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief integrate stress, state with 1st order Euler method with adaptive step size !> @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(3,3) :: F_0,F,Fp0,Fi0
real(pREAL), intent(in),dimension(:) :: state0 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) :: & integer, intent(in) :: &
ph, & ph, &
en en
logical :: & integer(kind(STATUS_OK)) :: &
broken status
integer :: & integer :: &
sizeDotState sizeDotState
@ -724,7 +724,7 @@ function integrateStateAdaptiveEuler(F_0,F,Fp0,Fi0,state0,Delta_t,ph,en) result(
dotState dotState
broken = .true. status = STATUS_FAIL_PHASE_MECHANICAL_STATE
dotState = plastic_dotState(Delta_t,ph,en) dotState = plastic_dotState(Delta_t,ph,en)
if (any(IEEE_is_NaN(dotState))) return 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 r = - dotState * 0.5_pREAL * Delta_t
plasticState(ph)%state(1:sizeDotState,en) = state0 + dotState*Delta_t plasticState(ph)%state(1:sizeDotState,en) = state0 + dotState*Delta_t
broken = plastic_deltaState(ph,en) status = plastic_deltaState(ph,en)
if (broken) return if (status /= STATUS_OK) return
broken = integrateStress(F,Fp0,Fi0,Delta_t,ph,en) status = integrateStress(F,Fp0,Fi0,Delta_t,ph,en)
if (broken) return if (status /= STATUS_OK) return
dotState = plastic_dotState(Delta_t,ph,en) dotState = plastic_dotState(Delta_t,ph,en)
! ToDo: MD: need to set status to failed
if (any(IEEE_is_NaN(dotState))) return 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)%state(1:sizeDotState,en), &
plasticState(ph)%atol(1:sizeDotState)) plasticState(ph)%atol(1:sizeDotState)))
end function integrateStateAdaptiveEuler end function integrateStateAdaptiveEuler
@ -753,13 +756,13 @@ end function integrateStateAdaptiveEuler
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> @brief Integrate state (including stress integration) with the classic Runge Kutta method !> @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(3,3) :: F_0,F,Fp0,Fi0
real(pREAL), intent(in),dimension(:) :: state0 real(pREAL), intent(in),dimension(:) :: state0
real(pREAL), intent(in) :: Delta_t real(pREAL), intent(in) :: Delta_t
integer, intent(in) :: ph, en integer, intent(in) :: ph, en
logical :: broken integer(kind(STATUS_OK)) :: status
real(pREAL), dimension(3,3), parameter :: & real(pREAL), dimension(3,3), parameter :: &
A = reshape([& 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) 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 end function integrateStateRK4
@ -781,13 +784,13 @@ end function integrateStateRK4
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> @brief Integrate state (including stress integration) with the Cash-Carp method !> @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(3,3) :: F_0,F,Fp0,Fi0
real(pREAL), intent(in),dimension(:) :: state0 real(pREAL), intent(in),dimension(:) :: state0
real(pREAL), intent(in) :: Delta_t real(pREAL), intent(in) :: Delta_t
integer, intent(in) :: ph, en integer, intent(in) :: ph, en
logical :: broken integer(kind(STATUS_OK)) :: status
real(pREAL), dimension(5,5), parameter :: & real(pREAL), dimension(5,5), parameter :: &
A = reshape([& 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] 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 end function integrateStateRKCK45
@ -817,7 +820,7 @@ end function integrateStateRKCK45
!> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an !> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an
!! embedded explicit Runge-Kutta method !! 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(3,3) :: F_0,F,Fp0,Fi0
real(pREAL), intent(in),dimension(:) :: state0 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) :: & integer, intent(in) :: &
ph, & ph, &
en en
logical :: broken integer(kind(STATUS_OK)) :: status
integer :: & integer :: &
stage, & ! stage index in integration stage loop 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 plastic_RKdotState
broken = .true. status = STATUS_FAIL_PHASE_MECHANICAL_STATE
dotState = plastic_dotState(Delta_t,ph,en) dotState = plastic_dotState(Delta_t,ph,en)
if (any(IEEE_is_NaN(dotState))) return 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 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) status = integrateStress(F_0+(F-F_0)*Delta_t*C(stage),Fp0,Fi0,Delta_t*C(stage), ph,en)
if (broken) exit if (status /= STATUS_OK) return
dotState = plastic_dotState(Delta_t*C(stage), ph,en) dotState = plastic_dotState(Delta_t*C(stage), ph,en)
! ToDo: MD: need to set status to failed
if (any(IEEE_is_NaN(dotState))) exit if (any(IEEE_is_NaN(dotState))) exit
end do end do
if (broken) return if (status /= STATUS_OK) return
plastic_RKdotState(1:sizeDotState,size(B)) = dotState 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 plasticState(ph)%state(1:sizeDotState,en) = state0 + dotState*Delta_t
if (present(DB)) & 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)%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) status = plastic_deltaState(ph,en)
if (broken) return 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 end function integrateStateRK
@ -984,13 +990,13 @@ end subroutine mechanical_forward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculate stress (P) !> @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 real(pREAL), intent(in) :: Delta_t
integer, intent(in) :: & integer, intent(in) :: &
co, & co, &
ce ce
logical :: converged_ integer(kind(STATUS_OK)) :: status
real(pREAL) :: & real(pREAL) :: &
formerStep 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) F0 = phase_mechanical_F0(ph)%data(1:3,1:3,en)
stepFrac = 0.0_pREAL stepFrac = 0.0_pREAL
todo = .true. todo = .true.
step = 1.0_pREAL/num%stepSizeCryst step = 1.0_pREAL/num%stepSizeCryst ! pretend failed step of 1/stepSizeCryst
converged_ = .false. ! pretend failed step of 1/stepSizeCryst status = STATUS_ITERATING
todo = .true. todo = .true.
cutbackLooping: do while (todo) cutbackLooping: do while (todo)
if (converged_) then if (status == STATUS_OK) then
formerStep = step formerStep = step
stepFrac = stepFrac + step stepFrac = stepFrac + step
step = min(1.0_pREAL - stepFrac, num%stepIncreaseCryst * 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 sizeDotState = plasticState(ph)%sizeDotState
F = F0 & F = F0 &
+ step * (phase_mechanical_F(ph)%data(1:3,1:3,en) - phase_mechanical_F0(ph)%data(1:3,1:3,en)) + 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 if
end do cutbackLooping end do cutbackLooping

10
src/phase_mechanical_plastic.f90
View File

@ -375,12 +375,12 @@ end subroutine plastic_dependentState
!> @brief for constitutive models that have an instantaneous change of state !> @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 !> 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) :: & integer, intent(in) :: &
ph, & ph, &
en en
logical :: broken integer(kind(STATUS_OK)) :: status
real(pREAL), dimension(3,3) :: & real(pREAL), dimension(3,3) :: &
Mp Mp
@ -388,7 +388,7 @@ module function plastic_deltaState(ph, en) result(broken)
mySize mySize
broken = .false. status = STATUS_OK
select case (mechanical_plasticity_type(ph)) select case (mechanical_plasticity_type(ph))
case (MECHANICAL_PLASTICITY_NONLOCAL,MECHANICAL_PLASTICITY_KINEHARDENING) case (MECHANICAL_PLASTICITY_NONLOCAL,MECHANICAL_PLASTICITY_KINEHARDENING)
@ -407,8 +407,8 @@ module function plastic_deltaState(ph, en) result(broken)
end select plasticType end select plasticType
broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,en))) if (any(IEEE_is_NaN(plasticState(ph)%deltaState(:,en)))) status = STATUS_FAIL_PHASE_MECHANICAL_DELTASTATE
if (.not. broken) then if (status == STATUS_OK) then
mySize = plasticState(ph)%sizeDeltaState mySize = plasticState(ph)%sizeDeltaState
plasticState(ph)%deltaState2(1:mySize,en) = plasticState(ph)%deltaState2(1:mySize,en) & plasticState(ph)%deltaState2(1:mySize,en) = plasticState(ph)%deltaState2(1:mySize,en) &
+ plasticState(ph)%deltaState(1:mySize,en) + plasticState(ph)%deltaState(1:mySize,en)

22
src/phase_thermal.f90
View File

@ -187,22 +187,22 @@ end function phase_f_T
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief tbd. !> @brief tbd.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function phase_thermal_collectDotState(ph,en) result(ok) function phase_thermal_collectDotState(ph,en) result(status)
integer, intent(in) :: ph, en integer, intent(in) :: ph, en
logical :: ok integer(kind(STATUS_OK)) :: status
integer :: i integer :: i
ok = .true. status = STATUS_OK
SourceLoop: do i = 1, thermal_Nsources(ph) SourceLoop: do i = 1, thermal_Nsources(ph)
if (thermal_source_type(i,ph) == THERMAL_SOURCE_EXTERNALHEAT) & if (thermal_source_type(i,ph) == THERMAL_SOURCE_EXTERNALHEAT) &
call source_externalheat_dotState(ph,en) 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 end do SourceLoop
@ -238,14 +238,14 @@ module function phase_K_T(co,ce) result(K)
end function phase_K_T 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 real(pREAL), intent(in) :: Delta_t
integer, intent(in) :: ph, en 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 end function phase_thermal_constitutive
@ -253,19 +253,19 @@ end function phase_thermal_constitutive
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Integrate state with 1st order explicit Euler method. !> @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 real(pREAL), intent(in) :: Delta_t
integer, intent(in) :: ph, en integer, intent(in) :: ph, en
logical :: converged integer(kind(STATUS_OK)) :: status
integer :: & integer :: &
so, & so, &
sizeDotState sizeDotState
converged = phase_thermal_collectDotState(ph,en) status = phase_thermal_collectDotState(ph,en)
if (converged) then if (status == STATUS_OK) then
do so = 1, thermal_Nsources(ph) do so = 1, thermal_Nsources(ph)
sizeDotState = thermalState(ph)%p(so)%sizeDotState sizeDotState = thermalState(ph)%p(so)%sizeDotState