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

simplified interface of CPFEM_general (removed P and F, made cs and dcs/dE optional) 

commented and cleaned up the marc interface.

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
For marc simulations, run
 ./code/setup/setup_code.sh
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
This commit is contained in:
Martin Diehl 2013-04-09 10:08:00 +00:00
parent 6625500cc7
commit 924d943edc
6 changed files with 294 additions and 303 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

60
code/CPFEM.f90
View File

@ -266,8 +266,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, ffn, ffn1, Temperature, dt, element, IP, cauchyStress,& subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchyStress, jacobian)
& jacobian, pstress, dPdF)
! note: cauchyStress = Cauchy stress cs(6) and jacobian = Consistent tangent dcs/dE ! note: cauchyStress = Cauchy stress cs(6) and jacobian = Consistent tangent dcs/dE
use numerics, only: defgradTolerance, & use numerics, only: defgradTolerance, &
iJacoStiffness iJacoStiffness
@ -343,18 +342,16 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
use DAMASK_interface use DAMASK_interface
implicit none implicit none
integer(pInt), intent(in) :: element, & ! FE element number integer(pInt), intent(in) :: element, & !< FE element number
IP ! FE integration point number IP !< FE integration point number
real(pReal), intent(inout) :: Temperature ! temperature real(pReal), intent(inout) :: Temperature !< temperature
real(pReal), intent(in) :: dt ! time increment real(pReal), intent(in) :: dt !< time increment
real(pReal), dimension (3,3), intent(in) :: ffn, & ! deformation gradient for t=t0 real(pReal), dimension (3,3), intent(in) :: ffn, & !< deformation gradient for t=t0
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), dimension(6), intent(out) :: cauchyStress ! stress vector in Mandel notation real(pReal), dimension(6), intent(out), optional :: cauchyStress !< stress vector in Mandel notation
real(pReal), dimension(6,6), intent(out) :: jacobian ! jacobian in Mandel notation real(pReal), dimension(6,6), intent(out), optional :: 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 real(pReal) J_inverse, & ! inverse of Jacobian
rnd rnd
@ -573,9 +570,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
endif endif
endif endif
! --+>> COLLECTION OF FEM INPUT WITH RETURNING OF RANDOMIZED ODD STRESS AND JACOBIAN <<+-- !--------------------------------------------------------------------------------------------------
! collection of FEM input with returning of randomize odd stress and jacobian
if (iand(mode, CPFEM_BACKUPJACOBIAN) /= 0_pInt) & if (iand(mode, CPFEM_BACKUPJACOBIAN) /= 0_pInt) &
CPFEM_dcsde_knownGood = CPFEM_dcsde CPFEM_dcsde_knownGood = CPFEM_dcsde
if (iand(mode, CPFEM_RESTOREJACOBIAN) /= 0_pInt) & if (iand(mode, CPFEM_RESTOREJACOBIAN) /= 0_pInt) &
@ -592,11 +588,6 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
CPFEM_calc_done = .false. CPFEM_calc_done = .false.
endif endif
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 if (theTime > 0.0_pReal) then
Temperature = materialpoint_Temperature(IP,cp_en) ! homogenized result except for potentially non-isothermal starting condition. Temperature = materialpoint_Temperature(IP,cp_en) ! homogenized result except for potentially non-isothermal starting condition.
endif endif
@ -605,24 +596,24 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
.and. ((debug_e == cp_en .and. debug_i == IP) & .and. ((debug_e == cp_en .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
!$OMP CRITICAL (write2out) !$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,/,12x,6(f10.3,1x)/)') '<< CPFEM >> stress/MPa at el ip ', &
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 cp_en, IP, CPFEM_cs(1:6,IP,cp_en)/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(CPFEM_dcsdE(1:6,1:6,IP,cp_en))/1.0e9_pReal
flush(6) flush(6)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
!--------------------------------------------------------------------------------------------------
!*** warn if stiffness close to zero ! warn if stiffness close to zero
if (all(abs(CPFEM_dcsdE(1:6,1:6,IP,cp_en)) < 1e-10_pReal)) then
if (all(abs(jacobian) < 1e-10_pReal)) then
call IO_warning(601,cp_en,IP) call IO_warning(601,cp_en,IP)
endif endif
!--------------------------------------------------------------------------------------------------
!*** remember extreme values of stress and jacobian ! remember extreme values of stress and jacobian
if (mode < 3) then if (mode < 3) then
cauchyStress33 = math_Mandel6to33(cauchyStress) cauchyStress33 = math_Mandel6to33(CPFEM_cs(1:6,IP,cp_en))
if (maxval(cauchyStress33) > debug_stressMax) then if (maxval(cauchyStress33) > debug_stressMax) then
debug_stressMaxLocation = [cp_en, IP] debug_stressMaxLocation = [cp_en, IP]
debug_stressMax = maxval(cauchyStress33) debug_stressMax = maxval(cauchyStress33)
@ -631,7 +622,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
debug_stressMinLocation = [cp_en, IP] debug_stressMinLocation = [cp_en, IP]
debug_stressMin = minval(cauchyStress33) debug_stressMin = minval(cauchyStress33)
endif endif
jacobian3333 = math_Mandel66to3333(jacobian) jacobian3333 = math_Mandel66to3333(CPFEM_dcsdE(1:6,1:6,IP,cp_en))
if (maxval(jacobian3333) > debug_jacobianMax) then if (maxval(jacobian3333) > debug_jacobianMax) then
debug_jacobianMaxLocation = [cp_en, IP] debug_jacobianMaxLocation = [cp_en, IP]
debug_jacobianMax = maxval(jacobian3333) debug_jacobianMax = maxval(jacobian3333)
@ -642,6 +633,11 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
endif endif
endif endif
!--------------------------------------------------------------------------------------------------
! copy to output if required (FEM solver)
if(present(cauchyStress)) cauchyStress = CPFEM_cs(1:6,IP,cp_en)
if(present(jacobian)) jacobian = CPFEM_dcsdE(1:6,1:6,IP,cp_en)
end subroutine CPFEM_general end subroutine CPFEM_general
end module CPFEM end module CPFEM

4
code/DAMASK_abaqus_exp.f
View File

@ -204,8 +204,6 @@ subroutine vumat(nBlock, nDir, nshr, nStateV, nFieldV, nProps, lAnneal, &
real(pReal), dimension(nblock(1),nstatev), intent(out) :: & real(pReal), dimension(nblock(1),nstatev), intent(out) :: &
stateNew !< state variables at each material point at the end of the increment stateNew !< state variables at each material point at the end of the increment
real(pReal), dimension(3,3) :: pstress ! not used, but needed for call of cpfem_general
real(pReal), dimension(3,3,3,3) :: dPdF ! not used, but needed for call of cpfem_general
real(pReal), dimension(3) :: coordinates real(pReal), dimension(3) :: coordinates
real(pReal), dimension(3,3) :: defgrd0,defgrd1 real(pReal), dimension(3,3) :: defgrd0,defgrd1
real(pReal), dimension(6) :: stress real(pReal), dimension(6) :: stress
@ -286,7 +284,7 @@ subroutine vumat(nBlock, nDir, nshr, nStateV, nFieldV, nProps, lAnneal, &
cp_en = mesh_FEasCP('elem',nBlock(4_pInt+n)) cp_en = mesh_FEasCP('elem',nBlock(4_pInt+n))
mesh_ipCoordinates(1:3,n,cp_en) = numerics_unitlength * coordMp(n,1:3) mesh_ipCoordinates(1:3,n,cp_en) = numerics_unitlength * coordMp(n,1:3)
call CPFEM_general(computationMode,defgrd0,defgrd1,temp,timeInc,cp_en,nBlock(2),stress,ddsdde, pstress, dPdF) call CPFEM_general(computationMode,defgrd0,defgrd1,temp,timeInc,cp_en,nBlock(2),stress,ddsdde)
! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 ! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13
! straight: 11, 22, 33, 12, 23, 13 ! straight: 11, 22, 33, 12, 23, 13

4
code/DAMASK_abaqus_std.f
View File

@ -213,8 +213,6 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,&
DDSDDE !< Jacobian matrix of the constitutive model DDSDDE !< Jacobian matrix of the constitutive model
real(pReal) :: temperature ! temp by Abaqus is intent(in) real(pReal) :: temperature ! temp by Abaqus is intent(in)
real(pReal), dimension (3,3) :: pstress ! not used, but needed for call of cpfem_general
real(pReal), dimension (3,3,3,3) :: dPdF ! not used, but needed for call of cpfem_general
real(pReal), dimension(6) :: stress_h real(pReal), dimension(6) :: stress_h
real(pReal), dimension(6,6) :: ddsdde_h real(pReal), dimension(6,6) :: ddsdde_h
integer(pInt) :: computationMode, i, cp_en integer(pInt) :: computationMode, i, cp_en
@ -313,7 +311,7 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,&
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
call CPFEM_general(computationMode,dfgrd0,dfgrd1,temperature,dtime,noel,npt,stress_h,ddsdde_h, pstress, dPdF) call CPFEM_general(computationMode,dfgrd0,dfgrd1,temperature,dtime,noel,npt,stress_h,ddsdde_h)
! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 ! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13
! straight: 11, 22, 33, 12, 23, 13 ! straight: 11, 22, 33, 12, 23, 13

275
code/DAMASK_marc.f90
View File

@ -30,21 +30,18 @@
!> @details - choose material as hypela2 !> @details - choose material as hypela2
!> @details - set statevariable 2 to index of homogenization !> @details - set statevariable 2 to index of homogenization
!> @details - set statevariable 3 to index of microstructure !> @details - set statevariable 3 to index of microstructure
!> @details - make sure the file "material.config" exists in the working !> @details - make sure the file "material.config" exists in the working directory
!> @details directory !> @details - make sure the file "numerics.config" exists in the working directory
!> @details - make sure the file "numerics.config" exists in the working !> @details - use nonsymmetric option for solver (e.g. direct profile or multifrontal sparse,
!> @details directory !> @details the latter seems to be faster!)
!> @details - use nonsymmetric option for solver (e.g. direct !> @details - in case of ddm (domain decomposition) a SYMMETRIC solver has to be used,
!> @details profile or multifrontal sparse, the latter seems !> @details i.e uncheck "non-symmetric"
!> @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 Marc subroutines used:
!> @details - hypela2 !> @details - hypela2
!> @details - plotv !> @details - plotv
!> @details - quit !> @details - quit
!> @details Marc common blocks included: !> @details Marc common blocks included:
!> @details - concom: lovl, ncycle, inc, incsub !> @details - concom: lovl, inc
!> @details - creeps: timinc !> @details - creeps: timinc
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -74,13 +71,13 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine DAMASK_interface_init subroutine DAMASK_interface_init
implicit none implicit none
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(/,a)') ' <<<+- DAMASK_marc init -+>>>' write(6,'(/,a)') ' <<<+- DAMASK_marc init -+>>>'
write(6,'(a)') ' $Id$' write(6,'(a)') ' $Id$'
#include "compilation_info.f90" #include "compilation_info.f90"
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
end subroutine DAMASK_interface_init end subroutine DAMASK_interface_init
@ -98,7 +95,6 @@ function getSolverWorkingDirectoryName()
inputName='' inputName=''
inquire(5, name=inputName) ! determine inputputfile inquire(5, name=inputName) ! determine inputputfile
getSolverWorkingDirectoryName=inputName(1:scan(inputName,pathSep,back=.true.)) getSolverWorkingDirectoryName=inputName(1:scan(inputName,pathSep,back=.true.))
! write(6,*) 'getSolverWorkingDirectoryName', getSolverWorkingDirectoryName
end function getSolverWorkingDirectoryName end function getSolverWorkingDirectoryName
@ -151,21 +147,11 @@ end module DAMASK_interface
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief This is the MSC.Marc user subroutine for defining material behavior !> @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 (1) F,R,U are only available for continuum and membrane elements (not for
!> @details shells and beams). !> @details shells and beams).
!> @details !> @details
!> @details (2) For total Lagrangian formulation use the -> 'Elasticity,1' card(= !> @details (2) Use the -> 'Plasticity,3' card(=update+finite+large disp+constant d)
!> @details total Lagrange with large disp) in the parameter section of input deck. !> @details in the parameter section of input deck (updated Lagrangian formulation).
!> @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
!> @details The following operation obtains U (stretch tensor) at t=n+1 : !> @details The following operation obtains U (stretch tensor) at t=n+1 :
!> @details !> @details
@ -178,47 +164,18 @@ end module DAMASK_interface
!> @details enddo !> @details enddo
!> @details enddo !> @details enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine hypela2(& subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
d,& !< stress strain law to be formed dispt,coord,ffn,frotn,strechn,eigvn,ffn1,frotn1, &
g,& !< change in stress due to temperature effects strechn1,eigvn1,ncrd,itel,ndeg,ndm,nnode, &
e,& !< total elastic strain jtype,lclass,ifr,ifu)
de,& !< increment of strain use prec, only: &
s,& !< stress - should be updated by user pReal, &
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, &
pInt pInt
use numerics, only: numerics_unitlength use numerics, only: &
use FEsolving, only: cycleCounter, & !$ DAMASK_NumThreadsInt, &
numerics_unitlength
use FEsolving, only: &
cycleCounter, &
theInc, & theInc, &
calcMode, & calcMode, &
lastMode, & lastMode, &
@ -229,10 +186,17 @@ subroutine hypela2(&
outdatedFFN1, & outdatedFFN1, &
terminallyIll, & terminallyIll, &
symmetricSolver symmetricSolver
use math, only: invnrmMandel use math, only: &
use debug, only: debug_info, & math_transpose33,&
invnrmMandel
use debug, only: &
debug_level, &
debug_LEVELBASIC, &
debug_MARC, &
debug_info, &
debug_reset debug_reset
use mesh, only: mesh_FEasCP, & use mesh, only: &
mesh_FEasCP, &
mesh_element, & mesh_element, &
mesh_node0, & mesh_node0, &
mesh_node, & mesh_node, &
@ -250,52 +214,89 @@ subroutine hypela2(&
CPFEM_RESTOREJACOBIAN, & CPFEM_RESTOREJACOBIAN, &
CPFEM_BACKUPJACOBIAN CPFEM_BACKUPJACOBIAN
!$ use numerics, only: DAMASK_NumThreadsInt ! number of threads set by DAMASK_NUM_THREADS
implicit none implicit none
!$ include "omp_lib.h" ! the openMP function library !$ include "omp_lib.h" ! the openMP function library
! ** Start of generated type statements ** integer(pInt), intent(in) :: & ! according to MSC.Marc 2012 Manual D
real(pReal) coord, d, de, disp, dispt, dt, e, eigvn, eigvn1, ffn, ffn1 ngens, & !< size of stress-strain law
real(pReal) frotn, frotn1, g nn, & !< integration point number
integer(pInt) ifr, ifu, itel, jtype, kcus, lclass, matus, n, ncrd, ndeg ndi, & !< number of direct components
integer(pInt) ndi, ndm, ngens, nn, nnode, nshear nshear, & !< number of shear components
real(pReal) s, strechn, strechn1, t ncrd, & !< number of coordinates
logical :: cutBack itel, & !< dimension of F and R, either 2 or 3
! ** End of generated type statements ** ndeg, & !< number of degrees of freedom
ndm, & !< not specified in MSC.Marc 2012 Manual D
dimension e(*),de(*),t(*),dt(*),g(*),d(ngens,*),s(*), n(2),coord(ncrd,*),disp(ndeg,*),matus(2),dispt(ndeg,*),ffn(itel,*),& nnode, & !< number of nodes per element
frotn(itel,*),strechn(itel),eigvn(itel,*),ffn1(itel,*),frotn1(itel,*),strechn1(itel),eigvn1(itel,*),kcus(2), lclass(2) jtype, & !< element type
ifr, & !< set to 1 if R has been calculated
ifu !< set to 1 if stretch has been calculated
integer(pInt), dimension(2), intent(in) :: & ! according to MSC.Marc 2012 Manual D
m, & !< (1) user element number, (2) internal element number
matus, & !< (1) user material identification number, (2) internal material identification number
kcus, & !< (1) layer number, (2) internal layer number
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(*)
e, & !< total elastic strain
de, & !< increment of strain
dt !< increment of state variables
real(pReal), dimension(itel), intent(in) :: & ! according to MSC.Marc 2012 Manual D
strechn, & !< square of principal stretch ratios, lambda(i) at t=n
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)
ffn, & !< deformation gradient at t=n
ffn1 !< deformation gradient at t=n+1
real(pReal), dimension(itel,*), intent(in) :: & ! according to MSC.Marc 2012 Manual D
frotn, & !< rotation tensor at t=n
eigvn, & !< i principal direction components for j eigenvalues at t=n
frotn1, & !< rotation tensor 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
disp, & !< incremental displacements
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
coord !< coordinates
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)
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
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,*)
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)
! Beware of changes in newer Marc versions ! Beware of changes in newer Marc versions
include "include/concom%%MARCVERSION%%" ! concom is needed for inc, lovl
include "include/concom%%MARCVERSION%%" ! concom is needed for inc, subinc, ncycle, lovl
include "include/creeps%%MARCVERSION%%" ! creeps is needed for timinc (time increment) include "include/creeps%%MARCVERSION%%" ! creeps is needed for timinc (time increment)
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(pInt) :: computationMode, i, cp_en, node, FEnodeID
!$ integer(pInt) :: defaultNumThreadsInt !< default value set by Marc
real(pReal), dimension (3,3) :: pstress ! dummy argument for call of cpfem_general (used by mpie_spectral) !$ defaultNumThreadsInt = omp_get_num_threads() ! remember number of threads set by Marc
real(pReal), dimension (3,3,3,3) :: dPdF ! dummy argument for call of cpfem_general (used by mpie_spectral) if (iand(debug_level(debug_MARC),debug_LEVELBASIC) /= 0_pInt) then
write(6,'(a,/,i8,i8,i2)') ' MSC.MARC information on shape of element(2), IP:', m, nn
write(6,'(a,2(1i))'), ' Jacobian: ', ngens,ngens
write(6,'(a,1i)'), ' Direct stress: ', ndi
write(6,'(a,1i)'), ' Shear stress: ', nshear
write(6,'(a,1i)'), ' DoF: ', ndeg
write(6,'(a,1i)'), ' Coordinates: ', ncrd
write(6,'(a,1i)'), ' Nodes: ', nnode
write(6,'(a,1i)'), ' Deformation gradient: ', itel
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n:', &
math_transpose33(ffn)
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n+1:', &
math_transpose33(ffn1)
endif
if (.not. CPFEM_init_done) call CPFEM_initAll(t(1),m(1),nn)
integer(pInt) computationMode, i, cp_en !$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS
integer(pInt) node, FEnodeID
! OpenMP variable
!$ integer(pInt) defaultNumThreadsInt ! default value set by Marc
!$ defaultNumThreadsInt = omp_get_num_threads() ! remember number of threads set by Marc
if (.not. CPFEM_init_done) call CPFEM_initAll(t(1),n(1),nn)
!$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS
if (lovl == 4 ) then if (lovl == 4 ) then
if(timinc < theDelta .and. theInc == inc ) & ! first after cutback if(timinc < theDelta .and. theInc == inc ) computationMode = CPFEM_RESTOREJACOBIAN ! first after cutback
computationMode = CPFEM_RESTOREJACOBIAN
else ! stress requested (lovl == 6) else ! stress requested (lovl == 6)
cp_en = mesh_FEasCP('elem',n(1)) cp_en = mesh_FEasCP('elem',m(1))
if (cptim > theTime .or. inc /= theInc) then ! reached "convergence" if (cptim > theTime .or. inc /= theInc) then ! reached "convergence"
terminallyIll = .false. terminallyIll = .false.
cycleCounter = -1 ! first calc step increments this to cycle = 0 cycleCounter = -1 ! first calc step increments this to cycle = 0
@ -305,7 +306,7 @@ subroutine hypela2(&
lastMode = .false. ! pretend last step was collection lastMode = .false. ! pretend last step was collection
calcMode = .false. ! pretend last step was collection calcMode = .false. ! pretend last step was collection
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> start of analysis..! ',n(1),nn write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> start of analysis..! ',m(1),nn
flush(6) flush(6)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
else if (inc - theInc > 1) then ! >> restart of broken analysis << else if (inc - theInc > 1) then ! >> restart of broken analysis <<
@ -314,7 +315,7 @@ subroutine hypela2(&
lastMode = .true. ! pretend last step was calculation lastMode = .true. ! pretend last step was calculation
calcMode = .true. ! pretend last step was calculation calcMode = .true. ! pretend last step was calculation
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> restart of analysis..! ',n(1),nn write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> restart of analysis..! ',m(1),nn
flush(6) flush(6)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
else ! >> just the next inc << else ! >> just the next inc <<
@ -323,7 +324,7 @@ subroutine hypela2(&
lastMode = .true. ! assure last step was calculation lastMode = .true. ! assure last step was calculation
calcMode = .true. ! assure last step was calculation calcMode = .true. ! assure last step was calculation
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> new increment..! ',n(1),nn write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> new increment..! ',m(1),nn
flush(6) flush(6)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
@ -332,7 +333,7 @@ subroutine hypela2(&
cycleCounter = -1 ! first calc step increments this to cycle = 0 cycleCounter = -1 ! first calc step increments this to cycle = 0
calcMode = .true. ! pretend last step was calculation calcMode = .true. ! pretend last step was calculation
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> cutback detected..! ',n(1),nn write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> cutback detected..! ',m(1),nn
flush(6) flush(6)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif ! convergence treatment end endif ! convergence treatment end
@ -354,10 +355,7 @@ subroutine hypela2(&
computationMode = CPFEM_CALCRESULTS computationMode = CPFEM_CALCRESULTS
endif endif
else ! now --- COLLECT --- else ! now --- COLLECT ---
if ( lastMode /= calcMode(nn,cp_en) .and. & if ( lastMode /= calcMode(nn,cp_en) .and. & .not. terminallyIll) call debug_info() ! first after ping pong reports (meaningful) debugging
.not. terminallyIll ) then
call debug_info() ! first after ping pong reports (meaningful) debugging
endif
if ( lastIncConverged ) then if ( lastIncConverged ) then
computationMode = ior(CPFEM_COLLECT,CPFEM_BACKUPJACOBIAN) ! collect and backup Jacobian after convergence computationMode = ior(CPFEM_COLLECT,CPFEM_BACKUPJACOBIAN) ! collect and backup Jacobian after convergence
lastIncConverged = .false. ! reset flag lastIncConverged = .false. ! reset flag
@ -376,17 +374,18 @@ subroutine hypela2(&
lastMode = calcMode(nn,cp_en) ! record calculationMode lastMode = calcMode(nn,cp_en) ! record calculationMode
endif endif
call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,stress,ddsdde, pstress, dPdF) call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,m(1),nn,stress,ddsdde)
! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 ! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13
! Marc: 11, 22, 33, 12, 23, 13 ! Marc: 11, 22, 33, 12, 23, 13
! Marc: 11, 22, 33, 12 ! Marc: 11, 22, 33, 12
forall(i=1:ngens) d(1:ngens,i) = invnrmMandel(i)*ddsdde(1:ngens,i)*invnrmMandel(1:ngens) forall(i=1:ngens) d(1:ngens,i) = invnrmMandel(i)*ddsdde(1:ngens,i)*invnrmMandel(1:ngens)
s(1:ngens) = stress(1:ngens)*invnrmMandel(1:ngens) s(1:ndi+nshear) = stress(1:ndi+nshear)*invnrmMandel(1:ndi+nshear)
if(symmetricSolver) d(1:ngens,1:ngens) = 0.5_pReal*(d(1:ngens,1:ngens)+transpose(d(1:ngens,1:ngens))) 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 !$ call omp_set_num_threads(defaultNumThreadsInt) ! reset number of threads to stored default value
end subroutine hypela2 end subroutine hypela2
@ -395,33 +394,37 @@ end subroutine hypela2
!> @brief sets user defined output variables for Marc !> @brief sets user defined output variables for Marc
!> @details select a variable contour plotting (user subroutine). !> @details select a variable contour plotting (user subroutine).
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine plotv(& subroutine plotv(v,s,sp,etot,eplas,ecreep,t,m,nn,layer,ndi,nshear,jpltcd)
v,& !< variable use prec, only: &
s,& !< stress array pReal, &
sp,& !< stresses in preferred direction pInt
etot,& !< total strain (generalized) use mesh, only: &
eplas,& !< total plastic strain mesh_FEasCP
ecreep,& !< total creep strain use IO, only: &
t,& !< current temperature IO_error
m,& !< element number use homogenization, only: &
nn,& !< integration point number materialpoint_results,&
layer,& !< layer number materialpoint_sizeResults
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
use IO, only: IO_error
use homogenization, only: materialpoint_results,materialpoint_sizeResults
implicit none implicit none
real(pReal) s(*),etot(*),eplas(*),ecreep(*),sp(*) integer(pInt), intent(in) :: &
real(pReal) v, t(*) m, & !< element number
integer(pInt) m, nn, layer, ndi, nshear, jpltcd nn, & !< integration point number
layer, & !< layer number
ndi, & !< number of direct stress components
nshear, & !< number of shear stress components
jpltcd !< user variable index
real(pReal), dimension(*), intent(in) :: &
s, & !< stress array
sp, & !< stresses in preferred direction
etot, & !< total strain (generalized)
eplas, & !< total plastic strain
ecreep, & !< total creep strain
t !< current temperature
real(pReal), intent(out) :: &
v !< variable
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)) v = materialpoint_results(jpltcd,nn,mesh_FEasCP('elem', m))
end subroutine plotv end subroutine plotv

40
code/DAMASK_spectral_utilities.f90
View File

@ -764,9 +764,6 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,&
integer(pInt) :: & integer(pInt) :: &
calcMode, & !< CPFEM mode for calculation calcMode, & !< CPFEM mode for calculation
collectMode !< CPFEM mode for collection collectMode !< CPFEM mode for collection
real(pReal), dimension(3,3,3,3) :: dPdF !< d P / d F
real(pReal), dimension(6) :: sigma !< cauchy stress in mandel notation
real(pReal), dimension(6,6) :: dsde !< d sigma / d Epsilon
real(pReal), dimension(3,3,3,3) :: max_dPdF, min_dPdF real(pReal), dimension(3,3,3,3) :: max_dPdF, min_dPdF
real(pReal) :: max_dPdF_norm, min_dPdF_norm, defgradDetMin, defgradDetMax, defgradDet real(pReal) :: max_dPdF_norm, min_dPdF_norm, defgradDetMin, defgradDetMax, defgradDet
integer(pInt) :: i,j,k integer(pInt) :: i,j,k
@ -784,24 +781,8 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,&
calcMode = iand(calcMode, not(CPFEM_AGERESULTS)) calcMode = iand(calcMode, not(CPFEM_AGERESULTS))
endif endif
!--------------------------------------------------------------------------------------------------
! calculate bounds of det(F) and report
if(debugGeneral) then
defgradDetMax = -huge(1.0_pReal)
defgradDetMin = +huge(1.0_pReal)
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
defgradDet = math_det33(F(1:3,1:3,i,j,k))
defgradDetMax = max(defgradDetMax,defgradDet)
defgradDetMin = min(defgradDetMin,defgradDet)
enddo; enddo; enddo
write(6,'(a,1x,es11.4)') ' max determinant of deformation =', defgradDetMax
write(6,'(a,1x,es11.4)') ' min determinant of deformation =', defgradDetMin
flush(6)
endif
call CPFEM_general(collectMode,F_lastInc(1:3,1:3,1,1,1),F(1:3,1:3,1,1,1), & ! collect mode handles Jacobian backup / restoration call CPFEM_general(collectMode,F_lastInc(1:3,1:3,1,1,1),F(1:3,1:3,1,1,1), & ! collect mode handles Jacobian backup / restoration
temperature,timeinc,1_pInt,1_pInt,sigma,dsde,P(1:3,1:3,1,1,1),dPdF) temperature,timeinc,1_pInt,1_pInt)
materialpoint_F0 = reshape(F_lastInc, [3,3,1,mesh_NcpElems]) materialpoint_F0 = reshape(F_lastInc, [3,3,1,mesh_NcpElems])
materialpoint_F = reshape(F, [3,3,1,mesh_NcpElems]) materialpoint_F = reshape(F, [3,3,1,mesh_NcpElems])
@ -809,8 +790,23 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,&
call debug_reset() call debug_reset()
!--------------------------------------------------------------------------------------------------
! calculate bounds of det(F) and report
if(debugGeneral) then
defgradDetMax = -huge(1.0_pReal)
defgradDetMin = +huge(1.0_pReal)
do j = 1_pInt, mesh_NcpElems
defgradDet = math_det33(materialpoint_F(1:3,1:3,1,j))
defgradDetMax = max(defgradDetMax,defgradDet)
defgradDetMin = min(defgradDetMin,defgradDet)
end do
write(6,'(a,1x,es11.4)') ' max determinant of deformation =', defgradDetMax
write(6,'(a,1x,es11.4)') ' min determinant of deformation =', defgradDetMin
flush(6)
endif
call CPFEM_general(calcMode,F_lastInc(1:3,1:3,1,1,1), F(1:3,1:3,1,1,1), & ! first call calculates everything call CPFEM_general(calcMode,F_lastInc(1:3,1:3,1,1,1), F(1:3,1:3,1,1,1), & ! first call calculates everything
temperature,timeinc,1_pInt,1_pInt,sigma,dsde,P(1:3,1:3,1,1,1),dPdF) temperature,timeinc,1_pInt,1_pInt)
max_dPdF = 0.0_pReal max_dPdF = 0.0_pReal
max_dPdF_norm = 0.0_pReal max_dPdF_norm = 0.0_pReal
@ -825,7 +821,7 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,&
min_dPdF = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k) min_dPdF = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k)
min_dPdF_norm = sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k)**2.0_pReal) min_dPdF_norm = sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k)**2.0_pReal)
endif endif
enddo end do
P = reshape(materialpoint_P, [3,3,res(1),res(2),res(3)]) P = reshape(materialpoint_P, [3,3,res(1),res(2),res(3)])
C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) * wgt C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) * wgt

38
code/debug.f90
View File

@ -214,7 +214,7 @@ subroutine debug_init
what = debug_CPFEM what = debug_CPFEM
case ('spectral') case ('spectral')
what = debug_SPECTRAL what = debug_SPECTRAL
case ('MARC') case ('marc')
what = debug_MARC what = debug_MARC
case ('abaqus') case ('abaqus')
what = debug_ABAQUS what = debug_ABAQUS
@ -268,42 +268,42 @@ subroutine debug_init
do i = 1_pInt, debug_MAXNTYPE do i = 1_pInt, debug_MAXNTYPE
select case(i) select case(i)
case (debug_DEBUG) case (debug_DEBUG)
tag = 'Debug' tag = ' Debug'
case (debug_MATH) case (debug_MATH)
tag = 'Math' tag = ' Math'
case (debug_FESOLVING) case (debug_FESOLVING)
tag = 'FEsolving' tag = ' FEsolving'
case (debug_MESH) case (debug_MESH)
tag = 'Mesh' tag = ' Mesh'
case (debug_MATERIAL) case (debug_MATERIAL)
tag = 'Material' tag = ' Material'
case (debug_LATTICE) case (debug_LATTICE)
tag = 'Lattice' tag = ' Lattice'
case (debug_CONSTITUTIVE) case (debug_CONSTITUTIVE)
tag = 'Constitutive' tag = ' Constitutive'
case (debug_CRYSTALLITE) case (debug_CRYSTALLITE)
tag = 'Crystallite' tag = ' Crystallite'
case (debug_HOMOGENIZATION) case (debug_HOMOGENIZATION)
tag = 'Homogenizaiton' tag = ' Homogenizaiton'
case (debug_CPFEM) case (debug_CPFEM)
tag = 'CPFEM' tag = ' CPFEM'
case (debug_SPECTRAL) case (debug_SPECTRAL)
tag = 'Spectral solver' tag = ' Spectral solver'
case (debug_MARC) case (debug_MARC)
tag = 'MSC.MARC FEM solver' tag = ' MSC.MARC FEM solver'
case (debug_ABAQUS) case (debug_ABAQUS)
tag = 'ABAQUS FEM solver' tag = ' ABAQUS FEM solver'
end select end select
if(debug_level(i) /= 0) then if(debug_level(i) /= 0) then
write(6,'(a,a)') tag,' debugging:' write(6,'(3a)') ' debug level for ', trim(tag), ':'
if(iand(debug_level(i),debug_LEVELBASIC) /= 0) write(6,'(a)') ' basic' if(iand(debug_level(i),debug_LEVELBASIC) /= 0) write(6,'(a)') ' basic'
if(iand(debug_level(i),debug_LEVELEXTENSIVE) /= 0) write(6,'(a)') ' extensive' if(iand(debug_level(i),debug_LEVELEXTENSIVE) /= 0) write(6,'(a)') ' extensive'
if(iand(debug_level(i),debug_LEVELSELECTIVE) /= 0) then if(iand(debug_level(i),debug_LEVELSELECTIVE) /= 0) then
write(6,'(a)') 'selective on:' write(6,'(a)') ' selective on:'
write(6,'(a24,1x,i8)') 'element: ',debug_e write(6,'(a24,1x,i8)') ' element: ',debug_e
write(6,'(a24,1x,i8)') 'ip: ',debug_i write(6,'(a24,1x,i8)') ' ip: ',debug_i
write(6,'(a24,1x,i8)') 'grain: ',debug_g write(6,'(a24,1x,i8)') ' grain: ',debug_g
endif endif
if(iand(debug_level(i),debug_SPECTRALRESTART) /= 0) write(6,'(a)') ' restart' if(iand(debug_level(i),debug_SPECTRALRESTART) /= 0) write(6,'(a)') ' restart'
if(iand(debug_level(i),debug_SPECTRALFFTW) /= 0) write(6,'(a)') ' FFTW' if(iand(debug_level(i),debug_SPECTRALFFTW) /= 0) write(6,'(a)') ' FFTW'