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

Merge branch 'Marc-no-pingpong' into 'development' 

Marc no pingpong

See merge request damask/DAMASK!180
This commit is contained in:
Franz Roters 2020-06-18 15:29:21 +02:00
parent 0e11d6f049 06f6e15123
commit e959aaab5d
7 changed files with 175 additions and 368 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

138
src/CPFEM.f90
View File

@ -35,29 +35,15 @@ module CPFEM
CPFEM_dcsdE !< Cauchy stress tangent CPFEM_dcsdE !< Cauchy stress tangent
real(pReal), dimension (:,:,:,:), allocatable, private :: & real(pReal), dimension (:,:,:,:), allocatable, private :: &
CPFEM_dcsdE_knownGood !< known good tangent CPFEM_dcsdE_knownGood !< known good tangent
integer(pInt), public :: &
cycleCounter = 0_pInt, & !< needs description
theInc = -1_pInt, & !< needs description
lastLovl = 0_pInt !< lovl in previous call to marc hypela2
real(pReal), public :: &
theTime = 0.0_pReal, & !< needs description
theDelta = 0.0_pReal
logical, public :: &
outdatedFFN1 = .false., & !< needs description
lastIncConverged = .false., & !< needs description
outdatedByNewInc = .false. !< needs description
logical, public, protected :: & integer(pInt), public :: &
CPFEM_init_done = .false. !< remember whether init has been done already cycleCounter = 0_pInt !< needs description
logical, private :: &
CPFEM_calc_done = .false. !< remember whether first ip has already calced the results
integer(pInt), parameter, public :: & integer(pInt), parameter, public :: &
CPFEM_COLLECT = 2_pInt**0_pInt, & CPFEM_CALCRESULTS = 2_pInt**0_pInt, &
CPFEM_CALCRESULTS = 2_pInt**1_pInt, & CPFEM_AGERESULTS = 2_pInt**1_pInt, &
CPFEM_AGERESULTS = 2_pInt**2_pInt, & CPFEM_BACKUPJACOBIAN = 2_pInt**2_pInt, &
CPFEM_BACKUPJACOBIAN = 2_pInt**3_pInt, & CPFEM_RESTOREJACOBIAN = 2_pInt**3_pInt
CPFEM_RESTOREJACOBIAN = 2_pInt**4_pInt
public :: & public :: &
CPFEM_general, & CPFEM_general, &
@ -68,14 +54,10 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief call (thread safe) all module initializations !> @brief call all module initializations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine CPFEM_initAll(el,ip) subroutine CPFEM_initAll
integer(pInt), intent(in) :: el, & !< FE el number
ip !< FE integration point number
CPFEM_init_done = .true.
call DAMASK_interface_init call DAMASK_interface_init
call prec_init call prec_init
call IO_init call IO_init
@ -88,7 +70,7 @@ subroutine CPFEM_initAll(el,ip)
call YAML_init call YAML_init
call HDF5_utilities_init call HDF5_utilities_init
call results_init(.false.) call results_init(.false.)
call discretization_marc_init(ip, el) call discretization_marc_init
call lattice_init call lattice_init
call material_init(.false.) call material_init(.false.)
call constitutive_init call constitutive_init
@ -124,7 +106,7 @@ end subroutine CPFEM_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief perform initialization at first call, update variables and call the actual material model !> @brief perform initialization at first call, update variables and call the actual material model
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyStress, jacobian) subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyStress, jacobian)
integer(pInt), intent(in) :: elFE, & !< FE element number integer(pInt), intent(in) :: elFE, & !< FE element number
ip !< integration point number ip !< integration point number
@ -133,24 +115,21 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
ffn1 !< deformation gradient for t=t1 ffn1 !< deformation gradient for t=t1
integer(pInt), intent(in) :: mode !< computation mode 1: regular computation plus aging of results integer(pInt), intent(in) :: mode !< computation mode 1: regular computation plus aging of results
real(pReal), intent(in) :: temperature_inp !< temperature real(pReal), intent(in) :: temperature_inp !< temperature
logical, intent(in) :: parallelExecution !< flag indicating parallel computation of requested IPs
real(pReal), dimension(6), intent(out) :: cauchyStress !< stress as 6 vector 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) real(pReal), dimension(6,6), intent(out) :: jacobian !< jacobian as 66 tensor (Consistent tangent dcs/dE)
real(pReal) J_inverse, & ! inverse of Jacobian real(pReal) J_inverse, & ! inverse of Jacobian
rnd rnd
real(pReal), dimension (3,3) :: Kirchhoff, & ! Piola-Kirchhoff stress real(pReal), dimension (3,3) :: Kirchhoff ! Piola-Kirchhoff stress
cauchyStress33 ! stress vector
real(pReal), dimension (3,3,3,3) :: H_sym, & real(pReal), dimension (3,3,3,3) :: H_sym, &
H, & H
jacobian3333 ! jacobian in Matrix notation
integer(pInt) elCP, & ! crystal plasticity element number integer(pInt) elCP, & ! crystal plasticity element number
i, j, k, l, m, n, ph, homog, mySource i, j, k, l, m, n, ph, homog, mySource
logical updateJaco ! flag indicating if Jacobian has to be updated logical updateJaco ! flag indicating if Jacobian has to be updated
real(pReal), parameter :: ODD_STRESS = 1e15_pReal, & !< return value for stress in case of ping pong dummy cycle real(pReal), parameter :: ODD_STRESS = 1e15_pReal, & !< return value for stress if terminallyIll
ODD_JACOBIAN = 1e50_pReal !< return value for jacobian in case of ping pong dummy cycle ODD_JACOBIAN = 1e50_pReal !< return value for jacobian if terminallyIll
elCP = mesh_FEM2DAMASK_elem(elFE) elCP = mesh_FEM2DAMASK_elem(elFE)
@ -159,9 +138,6 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
write(6,'(/,a)') '#############################################' write(6,'(/,a)') '#############################################'
write(6,'(a1,a22,1x,i8,a13)') '#','element', elCP, '#' write(6,'(a1,a22,1x,i8,a13)') '#','element', elCP, '#'
write(6,'(a1,a22,1x,i8,a13)') '#','ip', ip, '#' write(6,'(a1,a22,1x,i8,a13)') '#','ip', ip, '#'
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)') '#','cycleCounter', cycleCounter, '#'
write(6,'(a1,a22,1x,i8,a13)') '#','computationMode',mode, '#' write(6,'(a1,a22,1x,i8,a13)') '#','computationMode',mode, '#'
if (terminallyIll) & if (terminallyIll) &
@ -174,13 +150,8 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
if (iand(mode, CPFEM_RESTOREJACOBIAN) /= 0_pInt) & if (iand(mode, CPFEM_RESTOREJACOBIAN) /= 0_pInt) &
CPFEM_dcsde = CPFEM_dcsde_knownGood CPFEM_dcsde = CPFEM_dcsde_knownGood
!*** age results
if (iand(mode, CPFEM_AGERESULTS) /= 0_pInt) call CPFEM_forward if (iand(mode, CPFEM_AGERESULTS) /= 0_pInt) call CPFEM_forward
!*** collection of FEM input with returning of randomize odd stress and jacobian
!* If no parallel execution is required, there is no need to collect FEM input
if (.not. parallelExecution) then
chosenThermal1: select case (thermal_type(material_homogenizationAt(elCP))) chosenThermal1: select case (thermal_type(material_homogenizationAt(elCP)))
case (THERMAL_conduction_ID) chosenThermal1 case (THERMAL_conduction_ID) chosenThermal1
temperature(material_homogenizationAt(elCP))%p(thermalMapping(material_homogenizationAt(elCP))%p(ip,elCP)) = & temperature(material_homogenizationAt(elCP))%p(thermalMapping(material_homogenizationAt(elCP))%p(ip,elCP)) = &
@ -189,64 +160,22 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
materialpoint_F0(1:3,1:3,ip,elCP) = ffn materialpoint_F0(1:3,1:3,ip,elCP) = ffn
materialpoint_F(1:3,1:3,ip,elCP) = ffn1 materialpoint_F(1:3,1:3,ip,elCP) = ffn1
elseif (iand(mode, CPFEM_COLLECT) /= 0_pInt) then
call random_number(rnd)
if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal
CPFEM_cs(1:6,ip,elCP) = rnd * ODD_STRESS
CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_identity2nd(6)
chosenThermal2: select case (thermal_type(material_homogenizationAt(elCP)))
case (THERMAL_conduction_ID) chosenThermal2
temperature(material_homogenizationAt(elCP))%p(thermalMapping(material_homogenizationAt(elCP))%p(ip,elCP)) = &
temperature_inp
end select chosenThermal2
materialpoint_F0(1:3,1:3,ip,elCP) = ffn
materialpoint_F(1:3,1:3,ip,elCP) = ffn1
CPFEM_calc_done = .false.
endif
!*** calculation of stress and jacobian
if (iand(mode, CPFEM_CALCRESULTS) /= 0_pInt) 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 validCalculation: if (terminallyIll) then
validCalculation: if (terminallyIll &
.or. outdatedFFN1 &
.or. any(abs(ffn1 - materialpoint_F(1:3,1:3,ip,elCP)) > defgradTolerance)) then
if (any(abs(ffn1 - materialpoint_F(1:3,1:3,ip,elCP)) > defgradTolerance)) then
if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) then
write(6,'(a,1x,i8,1x,i2)') '<< CPFEM >> OUTDATED at elFE ip',elFE,ip
write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 old:',&
transpose(materialpoint_F(1:3,1:3,ip,elCP))
write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 now:',transpose(ffn1)
endif
outdatedFFN1 = .true.
endif
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
CPFEM_cs(1:6,ip,elCP) = ODD_STRESS * rnd CPFEM_cs(1:6,ip,elCP) = ODD_STRESS * rnd
CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_identity2nd(6) CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_identity2nd(6)
!*** deformation gradient is not outdated
else validCalculation else validCalculation
updateJaco = mod(cycleCounter,iJacoStiffness) == 0 updateJaco = mod(cycleCounter,iJacoStiffness) == 0
!* no parallel computation, so we use just one single elFE and ip for computation FEsolving_execElem = elCP
if (.not. parallelExecution) then FEsolving_execIP = ip
FEsolving_execElem = elCP if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) &
FEsolving_execIP = ip write(6,'(a,i8,1x,i2)') '<< CPFEM >> calculation for elFE ip ',elFE,ip
if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) & call materialpoint_stressAndItsTangent(updateJaco, dt)
write(6,'(a,i8,1x,i2)') '<< CPFEM >> calculation for elFE ip ',elFE,ip
call materialpoint_stressAndItsTangent(updateJaco, dt)
!* parallel computation and calulation not yet done
elseif (.not. CPFEM_calc_done) then
if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) &
write(6,'(a,i8,a,i8)') '<< CPFEM >> calculation for elements ',FEsolving_execElem(1),&
' to ',FEsolving_execElem(2)
call materialpoint_stressAndItsTangent(updateJaco, dt)
CPFEM_calc_done = .true.
endif
!* map stress and stiffness (or return odd values if terminally ill)
terminalIllness: if (terminallyIll) then terminalIllness: if (terminallyIll) then
call random_number(rnd) call random_number(rnd)
@ -256,7 +185,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
else terminalIllness else terminalIllness
! translate from P to CS ! translate from P to sigma
Kirchhoff = matmul(materialpoint_P(1:3,1:3,ip,elCP), transpose(materialpoint_F(1:3,1:3,ip,elCP))) Kirchhoff = matmul(materialpoint_P(1:3,1:3,ip,elCP), transpose(materialpoint_F(1:3,1:3,ip,elCP)))
J_inverse = 1.0_pReal / math_det33(materialpoint_F(1:3,1:3,ip,elCP)) J_inverse = 1.0_pReal / math_det33(materialpoint_F(1:3,1:3,ip,elCP))
CPFEM_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.) CPFEM_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.)
@ -280,7 +209,6 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
endif terminalIllness endif terminalIllness
endif validCalculation endif validCalculation
!* report stress and stiffness
if ((iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) & if ((iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) &
.and. ((debug_e == elCP .and. debug_i == ip) & .and. ((debug_e == elCP .and. debug_i == ip) &
.or. .not. iand(debug_level(debug_CPFEM), debug_levelSelective) /= 0_pInt)) then .or. .not. iand(debug_level(debug_CPFEM), debug_levelSelective) /= 0_pInt)) then
@ -293,35 +221,11 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
endif endif
!*** warn if stiffness close to zero
if (all(abs(CPFEM_dcsdE(1:6,1:6,ip,elCP)) < 1e-10_pReal)) call IO_warning(601,elCP,ip) if (all(abs(CPFEM_dcsdE(1:6,1:6,ip,elCP)) < 1e-10_pReal)) call IO_warning(601,elCP,ip)
!*** copy to output if using commercial FEM solver
cauchyStress = CPFEM_cs (1:6, ip,elCP) cauchyStress = CPFEM_cs (1:6, ip,elCP)
jacobian = CPFEM_dcsdE(1:6,1:6,ip,elCP) jacobian = CPFEM_dcsdE(1:6,1:6,ip,elCP)
!*** remember extreme values of stress ...
cauchyStress33 = math_6toSym33(CPFEM_cs(1:6,ip,elCP),weighted=.false.)
if (maxval(cauchyStress33) > debug_stressMax) then
debug_stressMaxLocation = [elCP, ip]
debug_stressMax = maxval(cauchyStress33)
endif
if (minval(cauchyStress33) < debug_stressMin) then
debug_stressMinLocation = [elCP, ip]
debug_stressMin = minval(cauchyStress33)
endif
!*** ... and Jacobian
jacobian3333 = math_66toSym3333(CPFEM_dcsdE(1:6,1:6,ip,elCP),weighted=.false.)
if (maxval(jacobian3333) > debug_jacobianMax) then
debug_jacobianMaxLocation = [elCP, ip]
debug_jacobianMax = maxval(jacobian3333)
endif
if (minval(jacobian3333) < debug_jacobianMin) then
debug_jacobianMinLocation = [elCP, ip]
debug_jacobianMin = minval(jacobian3333)
endif
end subroutine CPFEM_general end subroutine CPFEM_general

323
src/DAMASK_marc.f90
View File

@ -43,8 +43,6 @@ module DAMASK_interface
logical, protected, public :: symmetricSolver logical, protected, public :: symmetricSolver
character(len=*), parameter, public :: INPUTFILEEXTENSION = '.dat' character(len=*), parameter, public :: INPUTFILEEXTENSION = '.dat'
logical, dimension(:,:), public, allocatable :: &
calcMode !< calculate or collect (ping pong scheme)
public :: & public :: &
DAMASK_interface_init, & DAMASK_interface_init, &
@ -175,62 +173,62 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
dispt,coord,ffn,frotn,strechn,eigvn,ffn1,frotn1, & dispt,coord,ffn,frotn,strechn,eigvn,ffn1,frotn1, &
strechn1,eigvn1,ncrd,itel,ndeg,ndm,nnode, & strechn1,eigvn1,ncrd,itel,ndeg,ndm,nnode, &
jtype,lclass,ifr,ifu) jtype,lclass,ifr,ifu)
use prec use prec
use DAMASK_interface use DAMASK_interface
use numerics use numerics
use FEsolving use FEsolving
use debug use debug
use discretization_marc use discretization_marc
use homogenization use homogenization
use CPFEM use CPFEM
implicit none implicit none
!$ include "omp_lib.h" ! the openMP function library include "omp_lib.h" ! the openMP function library
integer, intent(in) :: & ! according to MSC.Marc 2012 Manual D integer, intent(in) :: & ! according to MSC.Marc 2012 Manual D
ngens, & !< size of stress-strain law ngens, & !< size of stress-strain law
nn, & !< integration point number nn, & !< integration point number
ndi, & !< number of direct components ndi, & !< number of direct components
nshear, & !< number of shear components nshear, & !< number of shear components
ncrd, & !< number of coordinates ncrd, & !< number of coordinates
itel, & !< dimension of F and R, either 2 or 3 itel, & !< dimension of F and R, either 2 or 3
ndeg, & !< number of degrees of freedom ndeg, & !< number of degrees of freedom
ndm, & !< not specified in MSC.Marc 2012 Manual D ndm, & !< not specified in MSC.Marc 2012 Manual D
nnode, & !< number of nodes per element nnode, & !< number of nodes per element
jtype, & !< element type jtype, & !< element type
ifr, & !< set to 1 if R has been calculated ifr, & !< set to 1 if R has been calculated
ifu !< set to 1 if stretch has been calculated ifu !< set to 1 if stretch has been calculated
integer, dimension(2), intent(in) :: & ! according to MSC.Marc 2012 Manual D integer, dimension(2), intent(in) :: & ! according to MSC.Marc 2012 Manual D
m, & !< (1) user element number, (2) internal element number m, & !< (1) user element number, (2) internal element number
matus, & !< (1) user material identification number, (2) internal material identification number matus, & !< (1) user material identification number, (2) internal material identification number
kcus, & !< (1) layer number, (2) internal layer number kcus, & !< (1) layer number, (2) internal layer number
lclass !< (1) element class, (2) 0: displacement, 1: low order Herrmann, 2: high order Herrmann lclass !< (1) element class, (2) 0: displacement, 1: low order Herrmann, 2: high order Herrmann
real(pReal), dimension(*), intent(in) :: & ! has dimension(1) according to MSC.Marc 2012 Manual D, but according to example hypela2.f dimension(*) real(pReal), dimension(*), intent(in) :: & ! has dimension(1) according to MSC.Marc 2012 Manual D, but according to example hypela2.f dimension(*)
e, & !< total elastic strain e, & !< total elastic strain
de, & !< increment of strain de, & !< increment of strain
dt !< increment of state variables dt !< increment of state variables
real(pReal), dimension(itel), intent(in) :: & ! according to MSC.Marc 2012 Manual D real(pReal), dimension(itel), intent(in) :: & ! according to MSC.Marc 2012 Manual D
strechn, & !< square of principal stretch ratios, lambda(i) at t=n strechn, & !< square of principal stretch ratios, lambda(i) at t=n
strechn1 !< square of principal stretch ratios, lambda(i) at t=n+1 strechn1 !< square of principal stretch ratios, lambda(i) at t=n+1
real(pReal), dimension(3,3), intent(in) :: & ! has dimension(itel,*) according to MSC.Marc 2012 Manual D, but we alway assume dimension(3,3) real(pReal), dimension(3,3), intent(in) :: & ! has dimension(itel,*) according to MSC.Marc 2012 Manual D, but we alway assume dimension(3,3)
ffn, & !< deformation gradient at t=n ffn, & !< deformation gradient at t=n
ffn1 !< deformation gradient at t=n+1 ffn1 !< deformation gradient at t=n+1
real(pReal), dimension(itel,*), intent(in) :: & ! according to MSC.Marc 2012 Manual D real(pReal), dimension(itel,*), intent(in) :: & ! according to MSC.Marc 2012 Manual D
frotn, & !< rotation tensor at t=n frotn, & !< rotation tensor at t=n
eigvn, & !< i principal direction components for j eigenvalues at t=n eigvn, & !< i principal direction components for j eigenvalues at t=n
frotn1, & !< rotation tensor at t=n+1 frotn1, & !< rotation tensor at t=n+1
eigvn1 !< i principal direction components for j eigenvalues at t=n+1 eigvn1 !< i principal direction components for j eigenvalues at t=n+1
real(pReal), dimension(ndeg,*), intent(in) :: & ! according to MSC.Marc 2012 Manual D real(pReal), dimension(ndeg,*), intent(in) :: & ! according to MSC.Marc 2012 Manual D
disp, & !< incremental displacements disp, & !< incremental displacements
dispt !< displacements at t=n (at assembly, lovl=4) and displacements at t=n+1 (at stress recovery, lovl=6) dispt !< displacements at t=n (at assembly, lovl=4) and displacements at t=n+1 (at stress recovery, lovl=6)
real(pReal), dimension(ncrd,*), intent(in) :: & ! according to MSC.Marc 2012 Manual D real(pReal), dimension(ncrd,*), intent(in) :: & ! according to MSC.Marc 2012 Manual D
coord !< coordinates coord !< coordinates
real(pReal), dimension(*), intent(inout) :: & ! according to MSC.Marc 2012 Manual D real(pReal), dimension(*), intent(inout) :: & ! according to MSC.Marc 2012 Manual D
t !< state variables (comes in at t=n, must be updated to have state variables at t=n+1) t !< state variables (comes in at t=n, must be updated to have state variables at t=n+1)
real(pReal), dimension(ndi+nshear), intent(out) :: & ! has dimension(*) according to MSC.Marc 2012 Manual D, but we need to loop over it real(pReal), dimension(ndi+nshear), intent(out) :: & ! has dimension(*) according to MSC.Marc 2012 Manual D, but we need to loop over it
s, & !< stress - should be updated by user s, & !< stress - should be updated by user
g !< change in stress due to temperature effects g !< change in stress due to temperature effects
real(pReal), dimension(ngens,ngens), intent(out) :: & ! according to MSC.Marc 2012 Manual D, but according to example hypela2.f dimension(ngens,*) real(pReal), dimension(ngens,ngens), intent(out) :: & ! according to MSC.Marc 2012 Manual D, but according to example hypela2.f dimension(ngens,*)
d !< stress-strain law to be formed d !< stress-strain law to be formed
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Marc common blocks are in fixed format so they have to be reformated to free format (f90) ! Marc common blocks are in fixed format so they have to be reformated to free format (f90)
@ -239,137 +237,108 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
#include QUOTE(PASTE(./marc/include/concom,Marc4DAMASK)) ! concom is needed for inc, lovl #include QUOTE(PASTE(./marc/include/concom,Marc4DAMASK)) ! concom is needed for inc, lovl
#include QUOTE(PASTE(./marc/include/creeps,Marc4DAMASK)) ! creeps is needed for timinc (time increment) #include QUOTE(PASTE(./marc/include/creeps,Marc4DAMASK)) ! creeps is needed for timinc (time increment)
logical :: cutBack logical :: cutBack
real(pReal), dimension(6) :: stress real(pReal), dimension(6) :: stress
real(pReal), dimension(6,6) :: ddsdde real(pReal), dimension(6,6) :: ddsdde
integer :: computationMode, i, cp_en, node, CPnodeID integer :: computationMode, i, cp_en, node, CPnodeID
!$ integer(4) :: defaultNumThreadsInt !< default value set by Marc integer(4) :: defaultNumThreadsInt !< default value set by Marc
if (iand(debug_level(debug_MARC),debug_LEVELBASIC) /= 0) then integer(pInt), save :: &
write(6,'(a,/,i8,i8,i2)') ' MSC.MARC information on shape of element(2), IP:', m, nn theInc = -1_pInt, & !< needs description
write(6,'(a,2(i1))') ' Jacobian: ', ngens,ngens lastLovl = 0_pInt !< lovl in previous call to marc hypela2
write(6,'(a,i1)') ' Direct stress: ', ndi real(pReal), save :: &
write(6,'(a,i1)') ' Shear stress: ', nshear theTime = 0.0_pReal, & !< needs description
write(6,'(a,i2)') ' DoF: ', ndeg theDelta = 0.0_pReal
write(6,'(a,i2)') ' Coordinates: ', ncrd logical, save :: &
write(6,'(a,i12)') ' Nodes: ', nnode lastIncConverged = .false., & !< needs description
write(6,'(a,i1)') ' Deformation gradient: ', itel outdatedByNewInc = .false., & !< needs description
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n:', & CPFEM_init_done = .false. !< remember whether init has been done already
transpose(ffn)
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n+1:', &
transpose(ffn1)
endif
!$ defaultNumThreadsInt = omp_get_num_threads() ! remember number of threads set by Marc if (iand(debug_level(debug_MARC),debug_LEVELBASIC) /= 0) then
!$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS write(6,'(a,/,i8,i8,i2)') ' MSC.MARC information on shape of element(2), IP:', m, nn
write(6,'(a,2(i1))') ' Jacobian: ', ngens,ngens
write(6,'(a,i1)') ' Direct stress: ', ndi
write(6,'(a,i1)') ' Shear stress: ', nshear
write(6,'(a,i2)') ' DoF: ', ndeg
write(6,'(a,i2)') ' Coordinates: ', ncrd
write(6,'(a,i12)') ' Nodes: ', nnode
write(6,'(a,i1)') ' Deformation gradient: ', itel
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n:', &
transpose(ffn)
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n+1:', &
transpose(ffn1)
endif
if (.not. CPFEM_init_done) call CPFEM_initAll(m(1),nn) defaultNumThreadsInt = omp_get_num_threads() ! remember number of threads set by Marc
call omp_set_num_threads(1) ! no openMP
computationMode = 0 ! save initialization value, since it does not result in any calculation if (.not. CPFEM_init_done) then
if (lovl == 4 ) then ! jacobian requested by marc CPFEM_init_done = .true.
if (timinc < theDelta .and. theInc == inc .and. lastLovl /= lovl) & ! first after cutback call CPFEM_initAll
computationMode = CPFEM_RESTOREJACOBIAN endif
elseif (lovl == 6) then ! stress requested by marc
cp_en = mesh_FEM2DAMASK_elem(m(1))
if (cptim > theTime .or. inc /= theInc) then ! reached "convergence"
terminallyIll = .false.
cycleCounter = -1 ! first calc step increments this to cycle = 0
if (inc == 0) then ! >> start of analysis <<
lastIncConverged = .false. ! no Jacobian backup
outdatedByNewInc = .false. ! no aging of state
calcMode = .false. ! pretend last step was collection
lastLovl = lovl ! pretend that this is NOT the first after a lovl change
write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> start of analysis..! ',m(1),nn
flush(6)
else if (inc - theInc > 1) then ! >> restart of broken analysis <<
lastIncConverged = .false. ! no Jacobian backup
outdatedByNewInc = .false. ! no aging of state
calcMode = .true. ! pretend last step was calculation
write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> restart of analysis..! ',m(1),nn
flush(6)
else ! >> just the next inc <<
lastIncConverged = .true. ! request Jacobian backup
outdatedByNewInc = .true. ! request aging of state
calcMode = .true. ! assure last step was calculation
write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> new increment..! ',m(1),nn
flush(6)
endif
else if ( timinc < theDelta ) then ! >> cutBack <<
lastIncConverged = .false. ! no Jacobian backup
outdatedByNewInc = .false. ! no aging of state
terminallyIll = .false.
cycleCounter = -1 ! first calc step increments this to cycle = 0
calcMode = .true. ! pretend last step was calculation
write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> cutback detected..! ',m(1),nn
flush(6)
endif ! convergence treatment end
computationMode = 0 ! save initialization value, since it does not result in any calculation
if (lovl == 4 ) then ! jacobian requested by marc
if (timinc < theDelta .and. theInc == inc .and. lastLovl /= lovl) & ! first after cutback
computationMode = CPFEM_RESTOREJACOBIAN
elseif (lovl == 6) then ! stress requested by marc
computationMode = CPFEM_CALCRESULTS
cp_en = mesh_FEM2DAMASK_elem(m(1))
if (cptim > theTime .or. inc /= theInc) then ! reached "convergence"
terminallyIll = .false.
cycleCounter = -1 ! first calc step increments this to cycle = 0
if (inc == 0) then ! >> start of analysis <<
lastIncConverged = .false.
outdatedByNewInc = .false.
lastLovl = lovl ! pretend that this is NOT the first after a lovl change
write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> start of analysis..! ',m(1),nn
else if (inc - theInc > 1) then ! >> restart of broken analysis <<
lastIncConverged = .false.
outdatedByNewInc = .false.
write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> restart of analysis..! ',m(1),nn
else ! >> just the next inc <<
lastIncConverged = .true.
outdatedByNewInc = .true.
write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> new increment..! ',m(1),nn
endif
else if ( timinc < theDelta ) then ! >> cutBack <<
lastIncConverged = .false.
outdatedByNewInc = .false.
terminallyIll = .false.
cycleCounter = -1 ! first calc step increments this to cycle = 0
write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> cutback detected..! ',m(1),nn
endif ! convergence treatment end
flush(6)
if (usePingPong) then if (lastLovl /= lovl) then
calcMode(nn,cp_en) = .not. calcMode(nn,cp_en) ! ping pong (calc <--> collect) cycleCounter = cycleCounter + 1
if (calcMode(nn,cp_en)) then ! now --- CALC --- !mesh_cellnode = mesh_build_cellnodes() ! update cell node coordinates
computationMode = CPFEM_CALCRESULTS !call mesh_build_ipCoordinates() ! update ip coordinates
if (lastLovl /= lovl) then ! first after ping pong endif
call debug_reset() ! resets debugging if (outdatedByNewInc) then
outdatedFFN1 = .false. computationMode = ior(computationMode,CPFEM_AGERESULTS)
cycleCounter = cycleCounter + 1 outdatedByNewInc = .false.
!mesh_cellnode = mesh_build_cellnodes() ! update cell node coordinates endif
!call mesh_build_ipCoordinates() ! update ip coordinates if (lastIncConverged) then
endif computationMode = ior(computationMode,CPFEM_BACKUPJACOBIAN)
if (outdatedByNewInc) then lastIncConverged = .false.
computationMode = ior(computationMode,CPFEM_AGERESULTS) ! calc and age results endif
outdatedByNewInc = .false. ! reset flag
endif
else ! now --- COLLECT ---
computationMode = CPFEM_COLLECT ! plain collect
if (lastLovl /= lovl .and. & .not. terminallyIll) &
call debug_info() ! first after ping pong reports (meaningful) debugging
if (lastIncConverged) then
computationMode = ior(computationMode,CPFEM_BACKUPJACOBIAN) ! collect and backup Jacobian after convergence
lastIncConverged = .false. ! reset flag
endif
!do node = 1,theMesh%elem%nNodes
!CPnodeID = mesh_element(4+node,cp_en)
!mesh_node(1:ndeg,CPnodeID) = mesh_node0(1:ndeg,CPnodeID) + numerics_unitlength * dispt(1:ndeg,node)
!enddo
endif
else ! --- PLAIN MODE --- theTime = cptim
computationMode = CPFEM_CALCRESULTS ! always calc theDelta = timinc
if (lastLovl /= lovl) then theInc = inc
if (.not. terminallyIll) &
call debug_info() ! first reports (meaningful) debugging
call debug_reset() ! and resets debugging
outdatedFFN1 = .false.
cycleCounter = cycleCounter + 1
!mesh_cellnode = mesh_build_cellnodes() ! update cell node coordinates
!call mesh_build_ipCoordinates() ! update ip coordinates
endif
if (outdatedByNewInc) then
computationMode = ior(computationMode,CPFEM_AGERESULTS)
outdatedByNewInc = .false. ! reset flag
endif
if (lastIncConverged) then
computationMode = ior(computationMode,CPFEM_BACKUPJACOBIAN) ! backup Jacobian after convergence
lastIncConverged = .false. ! reset flag
endif
endif
theTime = cptim ! record current starting time endif
theDelta = timinc ! record current time increment lastLovl = lovl
theInc = inc ! record current increment number
endif call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,m(1),nn,stress,ddsdde)
lastLovl = lovl ! record lovl
call CPFEM_general(computationMode,usePingPong,ffn,ffn1,t(1),timinc,m(1),nn,stress,ddsdde) d = ddsdde(1:ngens,1:ngens)
s = stress(1:ndi+nshear)
g = 0.0_pReal
if(symmetricSolver) d = 0.5_pReal*(d+transpose(d))
d = ddsdde(1:ngens,1:ngens) call omp_set_num_threads(defaultNumThreadsInt) ! reset number of threads to stored default value
s = stress(1:ndi+nshear)
g = 0.0_pReal
if(symmetricSolver) d = 0.5_pReal*(d+transpose(d))
!$ call omp_set_num_threads(defaultNumThreadsInt) ! reset number of threads to stored default value
end subroutine hypela2 end subroutine hypela2

5
src/crystallite.f90
View File

@ -245,7 +245,7 @@ subroutine crystallite_init
enddo enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
if(any(plasticState%nonlocal) .and. .not. usePingPong) call IO_error(601) ! exit if nonlocal but no ping-pong ToDo: Why not check earlier? or in nonlocal? !if(any(plasticState%nonlocal) .and. .not. usePingPong) call IO_error(601)
crystallite_partionedFp0 = crystallite_Fp0 crystallite_partionedFp0 = crystallite_Fp0
crystallite_partionedFi0 = crystallite_Fi0 crystallite_partionedFi0 = crystallite_Fi0
@ -276,9 +276,6 @@ subroutine crystallite_init
write(6,'(a42,1x,i10)') 'max # of constituents/integration point: ', cMax write(6,'(a42,1x,i10)') 'max # of constituents/integration point: ', cMax
flush(6) flush(6)
endif endif
call debug_info
call debug_reset
#endif #endif
end subroutine crystallite_init end subroutine crystallite_init

55
src/debug.f90
View File

@ -49,26 +49,11 @@ module debug
debug_i = 1, & debug_i = 1, &
debug_g = 1 debug_g = 1
integer, dimension(2), public :: &
debug_stressMaxLocation = 0, &
debug_stressMinLocation = 0, &
debug_jacobianMaxLocation = 0, &
debug_jacobianMinLocation = 0
real(pReal), public :: &
debug_stressMax = -huge(1.0_pReal), &
debug_stressMin = huge(1.0_pReal), &
debug_jacobianMax = -huge(1.0_pReal), &
debug_jacobianMin = huge(1.0_pReal)
#ifdef PETSc #ifdef PETSc
character(len=1024), parameter, public :: & character(len=1024), parameter, public :: &
PETSCDEBUG = ' -snes_view -snes_monitor ' PETSCDEBUG = ' -snes_view -snes_monitor '
#endif #endif
public :: debug_init, & public :: debug_init
debug_reset, &
debug_info
contains contains
@ -230,42 +215,4 @@ subroutine debug_init
end subroutine debug_init end subroutine debug_init
!--------------------------------------------------------------------------------------------------
!> @brief resets all debug values
!--------------------------------------------------------------------------------------------------
subroutine debug_reset
debug_stressMaxLocation = 0
debug_stressMinLocation = 0
debug_jacobianMaxLocation = 0
debug_jacobianMinLocation = 0
debug_stressMax = -huge(1.0_pReal)
debug_stressMin = huge(1.0_pReal)
debug_jacobianMax = -huge(1.0_pReal)
debug_jacobianMin = huge(1.0_pReal)
end subroutine debug_reset
!--------------------------------------------------------------------------------------------------
!> @brief writes debug statements to standard out
!--------------------------------------------------------------------------------------------------
subroutine debug_info
!$OMP CRITICAL (write2out)
debugOutputCPFEM: if (iand(debug_level(debug_CPFEM),debug_LEVELBASIC) /= 0 &
.and. any(debug_stressMinLocation /= 0) &
.and. any(debug_stressMaxLocation /= 0) ) then
write(6,'(2/,a,/)') ' Extreme values of returned stress and Jacobian'
write(6,'(a39)') ' value el ip'
write(6,'(a14,1x,e12.3,1x,i8,1x,i4)') ' stress min :', debug_stressMin, debug_stressMinLocation
write(6,'(a14,1x,e12.3,1x,i8,1x,i4)') ' max :', debug_stressMax, debug_stressMaxLocation
write(6,'(a14,1x,e12.3,1x,i8,1x,i4)') ' Jacobian min :', debug_jacobianMin, debug_jacobianMinLocation
write(6,'(a14,1x,e12.3,1x,i8,1x,i4,/)') ' max :', debug_jacobianMax, debug_jacobianMaxLocation
endif debugOutputCPFEM
!$OMP END CRITICAL (write2out)
end subroutine debug_info
end module debug end module debug

3
src/homogenization.f90
View File

@ -249,7 +249,8 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
NiterationHomog = 0 NiterationHomog = 0
cutBackLooping: do while (.not. terminallyIll .and. & cutBackLooping: do while (.not. terminallyIll .and. &
any(subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) > num%subStepMinHomog)) any(subStep(FEsolving_execIP(1):FEsolving_execIP(2),&
FEsolving_execElem(1):FEsolving_execElem(2)) > num%subStepMinHomog))
!$OMP PARALLEL DO PRIVATE(myNgrains) !$OMP PARALLEL DO PRIVATE(myNgrains)
elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2) elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2)

10
src/marc/discretization_marc.f90
View File

@ -45,9 +45,7 @@ contains
!> @brief initializes the mesh by calling all necessary private routines the mesh module !> @brief initializes the mesh by calling all necessary private routines the mesh module
!! Order and routines strongly depend on type of solver !! Order and routines strongly depend on type of solver
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine discretization_marc_init(ip,el) subroutine discretization_marc_init
integer, intent(in) :: el, ip
real(pReal), dimension(:,:), allocatable :: & real(pReal), dimension(:,:), allocatable :: &
node0_elem, & !< node x,y,z coordinates (initially!) node0_elem, & !< node x,y,z coordinates (initially!)
@ -70,7 +68,7 @@ subroutine discretization_marc_init(ip,el)
real(pReal), dimension(:,:,:,:),allocatable :: & real(pReal), dimension(:,:,:,:),allocatable :: &
unscaledNormals unscaledNormals
write(6,'(/,a)') ' <<<+- mesh init -+>>>'; flush(6) write(6,'(/,a)') ' <<<+- discretization_marc init -+>>>'; flush(6)
mesh_unitlength = numerics_unitlength ! set physical extent of a length unit in mesh mesh_unitlength = numerics_unitlength ! set physical extent of a length unit in mesh
@ -83,10 +81,6 @@ subroutine discretization_marc_init(ip,el)
FEsolving_execElem = [1,nElems] FEsolving_execElem = [1,nElems]
FEsolving_execIP = [1,elem%nIPs] FEsolving_execIP = [1,elem%nIPs]
allocate(calcMode(elem%nIPs,nElems),source=.false.) ! pretend to have collected what first call is asking (F = I)
calcMode(ip,mesh_FEM2DAMASK_elem(el)) = .true. ! first ip,el needs to be already pingponged to "calc"
allocate(cellNodeDefinition(elem%nNodes-1)) allocate(cellNodeDefinition(elem%nNodes-1))
allocate(connectivity_cell(elem%NcellNodesPerCell,elem%nIPs,nElems)) allocate(connectivity_cell(elem%NcellNodesPerCell,elem%nIPs,nElems))
call buildCells(connectivity_cell,cellNodeDefinition,& call buildCells(connectivity_cell,cellNodeDefinition,&

5
src/numerics.f90
View File

@ -28,8 +28,6 @@ module numerics
numerics_unitlength = 1.0_pReal, & !< determines the physical length of one computational length unit numerics_unitlength = 1.0_pReal, & !< determines the physical length of one computational length unit
charLength = 1.0_pReal, & !< characteristic length scale for gradient problems charLength = 1.0_pReal, & !< characteristic length scale for gradient problems
residualStiffness = 1.0e-6_pReal !< non-zero residual damage residualStiffness = 1.0e-6_pReal !< non-zero residual damage
logical, protected, public :: &
usePingPong = .true.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! field parameters: ! field parameters:
@ -133,8 +131,6 @@ subroutine numerics_init
defgradTolerance = IO_floatValue(line,chunkPos,2) defgradTolerance = IO_floatValue(line,chunkPos,2)
case ('ijacostiffness') case ('ijacostiffness')
iJacoStiffness = IO_intValue(line,chunkPos,2) iJacoStiffness = IO_intValue(line,chunkPos,2)
case ('usepingpong')
usepingpong = IO_intValue(line,chunkPos,2) > 0
case ('unitlength') case ('unitlength')
numerics_unitlength = IO_floatValue(line,chunkPos,2) numerics_unitlength = IO_floatValue(line,chunkPos,2)
@ -221,7 +217,6 @@ subroutine numerics_init
! writing parameters to output ! writing parameters to output
write(6,'(a24,1x,es8.1)') ' defgradTolerance: ',defgradTolerance write(6,'(a24,1x,es8.1)') ' defgradTolerance: ',defgradTolerance
write(6,'(a24,1x,i8)') ' iJacoStiffness: ',iJacoStiffness write(6,'(a24,1x,i8)') ' iJacoStiffness: ',iJacoStiffness
write(6,'(a24,1x,L8)') ' use ping pong scheme: ',usepingpong
write(6,'(a24,1x,es8.1,/)')' unitlength: ',numerics_unitlength write(6,'(a24,1x,es8.1,/)')' unitlength: ',numerics_unitlength
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------