DAMASK_EICMD/src/CPFEM.f90

282 lines
12 KiB
Fortran
Raw Normal View History

2013-03-06 20:11:15 +05:30
!--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief CPFEM engine
!--------------------------------------------------------------------------------------------------
module CPFEM
use DAMASK_interface
use prec
use IO
2020-04-28 13:59:13 +05:30
use YAML_types
2020-05-22 00:11:40 +05:30
use YAML_parse
use HDF5_utilities
use results
use config
use math
use rotations
2022-01-31 19:35:15 +05:30
use polynomials
use lattice
use material
use phase
use homogenization
use discretization
use discretization_marc
implicit none
private
2020-02-29 21:46:40 +05:30
real(pReal), dimension (:,:,:), allocatable, private :: &
CPFEM_cs !< Cauchy stress
real(pReal), dimension (:,:,:,:), allocatable, private :: &
CPFEM_dcsdE !< Cauchy stress tangent
real(pReal), dimension (:,:,:,:), allocatable, private :: &
CPFEM_dcsdE_knownGood !< known good tangent
2020-06-16 10:32:29 +05:30
integer, public :: &
cycleCounter = 0 !< needs description
integer, parameter, public :: &
CPFEM_CALCRESULTS = 2**0, &
CPFEM_AGERESULTS = 2**1, &
CPFEM_BACKUPJACOBIAN = 2**2, &
CPFEM_RESTOREJACOBIAN = 2**3
type, private :: tNumerics
integer :: &
iJacoStiffness !< frequency of stiffness update
end type tNumerics
2020-06-29 18:39:13 +05:30
type(tNumerics), private :: num
2020-09-13 13:49:38 +05:30
2020-07-01 23:24:14 +05:30
type, private :: tDebugOptions
logical :: &
basic, &
extensive, &
selective
2020-07-01 23:24:14 +05:30
integer:: &
element, &
ip
end type tDebugOptions
2020-09-13 13:49:38 +05:30
type(tDebugOptions), private :: debugCPFEM
2020-09-13 13:49:38 +05:30
public :: &
CPFEM_general, &
CPFEM_initAll, &
CPFEM_results
contains
2013-03-06 20:11:15 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Initialize all modules.
2013-03-06 20:11:15 +05:30
!--------------------------------------------------------------------------------------------------
2020-06-11 12:21:17 +05:30
subroutine CPFEM_initAll
2020-04-27 17:10:22 +05:30
call DAMASK_interface_init
2020-09-13 13:49:38 +05:30
call prec_init
2020-04-27 17:10:22 +05:30
call IO_init
2020-09-29 23:37:33 +05:30
call YAML_types_init
call YAML_parse_init
call HDF5_utilities_init
call results_init(.false.)
2020-04-27 17:10:22 +05:30
call config_init
call math_init
call rotations_init
2022-01-31 19:35:15 +05:30
call polynomials_init
2020-04-27 17:10:22 +05:30
call lattice_init
call discretization_marc_init
call material_init(.false.)
2021-02-09 03:51:53 +05:30
call phase_init
2020-04-27 17:10:22 +05:30
call homogenization_init
call CPFEM_init
2020-08-15 19:32:10 +05:30
call config_deallocate
end subroutine CPFEM_initAll
2013-03-06 20:11:15 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief allocate the arrays defined in module CPFEM and initialize them
!--------------------------------------------------------------------------------------------------
subroutine CPFEM_init
class(tNode), pointer :: &
debug_CPFEM
print'(/,1x,a)', '<<<+- CPFEM init -+>>>'; flush(IO_STDOUT)
allocate(CPFEM_cs( 6,discretization_nIPs,discretization_Nelems), source= 0.0_pReal)
allocate(CPFEM_dcsdE( 6,6,discretization_nIPs,discretization_Nelems), source= 0.0_pReal)
allocate(CPFEM_dcsdE_knownGood(6,6,discretization_nIPs,discretization_Nelems), source= 0.0_pReal)
!------------------------------------------------------------------------------
2020-09-13 13:49:38 +05:30
! read debug options
2021-03-26 22:07:38 +05:30
debug_CPFEM => config_debug%get('CPFEM',defaultVal=emptyList)
debugCPFEM%basic = debug_CPFEM%contains('basic')
debugCPFEM%extensive = debug_CPFEM%contains('extensive')
debugCPFEM%selective = debug_CPFEM%contains('selective')
debugCPFEM%element = config_debug%get_asInt('element',defaultVal = 1)
debugCPFEM%ip = config_debug%get_asInt('integrationpoint',defaultVal = 1)
2020-07-01 23:24:14 +05:30
if(debugCPFEM%basic) then
print'(a32,1x,6(i8,1x))', 'CPFEM_cs: ', shape(CPFEM_cs)
print'(a32,1x,6(i8,1x))', 'CPFEM_dcsdE: ', shape(CPFEM_dcsdE)
print'(a32,1x,6(i8,1x),/)', 'CPFEM_dcsdE_knownGood: ', shape(CPFEM_dcsdE_knownGood)
2020-09-22 16:39:12 +05:30
flush(IO_STDOUT)
endif
end subroutine CPFEM_init
2013-03-06 20:11:15 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief perform initialization at first call, update variables and call the actual material model
!--------------------------------------------------------------------------------------------------
2020-06-11 12:06:21 +05:30
subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyStress, jacobian)
integer, intent(in) :: elFE, & !< FE element number
2020-03-17 01:28:40 +05:30
ip !< integration point number
real(pReal), intent(in) :: dt !< time increment
real(pReal), dimension (3,3), intent(in) :: ffn, & !< deformation gradient for t=t0
ffn1 !< deformation gradient for t=t1
integer, intent(in) :: mode !< computation mode 1: regular computation plus aging of results
2020-03-17 01:28:40 +05:30
real(pReal), intent(in) :: temperature_inp !< temperature
real(pReal), dimension(6), intent(out) :: cauchyStress !< stress as 6 vector
real(pReal), dimension(6,6), intent(out) :: jacobian !< jacobian as 66 tensor (Consistent tangent dcs/dE)
2020-03-17 01:28:40 +05:30
real(pReal) J_inverse, & ! inverse of Jacobian
rnd
2020-06-16 10:11:53 +05:30
real(pReal), dimension (3,3) :: Kirchhoff ! Piola-Kirchhoff stress
2020-03-17 01:28:40 +05:30
real(pReal), dimension (3,3,3,3) :: H_sym, &
2020-06-11 12:06:21 +05:30
H
integer elCP, & ! crystal plasticity element number
2021-04-11 12:55:45 +05:30
i, j, k, l, m, n, ph, homog, mySource,ce
2020-06-16 10:11:53 +05:30
real(pReal), parameter :: ODD_STRESS = 1e15_pReal, & !< return value for stress if terminallyIll
ODD_JACOBIAN = 1e50_pReal !< return value for jacobian if terminallyIll
2020-06-17 18:51:51 +05:30
2021-04-11 00:22:46 +05:30
elCP = discretization_Marc_FEM2DAMASK_elem(elFE)
2021-04-11 12:55:45 +05:30
ce = discretization_Marc_FEM2DAMASK_cell(ip,elFE)
2020-07-03 20:15:11 +05:30
if (debugCPFEM%basic .and. elCP == debugCPFEM%element .and. ip == debugCPFEM%ip) then
print'(/,a)', '#############################################'
print'(a1,a22,1x,i8,a13)', '#','element', elCP, '#'
print'(a1,a22,1x,i8,a13)', '#','ip', ip, '#'
print'(a1,a22,1x,i8,a13)', '#','cycleCounter', cycleCounter, '#'
print'(a1,a22,1x,i8,a13)', '#','computationMode',mode, '#'
2020-03-17 01:28:40 +05:30
if (terminallyIll) &
print'(a,/)', '# --- terminallyIll --- #'
print'(a,/)', '#############################################'; flush (6)
2020-03-17 01:28:40 +05:30
endif
if (iand(mode, CPFEM_BACKUPJACOBIAN) /= 0) &
2020-03-17 01:28:40 +05:30
CPFEM_dcsde_knownGood = CPFEM_dcsde
if (iand(mode, CPFEM_RESTOREJACOBIAN) /= 0) &
2020-03-17 01:28:40 +05:30
CPFEM_dcsde = CPFEM_dcsde_knownGood
if (iand(mode, CPFEM_AGERESULTS) /= 0) call CPFEM_forward
2021-04-11 12:55:45 +05:30
homogenization_F0(1:3,1:3,ce) = ffn
homogenization_F(1:3,1:3,ce) = ffn1
if (iand(mode, CPFEM_CALCRESULTS) /= 0) then
2020-06-11 12:21:17 +05:30
validCalculation: if (terminallyIll) then
2020-03-17 01:28:40 +05:30
call random_number(rnd)
if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal
CPFEM_cs(1:6,ip,elCP) = ODD_STRESS * rnd
2020-09-13 14:48:09 +05:30
CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_eye(6)
2020-03-17 01:28:40 +05:30
else validCalculation
if (debugCPFEM%extensive) print'(a,i8,1x,i2)', '<< CPFEM >> calculation for elFE ip ',elFE,ip
2021-07-17 15:20:21 +05:30
call homogenization_mechanical_response(dt,[ip,ip],[elCP,elCP])
2021-07-17 00:13:08 +05:30
if (.not. terminallyIll) &
2021-07-17 15:20:21 +05:30
call homogenization_mechanical_response2(dt,[ip,ip],[elCP,elCP])
2020-03-17 01:28:40 +05:30
terminalIllness: if (terminallyIll) then
2020-03-17 01:28:40 +05:30
call random_number(rnd)
if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal
CPFEM_cs(1:6,ip,elCP) = ODD_STRESS * rnd
2020-09-13 14:48:09 +05:30
CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_eye(6)
2020-03-17 01:28:40 +05:30
else terminalIllness
2020-06-16 10:11:53 +05:30
! translate from P to sigma
2021-04-11 12:55:45 +05:30
Kirchhoff = matmul(homogenization_P(1:3,1:3,ce), transpose(homogenization_F(1:3,1:3,ce)))
J_inverse = 1.0_pReal / math_det33(homogenization_F(1:3,1:3,ce))
2020-03-17 01:28:40 +05:30
CPFEM_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.)
2020-03-17 01:28:40 +05:30
! translate from dP/dF to dCS/dE
H = 0.0_pReal
do i=1,3; do j=1,3; do k=1,3; do l=1,3; do m=1,3; do n=1,3
H(i,j,k,l) = H(i,j,k,l) &
2021-04-11 12:55:45 +05:30
+ homogenization_F(j,m,ce) * homogenization_F(l,n,ce) &
* homogenization_dPdF(i,m,k,n,ce) &
- math_delta(j,l) * homogenization_F(i,m,ce) * homogenization_P(k,m,ce) &
2020-03-17 01:28:40 +05:30
+ 0.5_pReal * ( Kirchhoff(j,l)*math_delta(i,k) + Kirchhoff(i,k)*math_delta(j,l) &
+ Kirchhoff(j,k)*math_delta(i,l) + Kirchhoff(i,l)*math_delta(j,k))
enddo; enddo; enddo; enddo; enddo; enddo
2020-03-17 01:28:40 +05:30
forall(i=1:3, j=1:3,k=1:3,l=1:3) &
H_sym(i,j,k,l) = 0.25_pReal * (H(i,j,k,l) + H(j,i,k,l) + H(i,j,l,k) + H(j,i,l,k))
2020-03-17 01:28:40 +05:30
CPFEM_dcsde(1:6,1:6,ip,elCP) = math_sym3333to66(J_inverse * H_sym,weighted=.false.)
2020-03-17 01:28:40 +05:30
endif terminalIllness
endif validCalculation
if (debugCPFEM%extensive &
2020-07-15 04:24:16 +05:30
.and. ((debugCPFEM%element == elCP .and. debugCPFEM%ip == ip) .or. .not. debugCPFEM%selective)) then
print'(a,i8,1x,i2,/,12x,6(f10.3,1x)/)', &
2020-03-17 01:28:40 +05:30
'<< CPFEM >> stress/MPa at elFE ip ', elFE, ip, CPFEM_cs(1:6,ip,elCP)*1.0e-6_pReal
print'(a,i8,1x,i2,/,6(12x,6(f10.3,1x)/))', &
2020-03-17 01:28:40 +05:30
'<< CPFEM >> Jacobian/GPa at elFE ip ', elFE, ip, transpose(CPFEM_dcsdE(1:6,1:6,ip,elCP))*1.0e-9_pReal
2020-09-22 16:39:12 +05:30
flush(IO_STDOUT)
2020-03-17 01:28:40 +05:30
endif
2020-03-17 01:28:40 +05:30
endif
2020-03-17 01:28:40 +05:30
if (all(abs(CPFEM_dcsdE(1:6,1:6,ip,elCP)) < 1e-10_pReal)) call IO_warning(601,elCP,ip)
2020-03-17 01:28:40 +05:30
cauchyStress = CPFEM_cs (1:6, ip,elCP)
jacobian = CPFEM_dcsdE(1:6,1:6,ip,elCP)
end subroutine CPFEM_general
!--------------------------------------------------------------------------------------------------
2020-02-25 22:23:15 +05:30
!> @brief Forward data for new time increment.
!--------------------------------------------------------------------------------------------------
subroutine CPFEM_forward
2020-12-23 11:28:54 +05:30
call homogenization_forward
2021-02-09 03:51:53 +05:30
call phase_forward
2020-02-25 22:23:15 +05:30
end subroutine CPFEM_forward
!--------------------------------------------------------------------------------------------------
2020-02-25 22:23:15 +05:30
!> @brief Trigger writing of results.
!--------------------------------------------------------------------------------------------------
subroutine CPFEM_results(inc,time)
integer, intent(in) :: inc
real(pReal), intent(in) :: time
call results_openJobFile
call results_addIncrement(inc,time)
2021-02-09 03:51:53 +05:30
call phase_results
2019-10-12 19:20:10 +05:30
call homogenization_results
call discretization_results
2020-01-10 06:15:00 +05:30
call results_finalizeIncrement
call results_closeJobFile
end subroutine CPFEM_results
end module CPFEM