From ada2beb8b8283da419e2ea409c987322405df754 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 1 Mar 2013 11:48:29 +0000 Subject: [PATCH] reorganized calculation modes for CPFEM, now having better readable and cleaner structure --- code/CPFEM.f90 | 1107 ++++++++++++++-------------- code/DAMASK_abaqus_exp.f | 17 +- code/DAMASK_abaqus_std.f | 62 +- code/DAMASK_marc.f90 | 297 ++++---- code/DAMASK_spectral_utilities.f90 | 21 +- code/IO.f90 | 4 +- code/crystallite.f90 | 4 +- 7 files changed, 749 insertions(+), 763 deletions(-) diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index 4286791ba..75914b395 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -19,26 +19,32 @@ !############################################################## !* $Id$ !############################################################## -MODULE CPFEM +module CPFEM !############################################################## ! *** CPFEM engine *** -! -use prec, only: pReal -implicit none - -real(pReal), parameter :: CPFEM_odd_stress = 1e15_pReal, & - CPFEM_odd_jacobian = 1e50_pReal -real(pReal), dimension (:,:,:), allocatable :: CPFEM_cs !> Cauchy stress -real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_dcsdE !> Cauchy stress tangent -real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_dcsdE_knownGood !> known good tangent + use prec, only: pReal, pInt + implicit none + real(pReal), parameter :: CPFEM_odd_stress = 1e15_pReal, & + CPFEM_odd_jacobian = 1e50_pReal -logical :: CPFEM_init_done = .false., & !> remember whether init has been done already - CPFEM_init_inProgress = .false., & !> remember whether first IP is currently performing init - CPFEM_calc_done = .false. !> remember whether first IP has already calced the results - + real(pReal), dimension (:,:,:), allocatable :: CPFEM_cs !> Cauchy stress + real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_dcsdE !> Cauchy stress tangent + real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_dcsdE_knownGood !> known good tangent -CONTAINS + logical :: CPFEM_init_done = .false., & !> remember whether init has been done already + CPFEM_init_inProgress = .false., & !> remember whether first IP is currently performing init + CPFEM_calc_done = .false. !> remember whether first IP has already calced the results + logical, public, protected :: usePingPong = .false. + integer(pInt), parameter, public :: & + CPFEM_CALCRESULTS = 2_pInt**0_pInt, & + CPFEM_AGERESULTS = 2_pInt**1_pInt, & + CPFEM_BACKUPJACOBIAN = 2_pInt**2_pInt, & + CPFEM_RESTOREJACOBIAN = 2_pInt**3_pInt, & + CPFEM_COLLECT = 2_pInt**4_pInt, & + CPFEM_EXPLICIT = 2_pInt**5_pInt + +contains !********************************************************* !*** call (thread safe) all module initializations *** @@ -46,64 +52,63 @@ CONTAINS subroutine CPFEM_initAll(Temperature,element,IP) - use prec, only: prec_init, & - pInt - use numerics, only: numerics_init - use debug, only: debug_init - use FEsolving, only: FE_init - use math, only: math_init - use mesh, only: mesh_init - use lattice, only: lattice_init - use material, only: material_init - use constitutive, only: constitutive_init - use crystallite, only: crystallite_init - use homogenization, only: homogenization_init - use IO, only: IO_init - use DAMASK_interface + use prec, only: prec_init + use numerics, only: numerics_init + use debug, only: debug_init + use FEsolving, only: FE_init + use math, only: math_init + use mesh, only: mesh_init + use lattice, only: lattice_init + use material, only: material_init + use constitutive, only: constitutive_init + use crystallite, only: crystallite_init + use homogenization, only: homogenization_init + use IO, only: IO_init + use DAMASK_interface - implicit none - integer(pInt), intent(in) :: element, & ! FE element number - IP ! FE integration point number - real(pReal), intent(in) :: Temperature ! temperature - real(pReal) rnd - integer(pInt) i,n + implicit none + integer(pInt), intent(in) :: element, & ! FE element number + IP ! FE integration point number + real(pReal), intent(in) :: Temperature ! temperature + real(pReal) rnd + integer(pInt) i,n - ! initialization step (three dimensional stress state check missing?) + ! initialization step (three dimensional stress state check missing?) - if (.not. CPFEM_init_done) then - call random_number(rnd) - do i=1,int(256.0*rnd) - n = n+1_pInt ! wasting random amount of time... - enddo ! ...to break potential race in multithreading - n = n+1_pInt - if (.not. CPFEM_init_inProgress) then ! yes my thread won! - CPFEM_init_inProgress = .true. + if (.not. CPFEM_init_done) then + call random_number(rnd) + do i=1,int(256.0*rnd) + n = n+1_pInt ! wasting random amount of time... + enddo ! ...to break potential race in multithreading + n = n+1_pInt + if (.not. CPFEM_init_inProgress) then ! yes my thread won! + CPFEM_init_inProgress = .true. #ifdef Spectral - call DAMASK_interface_init() ! Spectral solver is interfacing to commandline + call DAMASK_interface_init() ! Spectral solver is interfacing to commandline #endif - call prec_init - call IO_init - call numerics_init - call debug_init - call math_init - call FE_init - call mesh_init(IP, element) ! pass on coordinates to alter calcMode of first ip - call lattice_init - call material_init - call constitutive_init - call crystallite_init(Temperature) ! (have to) use temperature of first IP for whole model - call homogenization_init(Temperature) - call CPFEM_init + call prec_init + call IO_init + call numerics_init + call debug_init + call math_init + call FE_init + call mesh_init(IP, element) ! pass on coordinates to alter calcMode of first ip + call lattice_init + call material_init + call constitutive_init + call crystallite_init(Temperature) ! (have to) use temperature of first IP for whole model + call homogenization_init(Temperature) + call CPFEM_init #ifndef Spectral - call DAMASK_interface_init() ! Spectral solver init is already done + call DAMASK_interface_init() ! Spectral solver init is already done #endif - CPFEM_init_done = .true. - CPFEM_init_inProgress = .false. - else ! loser, loser... - do while (CPFEM_init_inProgress) - enddo - endif - endif + CPFEM_init_done = .true. + CPFEM_init_inProgress = .false. + else ! loser, loser... + do while (CPFEM_init_inProgress) + enddo + endif + endif end subroutine CPFEM_initAll @@ -113,114 +118,130 @@ end subroutine CPFEM_initAll !********************************************************* subroutine CPFEM_init - - use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) - use prec, only: pInt - use debug, only: debug_level, & - debug_CPFEM, & - debug_levelBasic, & - debug_levelExtensive - use IO, only: IO_read_jobBinaryFile,& - IO_read_jobBinaryIntFile, & - IO_timeStamp - use FEsolving, only: parallelExecution, & - symmetricSolver, & - restartRead, & - modelName - use mesh, only: mesh_NcpElems, & - mesh_maxNips - use material, only: homogenization_maxNgrains, & - material_phase - use constitutive, only: constitutive_state0 - use crystallite, only: crystallite_F0, & - crystallite_Fp0, & - crystallite_Lp0, & - crystallite_dPdF0, & - crystallite_Tstar0_v - use homogenization, only: homogenization_sizeState, & - homogenization_state0 + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) + use prec, only: & + pInt + use IO, only: & + IO_read_jobBinaryFile,& + IO_read_jobBinaryIntFile, & + IO_timeStamp, & + IO_error + use numerics, only: & + DAMASK_NumThreadsInt + use debug, only: & + debug_level, & + debug_CPFEM, & + debug_levelBasic, & + debug_levelExtensive + use FEsolving, only: & + parallelExecution, & + symmetricSolver, & + restartRead, & + modelName + use mesh, only: & + mesh_NcpElems, & + mesh_Nelems, & + mesh_maxNips + use material, only: & + homogenization_maxNgrains, & + material_phase + use constitutive, only: & + constitutive_state0 + use crystallite, only: & + crystallite_F0, & + crystallite_Fp0, & + crystallite_Lp0, & + crystallite_dPdF0, & + crystallite_Tstar0_v, & + crystallite_localPlasticity + use homogenization, only: & + homogenization_sizeState, & + homogenization_state0 + implicit none + integer(pInt) i,j,k,l,m - implicit none - integer(pInt) i,j,k,l,m + if (any(.not. crystallite_localPlasticity) .and. (mesh_Nelems /= mesh_NcpElems)) call IO_error(600) + if ((DAMASK_NumThreadsInt > 1_pInt) .and. (mesh_Nelems /= mesh_NcpElems)) call IO_error(601) + usePingPong = (any(.not. crystallite_localPlasticity) .or. (DAMASK_NumThreadsInt > 1_pInt)) - ! initialize stress and jacobian to zero - allocate(CPFEM_cs(6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_cs = 0.0_pReal - allocate(CPFEM_dcsdE(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dcsdE = 0.0_pReal - allocate(CPFEM_dcsdE_knownGood(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dcsdE_knownGood = 0.0_pReal + ! initialize stress and jacobian to zero + allocate(CPFEM_cs(6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_cs = 0.0_pReal + allocate(CPFEM_dcsdE(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dcsdE = 0.0_pReal + allocate(CPFEM_dcsdE_knownGood(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dcsdE_knownGood = 0.0_pReal - ! *** restore the last converged values of each essential variable from the binary file - if (restartRead) then - if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) then - !$OMP CRITICAL (write2out) - write(6,'(a)') '<< CPFEM >> restored state variables of last converged step from binary files' - !$OMP END CRITICAL (write2out) - endif + ! *** restore the last converged values of each essential variable from the binary file + if (restartRead) then + if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) then + !$OMP CRITICAL (write2out) + write(6,'(a)') '<< CPFEM >> restored state variables of last converged step from binary files' + !$OMP END CRITICAL (write2out) + endif - call IO_read_jobBinaryIntFile(777,'recordedPhase',modelName,size(material_phase)) - read (777,rec=1) material_phase - close (777) + call IO_read_jobBinaryIntFile(777,'recordedPhase',modelName,size(material_phase)) + read (777,rec=1) material_phase + close (777) - call IO_read_jobBinaryFile(777,'convergedF',modelName,size(crystallite_F0)) - read (777,rec=1) crystallite_F0 - close (777) + call IO_read_jobBinaryFile(777,'convergedF',modelName,size(crystallite_F0)) + read (777,rec=1) crystallite_F0 + close (777) - call IO_read_jobBinaryFile(777,'convergedFp',modelName,size(crystallite_Fp0)) - read (777,rec=1) crystallite_Fp0 - close (777) + call IO_read_jobBinaryFile(777,'convergedFp',modelName,size(crystallite_Fp0)) + read (777,rec=1) crystallite_Fp0 + close (777) - call IO_read_jobBinaryFile(777,'convergedLp',modelName,size(crystallite_Lp0)) - read (777,rec=1) crystallite_Lp0 - close (777) + call IO_read_jobBinaryFile(777,'convergedLp',modelName,size(crystallite_Lp0)) + read (777,rec=1) crystallite_Lp0 + close (777) - call IO_read_jobBinaryFile(777,'convergeddPdF',modelName,size(crystallite_dPdF0)) - read (777,rec=1) crystallite_dPdF0 - close (777) + call IO_read_jobBinaryFile(777,'convergeddPdF',modelName,size(crystallite_dPdF0)) + read (777,rec=1) crystallite_dPdF0 + close (777) - call IO_read_jobBinaryFile(777,'convergedTstar',modelName,size(crystallite_Tstar0_v)) - read (777,rec=1) crystallite_Tstar0_v - close (777) + call IO_read_jobBinaryFile(777,'convergedTstar',modelName,size(crystallite_Tstar0_v)) + read (777,rec=1) crystallite_Tstar0_v + close (777) - call IO_read_jobBinaryFile(777,'convergedStateConst',modelName) - m = 0_pInt - do i = 1,homogenization_maxNgrains; do j = 1,mesh_maxNips; do k = 1,mesh_NcpElems - do l = 1,size(constitutive_state0(i,j,k)%p) - m = m+1_pInt - read(777,rec=m) constitutive_state0(i,j,k)%p(l) - enddo - enddo; enddo; enddo - close (777) + call IO_read_jobBinaryFile(777,'convergedStateConst',modelName) + m = 0_pInt + do i = 1,homogenization_maxNgrains; do j = 1,mesh_maxNips; do k = 1,mesh_NcpElems + do l = 1,size(constitutive_state0(i,j,k)%p) + m = m+1_pInt + read(777,rec=m) constitutive_state0(i,j,k)%p(l) + enddo + enddo; enddo; enddo + close (777) - call IO_read_jobBinaryFile(777,'convergedStateHomog',modelName) - m = 0_pInt - do k = 1,mesh_NcpElems; do j = 1,mesh_maxNips - do l = 1,homogenization_sizeState(j,k) - m = m+1_pInt - read(777,rec=m) homogenization_state0(j,k)%p(l) - enddo - enddo; enddo - close (777) + call IO_read_jobBinaryFile(777,'convergedStateHomog',modelName) + m = 0_pInt + do k = 1,mesh_NcpElems; do j = 1,mesh_maxNips + do l = 1,homogenization_sizeState(j,k) + m = m+1_pInt + read(777,rec=m) homogenization_state0(j,k)%p(l) + enddo + enddo; enddo + close (777) - call IO_read_jobBinaryFile(777,'convergeddcsdE',modelName,size(CPFEM_dcsdE)) - read (777,rec=1) CPFEM_dcsdE - close (777) - restartRead = .false. - endif - ! *** end of restoring + call IO_read_jobBinaryFile(777,'convergeddcsdE',modelName,size(CPFEM_dcsdE)) + read (777,rec=1) CPFEM_dcsdE + close (777) + restartRead = .false. + endif + ! *** end of restoring - write(6,'(/,a)') '<<<+- CPFEM init -+>>>' - write(6,'(a)') '$Id$' - write(6,'(a16,a)') ' Current time : ',IO_timeStamp() + write(6,'(/,a)') '<<<+- CPFEM init -+>>>' + write(6,'(a)') '$Id$' + write(6,'(a16,a)') ' Current time : ',IO_timeStamp() #include "compilation_info.f90" - if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0) then - write(6,'(a32,1x,6(i8,1x))') 'CPFEM_cs: ', shape(CPFEM_cs) - write(6,'(a32,1x,6(i8,1x))') 'CPFEM_dcsdE: ', shape(CPFEM_dcsdE) - write(6,'(a32,1x,6(i8,1x))') 'CPFEM_dcsdE_knownGood: ', shape(CPFEM_dcsdE_knownGood) - write(6,*) - write(6,*) 'parallelExecution: ', parallelExecution - write(6,*) 'symmetricSolver: ', symmetricSolver - endif - flush(6) + if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0) then + write(6,'(a32,1x,6(i8,1x))') 'CPFEM_cs: ', shape(CPFEM_cs) + write(6,'(a32,1x,6(i8,1x))') 'CPFEM_dcsdE: ', shape(CPFEM_dcsdE) + write(6,'(a32,1x,6(i8,1x))') 'CPFEM_dcsdE_knownGood: ', shape(CPFEM_dcsdE_knownGood) + write(6,*) + write(6,*) 'parallelExecution: ', parallelExecution + write(6,*) 'symmetricSolver: ', symmetricSolver + endif + flush(6) end subroutine CPFEM_init @@ -230,425 +251,381 @@ end subroutine CPFEM_init !*** call the actual material model *** !*********************************************************************** subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchyStress,& - & jacobian, pstress, dPdF) - ! note: cauchyStress = Cauchy stress cs(6) and jacobian = Consistent tangent dcs/dE + & jacobian, pstress, dPdF) + ! note: cauchyStress = Cauchy stress cs(6) and jacobian = Consistent tangent dcs/dE + use numerics, only: defgradTolerance, & + iJacoStiffness + use debug, only: debug_level, & + debug_CPFEM, & + debug_levelBasic, & + debug_levelExtensive, & + debug_levelSelective, & + debug_e, & + debug_i, & + debug_stressMaxLocation, & + debug_stressMinLocation, & + debug_jacobianMaxLocation, & + debug_jacobianMinLocation, & + debug_stressMax, & + debug_stressMin, & + debug_jacobianMax, & + debug_jacobianMin + use FEsolving, only: parallelExecution, & + outdatedFFN1, & + terminallyIll, & + cycleCounter, & + theInc, & + theTime, & + theDelta, & + FEsolving_execElem, & + FEsolving_execIP, & + restartWrite + use math, only: math_identity2nd, & + math_mul33x33, & + math_det33, & + math_transpose33, & + math_I3, & + math_Mandel3333to66, & + math_Mandel66to3333, & + math_Mandel33to6, & + math_Mandel6to33 + use mesh, only: mesh_FEasCP, & + mesh_NcpElems, & + mesh_maxNips, & + mesh_element, & + FE_Nips, & + FE_Nnodes, & + FE_geomtype + use material, only: homogenization_maxNgrains, & + microstructure_elemhomo, & + material_phase + use constitutive, only: constitutive_state0,constitutive_state + use crystallite, only: crystallite_partionedF,& + crystallite_F0, & + crystallite_Fp0, & + crystallite_Fp, & + crystallite_Lp0, & + crystallite_Lp, & + crystallite_dPdF0, & + crystallite_dPdF, & + crystallite_Tstar0_v, & + crystallite_Tstar_v + use homogenization, only: homogenization_sizeState, & + homogenization_state, & + homogenization_state0, & + materialpoint_F, & + materialpoint_F0, & + materialpoint_P, & + materialpoint_dPdF, & + materialpoint_results, & + materialpoint_sizeResults, & + materialpoint_Temperature, & + materialpoint_stressAndItsTangent, & + materialpoint_postResults + use IO, only: IO_write_jobBinaryFile, & + IO_warning + use DAMASK_interface + + implicit none + integer(pInt), intent(in) :: element, & ! FE element number + IP ! FE integration point number + real(pReal), intent(inout) :: Temperature ! temperature + 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(pInt), intent(in) :: mode ! computation mode 1: regular computation plus aging of results + + real(pReal), dimension(6), intent(out) :: cauchyStress ! stress vector in Mandel notation + real(pReal), dimension(6,6), intent(out) :: jacobian ! jacobian in Mandel notation + real(pReal), dimension (3,3), intent(out) :: pstress ! Piola-Kirchhoff stress in Matrix notation + real(pReal), dimension (3,3,3,3), intent(out) :: dPdF ! + + real(pReal) J_inverse, & ! inverse of Jacobian + rnd + real(pReal), dimension (3,3) :: Kirchhoff, & ! Piola-Kirchhoff stress in Matrix notation + cauchyStress33 ! stress vector in Matrix notation + real(pReal), dimension (3,3,3,3) :: H_sym, & + H, & + jacobian3333 ! jacobian in Matrix notation + integer(pInt) cp_en, & ! crystal plasticity element number + i, j, k, l, m, n, e + logical updateJaco ! flag indicating if JAcobian has to be updated + + + cp_en = mesh_FEasCP('elem',element) + + if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt .and. cp_en == 1 .and. IP == 1) then + !$OMP CRITICAL (write2out) + write(6,'(/,a)') '#############################################' + write(6,'(a1,a22,1x,f15.7,a6)') '#','theTime', theTime, '#' + write(6,'(a1,a22,1x,f15.7,a6)') '#','theDelta', theDelta, '#' + write(6,'(a1,a22,1x,i8,a13)') '#','theInc', theInc, '#' + write(6,'(a1,a22,1x,i8,a13)') '#','cycleCounter', cycleCounter, '#' + write(6,'(a1,a22,1x,i8,a13)') '#','computationMode',mode, '#' + write(6,'(a,/)') '#############################################' + flush (6) + !$OMP END CRITICAL (write2out) + endif - !*** variables and functions from other modules ***! - use prec, only: pInt - use numerics, only: defgradTolerance, & - iJacoStiffness - use debug, only: debug_level, & - debug_CPFEM, & - debug_levelBasic, & - debug_levelExtensive, & - debug_levelSelective, & - debug_e, & - debug_i, & - debug_stressMaxLocation, & - debug_stressMinLocation, & - debug_jacobianMaxLocation, & - debug_jacobianMinLocation, & - debug_stressMax, & - debug_stressMin, & - debug_jacobianMax, & - debug_jacobianMin - use FEsolving, only: parallelExecution, & - outdatedFFN1, & - terminallyIll, & - cycleCounter, & - theInc, & - theTime, & - theDelta, & - FEsolving_execElem, & - FEsolving_execIP, & - restartWrite - use math, only: math_identity2nd, & - math_mul33x33, & - math_det33, & - math_transpose33, & - math_I3, & - math_Mandel3333to66, & - math_Mandel66to3333, & - math_Mandel33to6, & - math_Mandel6to33 - use mesh, only: mesh_FEasCP, & - mesh_NcpElems, & - mesh_maxNips, & - mesh_element, & - FE_Nips, & - FE_Nnodes, & - FE_geomtype - use material, only: homogenization_maxNgrains, & - microstructure_elemhomo, & - material_phase - use constitutive, only: constitutive_state0,constitutive_state - use crystallite, only: crystallite_partionedF,& - crystallite_F0, & - crystallite_Fp0, & - crystallite_Fp, & - crystallite_Lp0, & - crystallite_Lp, & - crystallite_dPdF0, & - crystallite_dPdF, & - crystallite_Tstar0_v, & - crystallite_Tstar_v - use homogenization, only: homogenization_sizeState, & - homogenization_state, & - homogenization_state0, & - materialpoint_F, & - materialpoint_F0, & - materialpoint_P, & - materialpoint_dPdF, & - materialpoint_results, & - materialpoint_sizeResults, & - materialpoint_Temperature, & - materialpoint_stressAndItsTangent, & - materialpoint_postResults - use IO, only: IO_write_jobBinaryFile, & - IO_warning - use DAMASK_interface - - implicit none - - - !*** input variables ***! - - integer(pInt), intent(in) :: element, & ! FE element number - IP ! FE integration point number - real(pReal), intent(inout) :: Temperature ! temperature - 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(pInt), intent(in) :: mode ! computation mode 1: regular computation plus aging of results - ! 2: regular computation - ! 3: collection of FEM data - ! 4: backup tangent from former converged inc - ! 5: restore tangent from former converged inc - ! 6: recycling of former results (MARC speciality) - - - !*** output variables ***! - - real(pReal), dimension(6), intent(out) :: cauchyStress ! stress vector in Mandel notation - real(pReal), dimension(6,6), intent(out) :: jacobian ! jacobian in Mandel notation - real(pReal), dimension (3,3), intent(out) :: pstress ! Piola-Kirchhoff stress in Matrix notation - real(pReal), dimension (3,3,3,3), intent(out) :: dPdF ! - - - !*** local variables ***! - - real(pReal) J_inverse, & ! inverse of Jacobian - rnd - real(pReal), dimension (3,3) :: Kirchhoff, & ! Piola-Kirchhoff stress in Matrix notation - cauchyStress33 ! stress vector in Matrix notation - real(pReal), dimension (3,3,3,3) :: H_sym, & - H, & - jacobian3333 ! jacobian in Matrix notation - integer(pInt) cp_en, & ! crystal plasticity element number - i, j, k, l, m, n, e - logical updateJaco ! flag indicating if JAcobian has to be updated - - - cp_en = mesh_FEasCP('elem',element) - - if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt .and. cp_en == 1 .and. IP == 1) then - !$OMP CRITICAL (write2out) - write(6,'(/,a)') '#############################################' - write(6,'(a1,a22,1x,f15.7,a6)') '#','theTime', theTime, '#' - write(6,'(a1,a22,1x,f15.7,a6)') '#','theDelta', theDelta, '#' - write(6,'(a1,a22,1x,i8,a13)') '#','theInc', theInc, '#' - write(6,'(a1,a22,1x,i8,a13)') '#','cycleCounter', cycleCounter, '#' - write(6,'(a1,a22,1x,i8,a13)') '#','computationMode',mode, '#' - write(6,'(a,/)') '#############################################' - flush (6) - !$OMP END CRITICAL (write2out) - endif + if (iand(mode, CPFEM_AGERESULTS) /= 0_pInt) then + crystallite_F0 = crystallite_partionedF ! crystallite deformation (_subF is perturbed...) + crystallite_Fp0 = crystallite_Fp ! crystallite plastic deformation + crystallite_Lp0 = crystallite_Lp ! crystallite plastic velocity + crystallite_dPdF0 = crystallite_dPdF ! crystallite stiffness + crystallite_Tstar0_v = crystallite_Tstar_v ! crystallite 2nd Piola Kirchhoff stress + forall ( i = 1:homogenization_maxNgrains, & + j = 1:mesh_maxNips, & + k = 1:mesh_NcpElems ) & + constitutive_state0(i,j,k)%p = constitutive_state(i,j,k)%p ! microstructure of crystallites + if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) then + !$OMP CRITICAL (write2out) + write(6,'(a)') '<< CPFEM >> aging states' + if (debug_e <= mesh_NcpElems .and. debug_i <= mesh_maxNips) then + write(6,'(a,1x,i8,1x,i2,1x,i4,/,(12x,6(e20.8,1x)))') '<< CPFEM >> aged state of el ip grain',& + debug_e, debug_i, 1, constitutive_state(1,debug_i,debug_e)%p + write(6,*) + endif + !$OMP END CRITICAL (write2out) + endif + !$OMP PARALLEL DO + do k = 1,mesh_NcpElems + do j = 1,mesh_maxNips + if (homogenization_sizeState(j,k) > 0_pInt) & + homogenization_state0(j,k)%p = homogenization_state(j,k)%p ! internal state of homogenization scheme + enddo + enddo + !$OMP END PARALLEL DO - - !*** according to our "mode" we decide what to do - - select case (mode) - - - ! --+>> REGULAR COMPUTATION (WITH AGING OF RESULTS IF MODE == 1) <<+-- - - case (1,2,8,9) - - !*** age results - - if (mode == 1 .or. mode == 8) then - crystallite_F0 = crystallite_partionedF ! crystallite deformation (_subF is perturbed...) - crystallite_Fp0 = crystallite_Fp ! crystallite plastic deformation - crystallite_Lp0 = crystallite_Lp ! crystallite plastic velocity - crystallite_dPdF0 = crystallite_dPdF ! crystallite stiffness - crystallite_Tstar0_v = crystallite_Tstar_v ! crystallite 2nd Piola Kirchhoff stress - forall ( i = 1:homogenization_maxNgrains, & - j = 1:mesh_maxNips, & - k = 1:mesh_NcpElems ) & - constitutive_state0(i,j,k)%p = constitutive_state(i,j,k)%p ! microstructure of crystallites - if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) then + ! * dump the last converged values of each essential variable to a binary file + + if (restartWrite) then + if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) then !$OMP CRITICAL (write2out) - write(6,'(a)') '<< CPFEM >> aging states' - if (debug_e <= mesh_NcpElems .and. debug_i <= mesh_maxNips) then - write(6,'(a,1x,i8,1x,i2,1x,i4,/,(12x,6(e20.8,1x)))') '<< CPFEM >> aged state of el ip grain',& - debug_e, debug_i, 1, constitutive_state(1,debug_i,debug_e)%p - write(6,*) - endif + write(6,'(a)') '<< CPFEM >> writing state variables of last converged step to binary files' !$OMP END CRITICAL (write2out) - endif - !$OMP PARALLEL DO - do k = 1,mesh_NcpElems - do j = 1,mesh_maxNips - if (homogenization_sizeState(j,k) > 0_pInt) & - homogenization_state0(j,k)%p = homogenization_state(j,k)%p ! internal state of homogenization scheme - enddo - enddo - !$OMP END PARALLEL DO + endif + + call IO_write_jobBinaryFile(777,'recordedPhase',size(material_phase)) + write (777,rec=1) material_phase + close (777) - ! * dump the last converged values of each essential variable to a binary file - - if (restartWrite) then - if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) then - !$OMP CRITICAL (write2out) - write(6,'(a)') '<< CPFEM >> writing state variables of last converged step to binary files' - !$OMP END CRITICAL (write2out) - endif - - call IO_write_jobBinaryFile(777,'recordedPhase',size(material_phase)) - write (777,rec=1) material_phase - close (777) + call IO_write_jobBinaryFile(777,'convergedF',size(crystallite_F0)) + write (777,rec=1) crystallite_F0 + close (777) + + call IO_write_jobBinaryFile(777,'convergedFp',size(crystallite_Fp0)) + write (777,rec=1) crystallite_Fp0 + close (777) + + call IO_write_jobBinaryFile(777,'convergedLp',size(crystallite_Lp0)) + write (777,rec=1) crystallite_Lp0 + close (777) + + call IO_write_jobBinaryFile(777,'convergeddPdF',size(crystallite_dPdF0)) + write (777,rec=1) crystallite_dPdF0 + close (777) + + call IO_write_jobBinaryFile(777,'convergedTstar',size(crystallite_Tstar0_v)) + write (777,rec=1) crystallite_Tstar0_v + close (777) + + call IO_write_jobBinaryFile(777,'convergedStateConst') + m = 0_pInt + do i = 1,homogenization_maxNgrains; do j = 1,mesh_maxNips; do k = 1,mesh_NcpElems + do l = 1,size(constitutive_state0(i,j,k)%p) + m = m+1_pInt + write(777,rec=m) constitutive_state0(i,j,k)%p(l) + enddo + enddo; enddo; enddo + close (777) + + call IO_write_jobBinaryFile(777,'convergedStateHomog') + m = 0_pInt + do k = 1,mesh_NcpElems; do j = 1,mesh_maxNips + do l = 1,homogenization_sizeState(j,k) + m = m+1_pInt + write(777,rec=m) homogenization_state0(j,k)%p(l) + enddo + enddo; enddo + close (777) + + call IO_write_jobBinaryFile(777,'convergeddcsdE',size(CPFEM_dcsdE)) + write (777,rec=1) CPFEM_dcsdE + close (777) + + endif + endif - call IO_write_jobBinaryFile(777,'convergedF',size(crystallite_F0)) - write (777,rec=1) crystallite_F0 - close (777) - - call IO_write_jobBinaryFile(777,'convergedFp',size(crystallite_Fp0)) - write (777,rec=1) crystallite_Fp0 - close (777) - - call IO_write_jobBinaryFile(777,'convergedLp',size(crystallite_Lp0)) - write (777,rec=1) crystallite_Lp0 - close (777) - - call IO_write_jobBinaryFile(777,'convergeddPdF',size(crystallite_dPdF0)) - write (777,rec=1) crystallite_dPdF0 - close (777) - - call IO_write_jobBinaryFile(777,'convergedTstar',size(crystallite_Tstar0_v)) - write (777,rec=1) crystallite_Tstar0_v - close (777) - - call IO_write_jobBinaryFile(777,'convergedStateConst') - m = 0_pInt - do i = 1,homogenization_maxNgrains; do j = 1,mesh_maxNips; do k = 1,mesh_NcpElems - do l = 1,size(constitutive_state0(i,j,k)%p) - m = m+1_pInt - write(777,rec=m) constitutive_state0(i,j,k)%p(l) - enddo - enddo; enddo; enddo - close (777) - - call IO_write_jobBinaryFile(777,'convergedStateHomog') - m = 0_pInt - do k = 1,mesh_NcpElems; do j = 1,mesh_maxNips - do l = 1,homogenization_sizeState(j,k) - m = m+1_pInt - write(777,rec=m) homogenization_state0(j,k)%p(l) - enddo - enddo; enddo - close (777) - - call IO_write_jobBinaryFile(777,'convergeddcsdE',size(CPFEM_dcsdE)) - write (777,rec=1) CPFEM_dcsdE - close (777) - - endif - ! * end of dumping - endif + if (iand(mode, CPFEM_EXPLICIT) /= 0_pInt) then ! Abaqus explicit skips collect + materialpoint_Temperature(IP,cp_en) = Temperature + materialpoint_F0(1:3,1:3,IP,cp_en) = ffn + materialpoint_F(1:3,1:3,IP,cp_en) = ffn1 + endif - if (mode == 8 .or. mode == 9) then ! Abaqus explicit skips collect - materialpoint_Temperature(IP,cp_en) = Temperature - materialpoint_F0(1:3,1:3,IP,cp_en) = ffn - materialpoint_F(1:3,1:3,IP,cp_en) = ffn1 - endif - - - !*** deformation gradient outdated or any actual deformation gradient differs more than relevantStrain from the stored one - - if (terminallyIll .or. outdatedFFN1 .or. any(abs(ffn1 - materialpoint_F(1:3,1:3,IP,cp_en)) > defgradTolerance)) then + if (iand(mode, CPFEM_CALCRESULTS) /= 0_pInt) then + !*** deformation gradient outdated or any actual deformation gradient differs more than relevantStrain from the stored one + + if (terminallyIll .or. outdatedFFN1 .or. any(abs(ffn1 - materialpoint_F(1:3,1:3,IP,cp_en)) > defgradTolerance)) then ! if (.not. terminallyIll .and. .not. outdatedFFN1) then - if (any(abs(ffn1 - materialpoint_F(1:3,1:3,IP,cp_en)) > defgradTolerance)) then - if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) then - !$OMP CRITICAL (write2out) - write(6,'(a,1x,i8,1x,i2)') '<< CPFEM >> OUTDATED at el ip',cp_en,IP - write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 old:',math_transpose33(materialpoint_F(1:3,1:3,IP,cp_en)) - write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 now:',math_transpose33(ffn1) - !$OMP END CRITICAL (write2out) - endif - outdatedFFN1 = .true. - endif - call random_number(rnd) - if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal - CPFEM_cs(1:6,IP,cp_en) = rnd*CPFEM_odd_stress - CPFEM_dcsde(1:6,1:6,IP,cp_en) = CPFEM_odd_jacobian*math_identity2nd(6) - - - !*** deformation gradient is not outdated - - else - updateJaco = mod(cycleCounter,iJacoStiffness) == 0 - - !* no parallel computation, so we use just one single element and IP for computation - - if (.not. parallelExecution) then - FEsolving_execElem(1) = cp_en - FEsolving_execElem(2) = cp_en - if (.not. microstructure_elemhomo(mesh_element(4,cp_en)) .or. & ! calculate unless homogeneous - (microstructure_elemhomo(mesh_element(4,cp_en)) .and. IP == 1_pInt)) then ! and then only first IP - FEsolving_execIP(1,cp_en) = IP - FEsolving_execIP(2,cp_en) = IP - if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) then - !$OMP CRITICAL (write2out) - write(6,'(a,i8,1x,i2)') '<< CPFEM >> calculation for el ip ',cp_en,IP - !$OMP END CRITICAL (write2out) - endif - call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent - call materialpoint_postResults(dt) ! post results - endif - - !* parallel computation and calulation not yet done - - elseif (.not. CPFEM_calc_done) then - if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) then - !$OMP CRITICAL (write2out) - write(6,'(a,i8,a,i8)') '<< CPFEM >> calculation for elements ',FEsolving_execElem(1),' to ',FEsolving_execElem(2) - !$OMP END CRITICAL (write2out) - endif - call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent (parallel execution inside) - call materialpoint_postResults(dt) ! post results - CPFEM_calc_done = .true. - endif - - !*** map stress and stiffness (or return odd values if terminally ill) - - if ( terminallyIll ) then - call random_number(rnd) - if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal - CPFEM_cs(1:6,IP,cp_en) = rnd * CPFEM_odd_stress - CPFEM_dcsde(1:6,1:6,IP,cp_en) = CPFEM_odd_jacobian * math_identity2nd(6) - else - if (microstructure_elemhomo(mesh_element(4,cp_en)) .and. IP > 1_pInt) then ! me homogenous? --> copy from first IP - materialpoint_P(1:3,1:3,IP,cp_en) = materialpoint_P(1:3,1:3,1,cp_en) - materialpoint_F(1:3,1:3,IP,cp_en) = materialpoint_F(1:3,1:3,1,cp_en) - materialpoint_dPdF(1:3,1:3,1:3,1:3,IP,cp_en) = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,cp_en) - materialpoint_results(1:materialpoint_sizeResults,IP,cp_en) = materialpoint_results(1:materialpoint_sizeResults,1,cp_en) - endif - ! translate from P to CS - Kirchhoff = math_mul33x33(materialpoint_P(1:3,1:3,IP,cp_en), math_transpose33(materialpoint_F(1:3,1:3,IP,cp_en))) - J_inverse = 1.0_pReal / math_det33(materialpoint_F(1:3,1:3,IP,cp_en)) - CPFEM_cs(1:6,IP,cp_en) = math_Mandel33to6(J_inverse * Kirchhoff) + if (any(abs(ffn1 - materialpoint_F(1:3,1:3,IP,cp_en)) > defgradTolerance)) then + if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) then + !$OMP CRITICAL (write2out) + write(6,'(a,1x,i8,1x,i2)') '<< CPFEM >> OUTDATED at el ip',cp_en,IP + write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 old:',math_transpose33(materialpoint_F(1:3,1:3,IP,cp_en)) + write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 now:',math_transpose33(ffn1) + !$OMP END CRITICAL (write2out) + endif + outdatedFFN1 = .true. + endif + call random_number(rnd) + if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal + CPFEM_cs(1:6,IP,cp_en) = rnd*CPFEM_odd_stress + CPFEM_dcsde(1:6,1:6,IP,cp_en) = CPFEM_odd_jacobian*math_identity2nd(6) + + + !*** deformation gradient is not outdated + + else + updateJaco = mod(cycleCounter,iJacoStiffness) == 0 + + !* no parallel computation, so we use just one single element and IP for computation + + if (.not. parallelExecution) then + FEsolving_execElem(1) = cp_en + FEsolving_execElem(2) = cp_en + if (.not. microstructure_elemhomo(mesh_element(4,cp_en)) .or. & ! calculate unless homogeneous + (microstructure_elemhomo(mesh_element(4,cp_en)) .and. IP == 1_pInt)) then ! and then only first IP + FEsolving_execIP(1,cp_en) = IP + FEsolving_execIP(2,cp_en) = IP + if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) then + !$OMP CRITICAL (write2out) + write(6,'(a,i8,1x,i2)') '<< CPFEM >> calculation for el ip ',cp_en,IP + !$OMP END CRITICAL (write2out) + endif + call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent + call materialpoint_postResults(dt) ! post results + endif + + !* parallel computation and calulation not yet done + + elseif (.not. CPFEM_calc_done) then + if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) then + !$OMP CRITICAL (write2out) + write(6,'(a,i8,a,i8)') '<< CPFEM >> calculation for elements ',FEsolving_execElem(1),' to ',FEsolving_execElem(2) + !$OMP END CRITICAL (write2out) + endif + call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent (parallel execution inside) + call materialpoint_postResults(dt) ! post results + CPFEM_calc_done = .true. + endif + + !*** map stress and stiffness (or return odd values if terminally ill) + + if ( terminallyIll ) then + call random_number(rnd) + if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal + CPFEM_cs(1:6,IP,cp_en) = rnd * CPFEM_odd_stress + CPFEM_dcsde(1:6,1:6,IP,cp_en) = CPFEM_odd_jacobian * math_identity2nd(6) + else + if (microstructure_elemhomo(mesh_element(4,cp_en)) .and. IP > 1_pInt) then ! me homogenous? --> copy from first IP + materialpoint_P(1:3,1:3,IP,cp_en) = materialpoint_P(1:3,1:3,1,cp_en) + materialpoint_F(1:3,1:3,IP,cp_en) = materialpoint_F(1:3,1:3,1,cp_en) + materialpoint_dPdF(1:3,1:3,1:3,1:3,IP,cp_en) = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,cp_en) + materialpoint_results(1:materialpoint_sizeResults,IP,cp_en) = materialpoint_results(1:materialpoint_sizeResults,1,cp_en) + endif + ! translate from P to CS + Kirchhoff = math_mul33x33(materialpoint_P(1:3,1:3,IP,cp_en), math_transpose33(materialpoint_F(1:3,1:3,IP,cp_en))) + J_inverse = 1.0_pReal / math_det33(materialpoint_F(1:3,1:3,IP,cp_en)) + CPFEM_cs(1:6,IP,cp_en) = math_Mandel33to6(J_inverse * Kirchhoff) - ! 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) + & - materialpoint_F(j,m,IP,cp_en) * & - materialpoint_F(l,n,IP,cp_en) * & - materialpoint_dPdF(i,m,k,n,IP,cp_en) - & - math_I3(j,l) * materialpoint_F(i,m,IP,cp_en) * materialpoint_P(k,m,IP,cp_en) + & - 0.5_pReal * (math_I3(i,k) * Kirchhoff(j,l) + math_I3(j,l) * Kirchhoff(i,k) + & - math_I3(i,l) * Kirchhoff(j,k) + math_I3(j,k) * Kirchhoff(i,l)) - enddo; enddo; enddo; enddo; enddo; enddo - do i=1,3; do j=1,3; do k=1,3; do 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)) - enddo; enddo; enddo; enddo - CPFEM_dcsde(1:6,1:6,IP,cp_en) = math_Mandel3333to66(J_inverse * H_sym) - endif - endif - + ! 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) + & + materialpoint_F(j,m,IP,cp_en) * & + materialpoint_F(l,n,IP,cp_en) * & + materialpoint_dPdF(i,m,k,n,IP,cp_en) - & + math_I3(j,l) * materialpoint_F(i,m,IP,cp_en) * materialpoint_P(k,m,IP,cp_en) + & + 0.5_pReal * (math_I3(i,k) * Kirchhoff(j,l) + math_I3(j,l) * Kirchhoff(i,k) + & + math_I3(i,l) * Kirchhoff(j,k) + math_I3(j,k) * Kirchhoff(i,l)) + enddo; enddo; enddo; enddo; enddo; enddo + do i=1,3; do j=1,3; do k=1,3; do 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)) + enddo; enddo; enddo; enddo + CPFEM_dcsde(1:6,1:6,IP,cp_en) = math_Mandel3333to66(J_inverse * H_sym) + endif + endif + endif - ! --+>> COLLECTION OF FEM INPUT WITH RETURNING OF RANDOMIZED ODD STRESS AND JACOBIAN <<+-- - - case (3,4,5) - if (mode == 4) then - CPFEM_dcsde_knownGood = CPFEM_dcsde ! --+>> BACKUP JACOBIAN FROM FORMER CONVERGED INC - else if (mode == 5) then - CPFEM_dcsde = CPFEM_dcsde_knownGood ! --+>> RESTORE CONSISTENT JACOBIAN FROM FORMER CONVERGED INC - end if - call random_number(rnd) - if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal - materialpoint_Temperature(IP,cp_en) = Temperature - materialpoint_F0(1:3,1:3,IP,cp_en) = ffn - materialpoint_F(1:3,1:3,IP,cp_en) = ffn1 - CPFEM_cs(1:6,IP,cp_en) = rnd * CPFEM_odd_stress - CPFEM_dcsde(1:6,1:6,IP,cp_en) = CPFEM_odd_jacobian * math_identity2nd(6) - CPFEM_calc_done = .false. - - ! --+>> RECYCLING OF FORMER RESULTS (MARC SPECIALTY) <<+-- - - case (6) - ! do nothing + ! --+>> COLLECTION OF FEM INPUT WITH RETURNING OF RANDOMIZED ODD STRESS AND JACOBIAN <<+-- - ! --+>> RESTORE CONSISTENT JACOBIAN FROM FORMER CONVERGED INC + if (iand(mode, CPFEM_BACKUPJACOBIAN) /= 0_pInt) & + CPFEM_dcsde_knownGood = CPFEM_dcsde + if (iand(mode, CPFEM_RESTOREJACOBIAN) /= 0_pInt) & + CPFEM_dcsde = CPFEM_dcsde_knownGood + + if (iand(mode, CPFEM_COLLECT) /= 0_pInt) then + call random_number(rnd) + if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal + materialpoint_Temperature(IP,cp_en) = Temperature + materialpoint_F0(1:3,1:3,IP,cp_en) = ffn + materialpoint_F(1:3,1:3,IP,cp_en) = ffn1 + CPFEM_cs(1:6,IP,cp_en) = rnd * CPFEM_odd_stress + CPFEM_dcsde(1:6,1:6,IP,cp_en) = CPFEM_odd_jacobian * math_identity2nd(6) + CPFEM_calc_done = .false. + endif - case (7) - CPFEM_dcsde = CPFEM_dcsde_knownGood - - end select + + cauchyStress = CPFEM_cs(1:6,IP,cp_en) + jacobian = CPFEM_dcsdE(1:6,1:6,IP,cp_en) + pstress = materialpoint_P(1:3,1:3,IP,cp_en) + dPdF = materialpoint_dPdF(1:3,1:3,1:3,1:3,IP,cp_en) + if (theTime > 0.0_pReal) then + Temperature = materialpoint_Temperature(IP,cp_en) ! homogenized result except for potentially non-isothermal starting condition. + endif - - !*** fill output variables with computed values - - cauchyStress = CPFEM_cs(1:6,IP,cp_en) - jacobian = CPFEM_dcsdE(1:6,1:6,IP,cp_en) - pstress = materialpoint_P(1:3,1:3,IP,cp_en) - dPdF = materialpoint_dPdF(1:3,1:3,1:3,1:3,IP,cp_en) - if (theTime > 0.0_pReal) then - Temperature = materialpoint_Temperature(IP,cp_en) ! homogenized result except for potentially non-isothermal starting condition. - endif - - if (mode < 3 .and. iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt & - .and. ((debug_e == cp_en .and. debug_i == IP) & - .or. .not. iand(debug_level(debug_CPFEM), debug_levelSelective) /= 0_pInt)) then - !$OMP CRITICAL (write2out) - write(6,'(a,i8,1x,i2,/,12x,6(f10.3,1x)/)') '<< CPFEM >> stress/MPa at el ip ', cp_en, IP, cauchyStress/1.0e6_pReal - write(6,'(a,i8,1x,i2,/,6(12x,6(f10.3,1x)/))') '<< CPFEM >> Jacobian/GPa at el ip ', cp_en, IP, transpose(jacobian)/1.0e9_pReal - flush(6) - !$OMP END CRITICAL (write2out) - endif - - - !*** warn if stiffness close to zero - - if (all(abs(jacobian) < 1e-10_pReal)) then - call IO_warning(601,cp_en,IP) - endif - - - !*** remember extreme values of stress and jacobian - - if (mode < 3) then - cauchyStress33 = math_Mandel6to33(cauchyStress) - if (maxval(cauchyStress33) > debug_stressMax) then - debug_stressMaxLocation = (/cp_en, IP/) - debug_stressMax = maxval(cauchyStress33) - endif - if (minval(cauchyStress33) < debug_stressMin) then - debug_stressMinLocation = (/cp_en, IP/) - debug_stressMin = minval(cauchyStress33) - endif - jacobian3333 = math_Mandel66to3333(jacobian) - if (maxval(jacobian3333) > debug_jacobianMax) then - debug_jacobianMaxLocation = (/cp_en, IP/) - debug_jacobianMax = maxval(jacobian3333) - endif - if (minval(jacobian3333) < debug_jacobianMin) then - debug_jacobianMinLocation = (/cp_en, IP/) - debug_jacobianMin = minval(jacobian3333) - endif - endif - + if (mode < 3 .and. iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt & + .and. ((debug_e == cp_en .and. debug_i == IP) & + .or. .not. iand(debug_level(debug_CPFEM), debug_levelSelective) /= 0_pInt)) then + !$OMP CRITICAL (write2out) + write(6,'(a,i8,1x,i2,/,12x,6(f10.3,1x)/)') '<< CPFEM >> stress/MPa at el ip ', cp_en, IP, cauchyStress/1.0e6_pReal + write(6,'(a,i8,1x,i2,/,6(12x,6(f10.3,1x)/))') '<< CPFEM >> Jacobian/GPa at el ip ', cp_en, IP, transpose(jacobian)/1.0e9_pReal + flush(6) + !$OMP END CRITICAL (write2out) + endif + + + !*** warn if stiffness close to zero + + if (all(abs(jacobian) < 1e-10_pReal)) then + call IO_warning(601,cp_en,IP) + endif + + + !*** remember extreme values of stress and jacobian + + if (mode < 3) then + cauchyStress33 = math_Mandel6to33(cauchyStress) + if (maxval(cauchyStress33) > debug_stressMax) then + debug_stressMaxLocation = (/cp_en, IP/) + debug_stressMax = maxval(cauchyStress33) + endif + if (minval(cauchyStress33) < debug_stressMin) then + debug_stressMinLocation = (/cp_en, IP/) + debug_stressMin = minval(cauchyStress33) + endif + jacobian3333 = math_Mandel66to3333(jacobian) + if (maxval(jacobian3333) > debug_jacobianMax) then + debug_jacobianMaxLocation = (/cp_en, IP/) + debug_jacobianMax = maxval(jacobian3333) + endif + if (minval(jacobian3333) < debug_jacobianMin) then + debug_jacobianMinLocation = (/cp_en, IP/) + debug_jacobianMin = minval(jacobian3333) + endif + endif + end subroutine CPFEM_general -END MODULE CPFEM +end module CPFEM diff --git a/code/DAMASK_abaqus_exp.f b/code/DAMASK_abaqus_exp.f index 95bb1caa8..50f97455f 100644 --- a/code/DAMASK_abaqus_exp.f +++ b/code/DAMASK_abaqus_exp.f @@ -187,7 +187,13 @@ subroutine vumat (jblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, & debug_abaqus use mesh, only: mesh_FEasCP, & mesh_ipCoordinates - use CPFEM, only: CPFEM_general,CPFEM_init_done, CPFEM_initAll + use CPFEM, only: & + CPFEM_general, & + CPFEM_init_done, & + CPFEM_initAll, & + CPFEM_CALCRESULTS, & + CPFEM_AGERESULTS, & + CPFEM_EXPLICIT use homogenization, only: materialpoint_sizeResults, materialpoint_results include 'vaba_param.inc' ! Abaqus exp initializes a first step in single prec. for this a two-step compilation is used. @@ -216,8 +222,8 @@ subroutine vumat (jblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, & real(pReal), dimension(6,6) :: ddsdde real(pReal) temp, timeInc integer(pInt) computationMode, n, i, cp_en - logical :: cutBack + computationMode = ior(CPFEM_CALCRESULTS,CPFEM_EXPLICIT) ! always calculate, always explicit do n = 1,nblock ! loop over vector of IPs temp = tempOld(n) @@ -249,16 +255,15 @@ subroutine vumat (jblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, & outdatedByNewInc = .false. call debug_info() ! first after new inc reports debugging call debug_reset() ! resets debugging - computationMode = 8 ! calc and age results with implicit collection - else - computationMode = 9 ! plain calc with implicit collection + computationMode = ior(computationMode, CPFEM_AGERESULTS) ! age results endif theTime = totalTime ! record current starting time if (iand(debug_level(debug_abaqus),debug_levelBasic) /= 0) then !$OMP CRITICAL (write2out) - write(6,'(a16,1x,i2,1x,a,i8,1x,i5,a)') 'computationMode',computationMode,'(',nElement(n),nMatPoint(n),')'; call flush(6) + write(6,'(i2,1x,a,i8,1x,i5,a)') '(',nElement(n),nMatPoint(n),')'; call flush(6) + write(6,'(a,l1)') 'Aging Results: ', iand(computationMode, CPFEM_AGERESULTS) /= 0_pInt !$OMP END CRITICAL (write2out) endif diff --git a/code/DAMASK_abaqus_std.f b/code/DAMASK_abaqus_std.f index f4cd8b0b7..5b2a26472 100644 --- a/code/DAMASK_abaqus_std.f +++ b/code/DAMASK_abaqus_std.f @@ -151,7 +151,16 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,& debug_abaqus use mesh, only: mesh_FEasCP, & mesh_ipCoordinates - use CPFEM, only: CPFEM_general,CPFEM_init_done, CPFEM_initAll + use CPFEM, only: & + CPFEM_general, & + CPFEM_init_done, & + CPFEM_initAll, & + CPFEM_CALCRESULTS, & + CPFEM_AGERESULTS, & + CPFEM_COLLECT, & + CPFEM_RESTOREJACOBIAN, & + CPFEM_BACKUPJACOBIAN + use homogenization, only: materialpoint_sizeResults, materialpoint_results @@ -173,7 +182,6 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,& real(pReal), dimension(6) :: stress_h real(pReal), dimension(6,6) :: ddsdde_h integer(pInt) computationMode, i, cp_en - logical :: cutBack if (iand(debug_level(debug_abaqus),debug_levelBasic) /= 0 .and. noel == 1 .and. npt == 1) then !$OMP CRITICAL (write2out) @@ -216,9 +224,7 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,& !$OMP END CRITICAL (write2out) endif - else if ( dtime < theDelta ) then ! >> cutBack << - - cutBack = .true. + else if ( dtime < theDelta ) then ! >> cutBack << terminallyIll = .false. cycleCounter = -1 ! first calc step increments this to cycle = 0 calcMode = .true. ! pretend last step was calculation @@ -230,33 +236,27 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,& calcMode(npt,cp_en) = .not. calcMode(npt,cp_en) ! ping pong (calc <--> collect) - if ( calcMode(npt,cp_en) ) then ! now calc - if ( lastMode .neqv. calcMode(npt,cp_en) ) then ! first after ping pong - call debug_reset() ! resets debugging - outdatedFFN1 = .false. - cycleCounter = cycleCounter + 1 - endif - if ( outdatedByNewInc ) then - outdatedByNewInc = .false. - computationMode = 1 ! calc and age results - else - computationMode = 2 ! plain calc - endif + if (calcMode(npt,cp_en)) then ! now calc + computationMode = CPFEM_CALCRESULTS + if ( lastMode .neqv. calcMode(npt,cp_en) ) then ! first after ping pong + call debug_reset() ! resets debugging + outdatedFFN1 = .false. + cycleCounter = cycleCounter + 1 + endif + if(outdatedByNewInc) then + outdatedByNewInc = .false. + computationMode = ior(computationMode,CPFEM_AGERESULTS) ! calc and age results + endif else ! now collect - if ( lastMode .neqv. calcMode(npt,cp_en) .and. & - .not. terminallyIll) then - call debug_info() ! first after ping pong reports debugging - endif - if ( lastIncConverged ) then - lastIncConverged = .false. - computationMode = 4 ! collect and backup Jacobian after convergence - elseif ( cutBack ) then - cutBack = .false. - computationMode = 5 ! collect and restore Jacobian after cutback - else - computationMode = 3 ! plain collect - endif - mesh_ipCoordinates(1:3,npt,cp_en) = numerics_unitlength * COORDS + computationMode = CPFEM_COLLECT + if(lastMode .neqv. calcMode(npt,cp_en) .and. .not. terminallyIll) then + call debug_info() ! first after ping pong reports debugging + endif + if (lastIncConverged) then + lastIncConverged = .false. + computationMode = ior(computationMode,CPFEM_BACKUPJACOBIAN) ! backup Jacobian after convergence + endif + mesh_ipCoordinates(1:3,npt,cp_en) = numerics_unitlength * COORDS endif theTime = time(2) ! record current starting time diff --git a/code/DAMASK_marc.f90 b/code/DAMASK_marc.f90 index d8c58989f..b9f4152a2 100644 --- a/code/DAMASK_marc.f90 +++ b/code/DAMASK_marc.f90 @@ -16,44 +16,37 @@ ! You should have received a copy of the GNU General Public License ! along with DAMASK. If not, see . ! -!############################################################## -!* $Id$ -!******************************************************************** -! Material subroutine for MSC.Marc -! -! written by P. Eisenlohr, -! F. Roters, -! L. Hantcherli, -! W.A. Counts -! D.D. Tjahjanto -! C. Kords -! -! MPI fuer Eisenforschung, Duesseldorf -! -!******************************************************************** -! Usage: -! - choose material as hypela2 -! - set statevariable 2 to index of homogenization -! - set statevariable 3 to index of microstructure -! - make sure the file "material.config" exists in the working -! directory -! - make sure the file "numerics.config" exists in the working -! directory -! - use nonsymmetric option for solver (e.g. direct -! profile or multifrontal sparse, the latter seems -! to be faster!) -! - in case of ddm (domain decomposition)a SYMMETRIC -! solver has to be used, i.e uncheck "non-symmetric" -!******************************************************************** -! Marc subroutines used: -! - hypela2 -! - plotv -! - quit -!******************************************************************** -! Marc common blocks included: -! - concom: lovl, ncycle, inc, incsub -! - creeps: timinc -!******************************************************************** +!-------------------------------------------------------------------------------------------------- +! $Id$ +!-------------------------------------------------------------------------------------------------- +!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH +!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH +!> @author Luc Hantcherli, Max-Planck-Institut für Eisenforschung GmbH +!> @author W.A. Counts +!> @author Denny Tjahjanto, Max-Planck-Institut für Eisenforschung GmbH +!> @author Christoph Kords, Max-Planck-Institut für Eisenforschung GmbH +!> @brief Material subroutine for MSC.Marc +!> @details Usage: +!> @details - choose material as hypela2 +!> @details - set statevariable 2 to index of homogenization +!> @details - set statevariable 3 to index of microstructure +!> @details - make sure the file "material.config" exists in the working +!> @details directory +!> @details - make sure the file "numerics.config" exists in the working +!> @details directory +!> @details - use nonsymmetric option for solver (e.g. direct +!> @details profile or multifrontal sparse, the latter seems +!> @details to be faster!) +!> @details - in case of ddm (domain decomposition)a SYMMETRIC +!> @details solver has to be used, i.e uncheck "non-symmetric" +!> @details Marc subroutines used: +!> @details - hypela2 +!> @details - plotv +!> @details - quit +!> @details Marc common blocks included: +!> @details - concom: lovl, ncycle, inc, incsub +!> @details - creeps: timinc +!-------------------------------------------------------------------------------------------------- #ifndef INT #define INT 4 @@ -68,7 +61,6 @@ #include "prec.f90" module DAMASK_interface - use prec, only: pInt implicit none character(len=4), parameter :: InputFileExtension = '.dat' @@ -76,24 +68,29 @@ module DAMASK_interface contains + +!-------------------------------------------------------------------------------------------------- +!> @brief only output of current version +!-------------------------------------------------------------------------------------------------- subroutine DAMASK_interface_init implicit none !$OMP CRITICAL (write2out) - write(6,*) - write(6,*) '<<<+- DAMASK_marc init -+>>>' - write(6,*) '$Id$' + write(6,'(/,a)') ' <<<+- DAMASK_marc init -+>>>' + write(6,'(a)') ' $Id$' #include "compilation_info.f90" !$OMP END CRITICAL (write2out) end subroutine DAMASK_interface_init +!-------------------------------------------------------------------------------------------------- +!> @brief returns the current workingDir +!-------------------------------------------------------------------------------------------------- function getSolverWorkingDirectoryName() implicit none - character(1024) getSolverWorkingDirectoryName, inputName character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forward and backward slash @@ -105,10 +102,16 @@ function getSolverWorkingDirectoryName() end function getSolverWorkingDirectoryName -function getSolverJobName() - - implicit none +!-------------------------------------------------------------------------------------------------- +!> @brief solver job name (no extension) as combination of geometry and load case name +!-------------------------------------------------------------------------------------------------- +function getSolverJobName() + use prec, only: & + pReal, & + pInt + + implicit none character(1024) :: getSolverJobName, inputName character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forward and backward slash integer(pInt) :: extPos @@ -118,10 +121,10 @@ function getSolverJobName() inquire(5, name=inputName) ! determine inputfile extPos = len_trim(inputName)-4 getSolverJobName=inputName(scan(inputName,pathSep,back=.true.)+1:extPos) -! write(6,*) 'getSolverJobName', getSolverJobName end function getSolverJobName + end module DAMASK_interface #include "IO.f90" @@ -146,77 +149,70 @@ end module DAMASK_interface #include "CPFEM.f90" -!******************************************************************** -! This is the Marc material routine -!******************************************************************** -! -! ************* user subroutine for defining material behavior ************** -! -! -! CAUTION : Due to calculation of the Deformation gradients, Stretch Tensors and -! Rotation tensors at previous and current states, the analysis can be -! computationally expensive. Please use the user subroutine -> hypela -! if these kinematic quantities are not needed in the constitutive model -! -! -! IMPORTANT NOTES : -! -! (1) F,R,U are only available for continuum and membrane elements (not for -! shells and beams). -! -! (2) For total Lagrangian formulation use the -> 'Elasticity,1' card(= -! total Lagrange with large disp) in the parameter section of input deck. -! For updated Lagrangian formulation use the -> 'Plasticity,3' card(= -! update+finite+large disp+constant d) in the parameter section of -! input deck. -! -! The following operation obtains U (stretch tensor) at t=n+1 : -! -! call scla(un1,0.d0,itel,itel,1) -! do 3 k=1,3 -! do 2 i=1,3 -! do 1 j=1,3 -! un1(i,j)=un1(i,j)+dsqrt(strechn1(k))*eigvn1(i,k)*eigvn1(j,k) -!1 continue -!2 continue -!3 continue -! -!******************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief This is the MSC.Marc user subroutine for defining material behavior +!> @details CAUTION : Due to calculation of the Deformation gradients, Stretch Tensors and +!> @details Rotation tensors at previous and current states, the analysis can be +!> @details computationally expensive. Please use the user subroutine -> hypela +!> @details if these kinematic quantities are not needed in the constitutive model +!> @details +!> @details IMPORTANT NOTES : +!> @details +!> @details (1) F,R,U are only available for continuum and membrane elements (not for +!> @details shells and beams). +!> @details +!> @details (2) For total Lagrangian formulation use the -> 'Elasticity,1' card(= +!> @details total Lagrange with large disp) in the parameter section of input deck. +!> @details For updated Lagrangian formulation use the -> 'Plasticity,3' card(= +!> @details update+finite+large disp+constant d) in the parameter section of +!> @details input deck. +!> @details +!> @details The following operation obtains U (stretch tensor) at t=n+1 : +!> @details +!> @details call scla(un1,0.d0,itel,itel,1) +!> @details do k=1,3 +!> @details do i=1,3 +!> @details do j=1,3 +!> @details un1(i,j)=un1(i,j)+dsqrt(strechn1(k))*eigvn1(i,k)*eigvn1(j,k) +!> @details enddo +!> @details enddo +!> @details enddo +!-------------------------------------------------------------------------------------------------- subroutine hypela2(& - d,& ! stress strain law to be formed - g,& ! change in stress due to temperature effects - e,& ! total elastic strain - de,& ! increment of strain - s,& ! stress - should be updated by user - t,& ! state variables (comes in at t=n, must be updated to have state variables at t=n+1) - dt,& ! increment of state variables - ngens,& ! size of stress - strain law - n,& ! element number - nn,& ! integration point number - kcus,& ! (1) layer number, (2) internal layer number - matus,& ! (1) user material identification number, (2) internal material identification number - ndi,& ! number of direct components - nshear,& ! number of shear components - disp,& ! incremental displacements - dispt,& ! displacements at t=n (at assembly, lovl=4) and displacements at t=n+1 (at stress recovery, lovl=6) - coord,& ! coordinates - ffn,& ! deformation gradient - frotn,& ! rotation tensor - strechn,& ! square of principal stretch ratios, lambda(i) - eigvn,& ! i principal direction components for j eigenvalues - ffn1,& ! deformation gradient - frotn1,& ! rotation tensor - strechn1,& ! square of principal stretch ratios, lambda(i) - eigvn1,& ! i principal direction components for j eigenvalues - ncrd,& ! number of coordinates - itel,& ! dimension of F and R, either 2 or 3 - ndeg,& ! number of degrees of freedom ==> is this at correct list position ?!? - ndm,& ! - nnode,& ! number of nodes per element - jtype,& ! element type - lclass,& ! element class - ifr,& ! set to 1 if R has been calculated - ifu & ! set to 1 if stretch has been calculated + d,& !< stress strain law to be formed + g,& !< change in stress due to temperature effects + e,& !< total elastic strain + de,& !< increment of strain + s,& !< stress - should be updated by user + t,& !< state variables (comes in at t=n, must be updated to have state variables at t=n+1) + dt,& !< increment of state variables + ngens,& !< size of stress - strain law + n,& !< element number + nn,& !< integration point number + kcus,& !< (1) layer number, (2) internal layer number + matus,& !< (1) user material identification number, (2) internal material identification number + ndi,& !< number of direct components + nshear,& !< number of shear components + disp,& !< incremental displacements + dispt,& !< displacements at t=n (at assembly, lovl=4) and displacements at t=n+1 (at stress recovery, lovl=6) + coord,& !< coordinates + ffn,& !< deformation gradient + frotn,& !< rotation tensor + strechn,& !< square of principal stretch ratios, lambda(i) + eigvn,& !< i principal direction components for j eigenvalues + ffn1,& !< deformation gradient + frotn1,& !< rotation tensor + strechn1,& !< square of principal stretch ratios, lambda(i) + eigvn1,& !< i principal direction components for j eigenvalues + ncrd,& !< number of coordinates + itel,& !< dimension of F and R, either 2 or 3 + ndeg,& !< number of degrees of freedom ==> is this at correct list position ?!? + ndm,& !< + nnode,& !< number of nodes per element + jtype,& !< element type + lclass,& !< element class + ifr,& !< set to 1 if R has been calculated + ifu & !< set to 1 if stretch has been calculated ) use prec, only: pReal, & @@ -244,7 +240,16 @@ subroutine hypela2(& mesh_build_ipCoordinates, & FE_Nnodes, & FE_geomtype - use CPFEM, only: CPFEM_initAll,CPFEM_general,CPFEM_init_done + use CPFEM, only: & + CPFEM_general, & + CPFEM_init_done, & + CPFEM_initAll, & + CPFEM_CALCRESULTS, & + CPFEM_AGERESULTS, & + CPFEM_COLLECT, & + CPFEM_RESTOREJACOBIAN, & + CPFEM_BACKUPJACOBIAN + !$ use numerics, only: DAMASK_NumThreadsInt ! number of threads set by DAMASK_NUM_THREADS implicit none @@ -286,12 +291,9 @@ subroutine hypela2(& !$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS - if (lovl == 4) then ! Marc requires stiffness in separate call (lovl == 4) - if ( timinc < theDelta .and. theInc == inc ) then ! first after cutback - computationMode = 7 ! --> restore tangent and return it - else - computationMode = 6 ! --> just return known tangent - endif + if (lovl == 4 ) then + if(timinc < theDelta .and. theInc == inc ) & ! first after cutback + computationMode = CPFEM_RESTOREJACOBIAN else ! stress requested (lovl == 6) cp_en = mesh_FEasCP('elem',n(1)) if (cptim > theTime .or. inc /= theInc) then ! reached "convergence" @@ -346,10 +348,10 @@ subroutine hypela2(& call mesh_build_ipCoordinates() ! update ip coordinates endif if ( outdatedByNewInc ) then + computationMode = ior(CPFEM_CALCRESULTS,CPFEM_AGERESULTS) outdatedByNewInc = .false. ! reset flag - computationMode = 1 ! calc and age results else - computationMode = 2 ! plain calc + computationMode = CPFEM_CALCRESULTS endif else ! now --- COLLECT --- if ( lastMode /= calcMode(nn,cp_en) .and. & @@ -357,10 +359,10 @@ subroutine hypela2(& call debug_info() ! first after ping pong reports (meaningful) debugging endif if ( lastIncConverged ) then + computationMode = ior(CPFEM_COLLECT,CPFEM_BACKUPJACOBIAN) ! collect and backup Jacobian after convergence lastIncConverged = .false. ! reset flag - computationMode = 4 ! collect and backup Jacobian after convergence else - computationMode = 3 ! plain collect + computationMode = CPFEM_COLLECT ! plain collect endif do node = 1,FE_Nnodes(FE_geomtype(mesh_element(2,cp_en))) FEnodeID = mesh_FEasCP('node',mesh_element(4+node,cp_en)) @@ -389,27 +391,24 @@ subroutine hypela2(& end subroutine hypela2 -!******************************************************************** -! This routine sets user defined output variables for Marc -!******************************************************************** -! -! select a variable contour plotting (user subroutine). -! -!******************************************************************** +!-------------------------------------------------------------------------------------------------- +!> @brief sets user defined output variables for Marc +!> @details select a variable contour plotting (user subroutine). +!-------------------------------------------------------------------------------------------------- subroutine plotv(& - v,& ! variable - s,& ! stress array - sp,& ! stresses in preferred direction - etot,& ! total strain (generalized) - eplas,& ! total plastic strain - ecreep,& ! total creep strain - t,& ! current temperature - m,& ! element number - nn,& ! integration point number - layer,& ! layer number - ndi,& ! number of direct stress components - nshear,& ! number of shear stress components - jpltcd & ! user variable index + v,& !< variable + s,& !< stress array + sp,& !< stresses in preferred direction + etot,& !< total strain (generalized) + eplas,& !< total plastic strain + ecreep,& !< total creep strain + t,& !< current temperature + m,& !< element number + nn,& !< integration point number + layer,& !< layer number + ndi,& !< number of direct stress components + nshear,& !< number of shear stress components + jpltcd & !< user variable index ) use prec, only: pReal,pInt use mesh, only: mesh_FEasCP @@ -421,7 +420,7 @@ subroutine plotv(& real(pReal) v, t(*) integer(pInt) m, nn, layer, ndi, nshear, jpltcd - if (jpltcd > materialpoint_sizeResults) call IO_error(700_pInt,jpltcd) ! complain about out of bounds error + if (jpltcd > materialpoint_sizeResults) call IO_error(700_pInt,jpltcd) ! complain about out of bounds error v = materialpoint_results(jpltcd,nn,mesh_FEasCP('elem', m)) diff --git a/code/DAMASK_spectral_utilities.f90 b/code/DAMASK_spectral_utilities.f90 index a4dca573c..da666df0a 100644 --- a/code/DAMASK_spectral_utilities.f90 +++ b/code/DAMASK_spectral_utilities.f90 @@ -691,7 +691,12 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,& wgt, & mesh_NcpElems use CPFEM, only: & - CPFEM_general + CPFEM_general, & + CPFEM_COLLECT, & + CPFEM_CALCRESULTS, & + CPFEM_AGERESULTS, & + CPFEM_BACKUPJACOBIAN, & + CPFEM_RESTOREJACOBIAN use homogenization, only: & materialpoint_F0, & materialpoint_F, & @@ -720,16 +725,16 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,& real(pReal), dimension(6,6) :: dsde !< d sigma / d Epsilon write(6,'(/,a)') ' ... evaluating constitutive response ......................................' + calcMode = CPFEM_CALCRESULTS + collectMode = CPFEM_COLLECT if (forwardData) then ! aging results - calcMode = 1_pInt - collectMode = 4_pInt - else ! normal calculation - calcMode = 2_pInt - collectMode = 3_pInt + calcMode = ior(calcMode, CPFEM_AGERESULTS) + collectMode = ior(collectMode, CPFEM_BACKUPJACOBIAN) endif if (cutBack) then ! restore saved variables - calcMode = 2_pInt - collectMode = 5_pInt + collectMode = ior(collectMode , CPFEM_RESTOREJACOBIAN) + collectMode = iand(collectMode, not(CPFEM_BACKUPJACOBIAN)) + calcMode = iand(calcMode, not(CPFEM_AGERESULTS)) endif !-------------------------------------------------------------------------------------------------- ! calculate bounds of det(F) and report diff --git a/code/IO.f90 b/code/IO.f90 index f44c4cc6c..e10633077 100644 --- a/code/IO.f90 +++ b/code/IO.f90 @@ -1515,9 +1515,9 @@ subroutine IO_error(error_ID,e,i,g,ext_msg) !-------------------------------------------------------------------------------------------------- ! user errors case (600_pInt) - msg = 'Non-local plasticity and non-CP elements in model' + msg = 'Cannot combine Non-local plasticity and non-DAMASK elements' case (601_pInt) - msg = 'OpenMP threads > 1 and using non-CP elements' + msg = 'Cannot combine OpenMP threading and non-DAMASK elements' !------------------------------------------------------------------------------------------------- ! DAMASK_marc errors diff --git a/code/crystallite.f90 b/code/crystallite.f90 index e843c872f..1636a6cf8 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -89,9 +89,9 @@ module crystallite logical, dimension (:,:,:), allocatable, public :: & crystallite_requested !< flag to request crystallite calculation logical, dimension (:,:,:), allocatable, public, protected :: & - crystallite_converged !< convergence flag + crystallite_converged, & !< convergence flag + crystallite_localPlasticity !< indicates this grain to have purely local constitutive law logical, dimension (:,:,:), allocatable, private :: & - crystallite_localPlasticity, & !< indicates this grain to have purely local constitutive law crystallite_todo !< flag to indicate need for further computation logical, dimension (:,:), allocatable, private :: & crystallite_clearToWindForward, &