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

openmp parallelization working again (at least for j2 and nonlocal constitutive model). 

In order to keep it like that, please follow these simple rules:

DON'T use implicit array subscripts:
example:    real, dimension(3,3) :: A,B
                  A(:,2) = B(:,1)               <--- DON'T USE
                  A(1:3,2) = B(1:3,1)       <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to  prevent memory leaks.

Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)

Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
This commit is contained in:
Christoph Kords 2011-03-17 10:46:17 +00:00
parent 6ac2b4cf88
commit 235266b169
19 changed files with 1427 additions and 1084 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

72
code/CPFEM.f90
View File

@ -264,6 +264,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
materialpoint_P, &
materialpoint_dPdF, &
materialpoint_results, &
materialpoint_sizeResults, &
materialpoint_Temperature, &
materialpoint_stressAndItsTangent, &
materialpoint_postResults
@ -320,7 +321,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
cp_en = mesh_FEasCP('elem',element)
if (selectiveDebugger .and. cp_en == debug_e .and. IP == debug_i) then
if (cp_en == debug_e .and. IP == debug_i) then
!$OMP CRITICAL (write2out)
write(6,*)
write(6,'(a)') '#######################################################'
@ -351,7 +352,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
j = 1:mesh_maxNips, &
k = 1:mesh_NcpElems ) &
constitutive_state0(i,j,k)%p = constitutive_state(i,j,k)%p ! microstructure of crystallites
if (selectiveDebugger .and. cp_en == debug_e .and. IP == debug_i) then
if (cp_en == debug_e .and. IP == debug_i) then
!$OMP CRITICAL (write2out)
write(6,'(a,x,i8,x,i2,/,4(3(e20.8,x),/))') '<< cpfem >> AGED state of grain 1, element ip',&
cp_en,IP, constitutive_state(1,IP,cp_en)%p
@ -425,25 +426,25 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
endif
if (mode == 8 .or. mode == 9) then ! Abaqus explicit skips collect
materialpoint_Temperature(IP,cp_en) = Temperature
materialpoint_F0(:,:,IP,cp_en) = ffn
materialpoint_F(:,:,IP,cp_en) = ffn1
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(:,:,IP,cp_en)) > defgradTolerance)) then
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
!$OMP CRITICAL (write2out)
write(6,'(a,x,i5,x,i2)') '<< cpfem >> OUTDATED at element ip',cp_en,IP
write(6,'(a,/,3(3(f10.6,x),/))') ' FFN1 old:',math_transpose3x3(materialpoint_F(:,:,IP,cp_en))
write(6,'(a,/,3(3(f10.6,x),/))') ' FFN1 now:',math_transpose3x3(ffn1(:,:))
write(6,'(a,/,3(3(f10.6,x),/))') ' FFN1 old:',math_transpose3x3(materialpoint_F(1:3,1:3,IP,cp_en))
write(6,'(a,/,3(3(f10.6,x),/))') ' FFN1 now:',math_transpose3x3(ffn1)
!$OMP END CRITICAL (write2out)
outdatedFFN1 = .true.
endif
call random_number(rnd)
rnd = 2.0_pReal * rnd - 1.0_pReal
CPFEM_cs(:,IP,cp_en) = rnd*CPFEM_odd_stress
CPFEM_dcsde(:,:,IP,cp_en) = CPFEM_odd_jacobian*math_identity2nd(6)
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
@ -467,10 +468,10 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! loop over all parallely processed elements
if (microstructure_elemhomo(mesh_element(4,e))) then ! dealing with homogeneous element?
forall (i = 2:FE_Nips(mesh_element(2,e))) ! copy results of first IP to all others
materialpoint_P(:,:,i,e) = materialpoint_P(:,:,1,e)
materialpoint_F(:,:,i,e) = materialpoint_F(:,:,1,e)
materialpoint_dPdF(:,:,:,:,i,e) = materialpoint_dPdF(:,:,:,:,1,e)
materialpoint_results(:,i,e) = materialpoint_results(:,1,e)
materialpoint_P(1:3,1:3,i,e) = materialpoint_P(1:3,1:3,1,e)
materialpoint_F(1:3,1:3,i,e) = materialpoint_F(1:3,1:3,1,e)
materialpoint_dPdF(1:3,1:3,1:3,1:3,i,e) = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e)
materialpoint_results(1:materialpoint_sizeResults,i,e) = materialpoint_results(1:materialpoint_sizeResults,1,e)
end forall
endif
enddo
@ -480,13 +481,13 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
if ( terminallyIll ) then
call random_number(rnd)
rnd = 2.0_pReal * rnd - 1.0_pReal
CPFEM_cs(:,IP,cp_en) = rnd*CPFEM_odd_stress
CPFEM_dcsde(:,:,IP,cp_en) = CPFEM_odd_jacobian*math_identity2nd(6)
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
! translate from P to CS
Kirchhoff = math_mul33x33(materialpoint_P(:,:,IP, cp_en),transpose(materialpoint_F(:,:,IP, cp_en)))
J_inverse = 1.0_pReal/math_det3x3(materialpoint_F(:,:,IP, cp_en))
CPFEM_cs(:,IP,cp_en) = math_Mandel33to6(J_inverse*Kirchhoff)
Kirchhoff = math_mul33x33(materialpoint_P(1:3,1:3,IP, cp_en), math_transpose3x3(materialpoint_F(1:3,1:3,IP,cp_en)))
J_inverse = 1.0_pReal / math_det3x3(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
@ -495,14 +496,14 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
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))
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))
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(:,:,IP,cp_en) = math_Mandel3333to66(J_inverse*H_sym)
CPFEM_dcsde(1:6,1:6,IP,cp_en) = math_Mandel3333to66(J_inverse * H_sym)
endif
endif
@ -515,11 +516,11 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
end if
call random_number(rnd)
rnd = 2.0_pReal * rnd - 1.0_pReal
materialpoint_Temperature(IP,cp_en) = Temperature
materialpoint_F0(:,:,IP,cp_en) = ffn
materialpoint_F(:,:,IP,cp_en) = ffn1
CPFEM_cs(:,IP,cp_en) = rnd*CPFEM_odd_stress
CPFEM_dcsde(:,:,IP,cp_en) = CPFEM_odd_jacobian*math_identity2nd(6)
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) <<+--
@ -532,29 +533,28 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
end select
! return the local stress and the jacobian from storage
cauchyStress(:) = CPFEM_cs(:,IP,cp_en)
jacobian(:,:) = CPFEM_dcsdE(:,:,IP,cp_en)
cauchyStress = CPFEM_cs(1:6,IP,cp_en)
jacobian = CPFEM_dcsdE(1:6,1:6,IP,cp_en)
! copy P and dPdF to the output variables
pstress(:,:) = materialpoint_P(:,:,IP,cp_en)
dPdF(:,:,:,:) = materialpoint_dPdF(:,:,:,:,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)
! warning for zero stiffness
if (all(abs(jacobian) < 1e-10_pReal)) then
call IO_warning(601,cp_en,IP)
endif
if (selectiveDebugger .and. cp_en == debug_e .and. IP == debug_i .and. mode < 6) then
if (cp_en == debug_e .and. IP == debug_i .and. mode < 6) then
!$OMP CRITICAL (write2out)
write(6,'(a,x,i2,x,a,x,i4,/,6(f10.3,x)/)') 'stress/MPa at ip', IP, 'el', cp_en, cauchyStress/1e6
write(6,'(a,x,i2,x,a,x,i4,/,6(6(f10.3,x)/))') 'jacobian/GPa at ip', IP, 'el', cp_en, transpose(jacobian(:,:))/1e9
write(6,'(a,x,i2,x,a,x,i4,/,6(6(f10.3,x)/))') 'jacobian/GPa at ip', IP, 'el', cp_en, transpose(jacobian)/1e9
call flush(6)
!$OMP END CRITICAL (write2out)
endif
! return temperature
if (theTime > 0.0_pReal) Temperature = materialpoint_Temperature(IP,cp_en) ! homogenized result except for potentially non-isothermal starting condition.
return
end subroutine

7
code/IO.f90
View File

@ -66,8 +66,10 @@ recursive function IO_abaqus_assembleInputFile(unit1,unit2) result(createSuccess
fname = trim(getSolverWorkingDirectoryName())//trim(line(9+scan(line(9:),'='):))
inquire(file=fname, exist=fexist)
if (.not.(fexist)) then
write(6,*)'ERROR: file does not exist error in IO_abaqus_assembleInputFile'
write(6,*)'filename: ', trim(fname)
!$OMP CRITICAL (write2out)
write(6,*)'ERROR: file does not exist error in IO_abaqus_assembleInputFile'
write(6,*)'filename: ', trim(fname)
!$OMP END CRITICAL (write2out)
createSuccess = .false.
return
endif
@ -1361,7 +1363,6 @@ endfunction
endif
endif
write(6,'(a38)') '+------------------------------------+'
call flush(6)
call quit(9000+ID)
!$OMP END CRITICAL (write2out)

336
code/constitutive.f90
View File

@ -402,123 +402,135 @@ return
endfunction
subroutine constitutive_microstructure(Temperature,Tstar_v,Fe,Fp,ipc,ip,el)
!*********************************************************************
!* This function calculates from state needed variables *
!* INPUT: *
!* - state : state variables *
!* - Tp : temperature *
!* - ipc : component-ID of current integration point *
!* - ip : current integration point *
!* - el : current element *
!*********************************************************************
use prec, only: pReal,pInt
use material, only: phase_constitution, &
subroutine constitutive_microstructure(Temperature,Tstar_v,Fe,Fp,ipc,ip,el)
use prec, only: pReal,pInt
use material, only: phase_constitution, &
material_phase, &
homogenization_maxNgrains
use mesh, only: mesh_NcpElems, &
use mesh, only: mesh_NcpElems, &
mesh_maxNips, &
mesh_maxNipNeighbors
use constitutive_j2
use constitutive_phenopowerlaw
use constitutive_titanmod
use constitutive_dislotwin
use constitutive_nonlocal
implicit none
use constitutive_j2, only: constitutive_j2_label, &
constitutive_j2_microstructure
use constitutive_phenopowerlaw, only: constitutive_phenopowerlaw_label, &
constitutive_phenopowerlaw_microstructure
use constitutive_titanmod, only: constitutive_titanmod_label, &
constitutive_titanmod_microstructure
use constitutive_dislotwin, only: constitutive_dislotwin_label, &
constitutive_dislotwin_microstructure
use constitutive_nonlocal, only: constitutive_nonlocal_label, &
constitutive_nonlocal_microstructure
implicit none
!* Definition of variables
integer(pInt), intent(in) :: ipc,ip,el
real(pReal), intent(in) :: Temperature
real(pReal), dimension(6) :: Tstar_v
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: Fe, Fp
!*** input variables ***!
integer(pInt), intent(in):: ipc, & ! component-ID of current integration point
ip, & ! current integration point
el ! current element
real(pReal), intent(in) :: Temperature
real(pReal), intent(in), dimension(6) :: Tstar_v ! 2nd Piola-Kirchhoff stress
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
Fe, & ! elastic deformation gradient
Fp ! plastic deformation gradient
select case (phase_constitution(material_phase(ipc,ip,el)))
!*** output variables ***!
!*** local variables ***!
select case (phase_constitution(material_phase(ipc,ip,el)))
case (constitutive_j2_label)
call constitutive_j2_microstructure(Temperature,constitutive_state,ipc,ip,el)
case (constitutive_j2_label)
call constitutive_j2_microstructure(Temperature,constitutive_state,ipc,ip,el)
case (constitutive_phenopowerlaw_label)
call constitutive_phenopowerlaw_microstructure(Temperature,constitutive_state,ipc,ip,el)
case (constitutive_titanmod_label)
call constitutive_titanmod_microstructure(Temperature,constitutive_state,ipc,ip,el)
case (constitutive_phenopowerlaw_label)
call constitutive_phenopowerlaw_microstructure(Temperature,constitutive_state,ipc,ip,el)
case (constitutive_titanmod_label)
call constitutive_titanmod_microstructure(Temperature,constitutive_state,ipc,ip,el)
case (constitutive_dislotwin_label)
call constitutive_dislotwin_microstructure(Temperature,constitutive_state,ipc,ip,el)
case (constitutive_nonlocal_label)
call constitutive_nonlocal_microstructure(constitutive_state, Temperature, Tstar_v, Fe, Fp, ipc, ip, el)
case (constitutive_dislotwin_label)
call constitutive_dislotwin_microstructure(Temperature,constitutive_state,ipc,ip,el)
case (constitutive_nonlocal_label)
call constitutive_nonlocal_microstructure(constitutive_state, Temperature, Tstar_v, Fe, Fp, ipc, ip, el)
end select
end select
endsubroutine
subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, Temperature, ipc, ip, el)
!*********************************************************************
!* This subroutine contains the constitutive equation for *
!* calculating the velocity gradient *
!* INPUT: *
!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) *
!* - ipc : component-ID of current integration point *
!* - ip : current integration point *
!* - el : current element *
!* OUTPUT: *
!* - Lp : plastic velocity gradient *
!* - dLp_dTstar : derivative of Lp (4th-order tensor) *
!*********************************************************************
use prec, only: pReal,pInt
use material, only: phase_constitution,material_phase
use constitutive_j2
use constitutive_phenopowerlaw
use constitutive_titanmod
use constitutive_dislotwin
use constitutive_nonlocal
implicit none
subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, Temperature, ipc, ip, el)
!* Definition of variables
integer(pInt) ipc,ip,el
real(pReal) Temperature
real(pReal), dimension(6) :: Tstar_v
real(pReal), dimension(3,3) :: Lp
real(pReal), dimension(9,9) :: dLp_dTstar
use prec, only: pReal,pInt
use material, only: phase_constitution, &
material_phase
use constitutive_j2, only: constitutive_j2_label, &
constitutive_j2_LpAndItsTangent
use constitutive_phenopowerlaw, only: constitutive_phenopowerlaw_label, &
constitutive_phenopowerlaw_LpAndItsTangent
use constitutive_titanmod, only: constitutive_titanmod_label, &
constitutive_titanmod_LpAndItsTangent
use constitutive_dislotwin, only: constitutive_dislotwin_label, &
constitutive_dislotwin_LpAndItsTangent
use constitutive_nonlocal, only: constitutive_nonlocal_label, &
constitutive_nonlocal_LpAndItsTangent
implicit none
select case (phase_constitution(material_phase(ipc,ip,el)))
case (constitutive_j2_label)
call constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_phenopowerlaw_label)
call constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_titanmod_label)
call constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_dislotwin_label)
call constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_nonlocal_label)
call constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, Temperature, constitutive_state, ipc, ip, el)
end select
!*** input variables ***!
integer(pInt), intent(in):: ipc, & ! component-ID of current integration point
ip, & ! current integration point
el ! current element
real(pReal), intent(in) :: Temperature
real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola-Kirchhoff stress
!*** output variables ***!
real(pReal), dimension(3,3), intent(out) :: Lp ! plastic velocity gradient
real(pReal), dimension(9,9), intent(out) :: dLp_dTstar ! derivative of Lp with respect to Tstar (4th-order tensor)
!*** local variables ***!
select case (phase_constitution(material_phase(ipc,ip,el)))
case (constitutive_j2_label)
call constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_phenopowerlaw_label)
call constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_titanmod_label)
call constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_dislotwin_label)
call constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_nonlocal_label)
call constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, Temperature, constitutive_state, ipc, ip, el)
end select
return
endsubroutine
subroutine constitutive_collectDotState(Tstar_v, Fe, Fp, Temperature, subdt, orientation, ipc, ip, el)
!*********************************************************************
!* This subroutine contains the constitutive equation for *
!* calculating the rate of change of microstructure *
!* INPUT: *
!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) *
!* - state : current microstructure *
!* - ipc : component-ID of current integration point *
!* - ip : current integration point *
!* - el : current element *
!* OUTPUT: *
!* - constitutive_dotState : evolution of state variable *
!*********************************************************************
subroutine constitutive_collectDotState(Tstar_v, Fe, Fp, Temperature, subdt, orientation, ipc, ip, el)
use prec, only: pReal, pInt
use debug, only: debug_cumDotStateCalls, &
debug_cumDotStateTicks
@ -542,16 +554,20 @@ use constitutive_nonlocal, only: constitutive_nonlocal_dotState, &
implicit none
!*** input variables
integer(pInt), intent(in) :: ipc, ip, el
integer(pInt), intent(in) :: ipc, & ! component-ID of current integration point
ip, & ! current integration point
el ! current element
real(pReal), intent(in) :: Temperature, &
subdt
subdt ! timestep
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
Fe, &
Fp
Fe, & ! elastic deformation gradient
Fp ! plastic deformation gradient
real(pReal), dimension(4,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
orientation
orientation ! crystal orientation (quaternion)
real(pReal), dimension(6), intent(in) :: &
Tstar_v
Tstar_v ! 2nd Piola Kirchhoff stress tensor (Mandel)
!*** output variables ***!
!*** local variables
integer(pLongInt) tick, tock, &
@ -587,27 +603,21 @@ call system_clock(count=tock,count_rate=tickrate,count_max=maxticks)
if (tock < tick) debug_cumDotStateTicks = debug_cumDotStateTicks + maxticks
!$OMP END CRITICAL (debugTimingDotState)
return
endsubroutine
function constitutive_dotTemperature(Tstar_v,Temperature,ipc,ip,el)
!*********************************************************************
!* This subroutine contains the constitutive equation for *
!* calculating the rate of change of microstructure *
!* INPUT: *
!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) *
!* - state : current microstructure *
!* - ipc : component-ID of current integration point *
!* - ip : current integration point *
!* - el : current element *
!* OUTPUT: *
!* - constitutive_dotTemperature : evolution of temperature *
!*********************************************************************
function constitutive_dotTemperature(Tstar_v,Temperature,ipc,ip,el)
use prec, only: pReal,pInt
use debug, only: debug_cumDotTemperatureCalls, &
debug_cumDotTemperatureTicks
use material, only: phase_constitution,material_phase
use material, only: phase_constitution, &
material_phase
use constitutive_j2, only: constitutive_j2_dotTemperature, &
constitutive_j2_label
use constitutive_phenopowerlaw, only: constitutive_phenopowerlaw_dotTemperature, &
@ -620,15 +630,22 @@ use constitutive_nonlocal, only: constitutive_nonlocal_dotTemperature, &
constitutive_nonlocal_label
implicit none
!* Definition of variables
integer(pInt) ipc, ip, el
real(pReal) Temperature
real(pReal) constitutive_dotTemperature
real(pReal), dimension(6) :: Tstar_v
integer(pLongInt) tick, &
tock, &
tickrate, &
maxticks
!*** input variables
integer(pInt), intent(in) :: ipc, & ! component-ID of current integration point
ip, & ! current integration point
el ! current element
real(pReal), intent(in) :: Temperature
real(pReal), dimension(6), intent(in) :: &
Tstar_v ! 2nd Piola Kirchhoff stress tensor (Mandel)
!*** output variables ***!
real(pReal) constitutive_dotTemperature ! evolution of temperature
!*** local variables
integer(pLongInt) tick, tock, &
tickrate, &
maxticks
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
@ -658,10 +675,10 @@ call system_clock(count=tock,count_rate=tickrate,count_max=maxticks)
if (tock < tick) debug_cumDotTemperatureTicks = debug_cumDotTemperatureTicks + maxticks
!$OMP END CRITICAL (debugTimingDotTemperature)
return
endfunction
function constitutive_postResults(Tstar_v, Fe, Temperature, dt, ipc, ip, el)
!*********************************************************************
!* return array of constitutive results *
@ -672,49 +689,62 @@ function constitutive_postResults(Tstar_v, Fe, Temperature, dt, ipc, ip, el)
!* - ip : current integration point *
!* - el : current element *
!*********************************************************************
use prec, only: pReal,pInt
use mesh, only: mesh_NcpElems, &
mesh_maxNips, &
mesh_maxNipNeighbors
use material, only: phase_constitution, &
material_phase, &
homogenization_maxNgrains
use constitutive_j2
use constitutive_phenopowerlaw
use constitutive_titanmod
use constitutive_dislotwin
use constitutive_nonlocal
implicit none
use prec, only: pReal,pInt
use mesh, only: mesh_NcpElems, &
mesh_maxNips, &
mesh_maxNipNeighbors
use material, only: phase_constitution, &
material_phase, &
homogenization_maxNgrains
use constitutive_j2, only: constitutive_j2_postResults, &
constitutive_j2_label
use constitutive_phenopowerlaw, only: constitutive_phenopowerlaw_postResults, &
constitutive_phenopowerlaw_label
use constitutive_titanmod, only: constitutive_titanmod_postResults, &
constitutive_titanmod_label
use constitutive_dislotwin, only: constitutive_dislotwin_postResults, &
constitutive_dislotwin_label
use constitutive_nonlocal, only: constitutive_nonlocal_postResults, &
constitutive_nonlocal_label
implicit none
!* Definition of variables
integer(pInt), intent(in) :: ipc,ip,el
real(pReal), intent(in) :: dt, Temperature
real(pReal), dimension(6), intent(in) :: Tstar_v
real(pReal), dimension(3,3), intent(in) :: Fe
real(pReal), dimension(constitutive_sizePostResults(ipc,ip,el)) :: constitutive_postResults
!*** input variables
integer(pInt), intent(in) :: ipc, & ! component-ID of current integration point
ip, & ! current integration point
el ! current element
real(pReal), intent(in) :: Temperature, &
dt ! timestep
real(pReal), dimension(3,3), intent(in) :: &
Fe ! elastic deformation gradient
real(pReal), dimension(6), intent(in) :: &
Tstar_v ! 2nd Piola Kirchhoff stress tensor (Mandel)
constitutive_postResults = 0.0_pReal
select case (phase_constitution(material_phase(ipc,ip,el)))
!*** output variables ***!
real(pReal), dimension(constitutive_sizePostResults(ipc,ip,el)) :: constitutive_postResults
!*** local variables
constitutive_postResults = 0.0_pReal
select case (phase_constitution(material_phase(ipc,ip,el)))
case (constitutive_j2_label)
constitutive_postResults = constitutive_j2_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el)
case (constitutive_phenopowerlaw_label)
constitutive_postResults = constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el)
case (constitutive_titanmod_label)
constitutive_postResults = constitutive_titanmod_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el)
case (constitutive_dislotwin_label)
constitutive_postResults = constitutive_dislotwin_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el)
case (constitutive_nonlocal_label)
constitutive_postResults = constitutive_nonlocal_postResults(Tstar_v, Fe, Temperature, dt, constitutive_state, &
constitutive_dotstate, ipc, ip, el)
end select
case (constitutive_j2_label)
constitutive_postResults = constitutive_j2_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el)
case (constitutive_phenopowerlaw_label)
constitutive_postResults = constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el)
case (constitutive_titanmod_label)
constitutive_postResults = constitutive_titanmod_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el)
case (constitutive_dislotwin_label)
constitutive_postResults = constitutive_dislotwin_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el)
case (constitutive_nonlocal_label)
constitutive_postResults = constitutive_nonlocal_postResults(Tstar_v, Fe, Temperature, dt, constitutive_state, &
constitutive_dotstate, ipc, ip, el)
end select
return
endfunction
END MODULE

10
code/constitutive_j2.f90
View File

@ -75,16 +75,20 @@ subroutine constitutive_j2_init(file)
character(len=64) tag
character(len=1024) line
!$OMP CRITICAL (write2out)
write(6,*)
write(6,'(a20,a20,a12)') '<<<+- constitutive_',constitutive_j2_label,' init -+>>>'
write(6,*) '$Id$'
write(6,*)
!$OMP END CRITICAL (write2out)
maxNinstance = count(phase_constitution == constitutive_j2_label)
if (maxNinstance == 0) return
!$OMP CRITICAL (write2out)
write(6,'(a16,x,i5)') '# instances:',maxNinstance
write(6,*)
!$OMP END CRITICAL (write2out)
allocate(constitutive_j2_sizeDotState(maxNinstance)) ; constitutive_j2_sizeDotState = 0_pInt
allocate(constitutive_j2_sizeState(maxNinstance)) ; constitutive_j2_sizeState = 0_pInt
@ -188,8 +192,8 @@ subroutine constitutive_j2_init(file)
constitutive_j2_Cslip_66(k,k,i) = constitutive_j2_C11(i)
constitutive_j2_Cslip_66(k+3,k+3,i) = 0.5_pReal*(constitutive_j2_C11(i)-constitutive_j2_C12(i))
end forall
constitutive_j2_Cslip_66(:,:,i) = &
math_Mandel3333to66(math_Voigt66to3333(constitutive_j2_Cslip_66(:,:,i)))
constitutive_j2_Cslip_66(1:6,1:6,i) = &
math_Mandel3333to66(math_Voigt66to3333(constitutive_j2_Cslip_66(1:6,1:6,i)))
enddo
@ -258,7 +262,7 @@ function constitutive_j2_homogenizedC(state,ipc,ip,el)
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
matID = phase_constitutionInstance(material_phase(ipc,ip,el))
constitutive_j2_homogenizedC = constitutive_j2_Cslip_66(:,:,matID)
constitutive_j2_homogenizedC = constitutive_j2_Cslip_66(1:6,1:6,matID)
return

118
code/constitutive_nonlocal.f90
View File

@ -170,10 +170,12 @@ character(len=64) tag
character(len=1024) line
write(6,*)
write(6,'(a20,a20,a12)') '<<<+- constitutive_',constitutive_nonlocal_label,' init -+>>>'
write(6,*) '$Id$'
write(6,*)
!$OMP CRITICAL (write2out)
write(6,*)
write(6,'(a20,a20,a12)') '<<<+- constitutive_',constitutive_nonlocal_label,' init -+>>>'
write(6,*) '$Id$'
write(6,*)
!$OMP END CRITICAL (write2out)
maxNinstance = count(phase_constitution == constitutive_nonlocal_label)
if (maxNinstance == 0) return ! we don't have to do anything if there's no instance for this constitutive law
@ -1081,6 +1083,8 @@ real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstan
tauThreshold, & ! threshold shear stress
tau, & ! resolved shear stress
rhoForest ! forest dislocation density
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: &
v ! velocity for the current element and ip
real(pReal) boltzmannProbability, &
tauRel, & ! relative thermally active resolved shear stress
wallFunc, & ! functions reflecting the shape of the obstacle wall (see PhD thesis Mohles p.53)
@ -1096,7 +1100,7 @@ tauThreshold = state%p(11*ns+1:12*ns)
Tdislocation_v = state%p(12*ns+1:12*ns+6)
tau = 0.0_pReal
constitutive_nonlocal_v(1:ns,1:4,g,ip,el) = 0.0_pReal
v = 0.0_pReal
if (present(dv_dtau)) dv_dtau = 0.0_pReal
@ -1124,27 +1128,26 @@ if (Temperature > 0.0_pReal) then
timeRatio = boltzmannProbability * constitutive_nonlocal_fattack(myInstance) &
/ (constitutive_nonlocal_vs(myInstance) * sqrt(rhoForest(s)))
constitutive_nonlocal_v(s,:,g,ip,el) = sign(constitutive_nonlocal_vs(myInstance),tau(s)) * timeRatio / (1.0_pReal + timeRatio)
v(s,1:4) = sign(constitutive_nonlocal_vs(myInstance),tau(s)) * timeRatio / (1.0_pReal + timeRatio)
if (present(dv_dtau)) then
dv_dtau(s) = abs(constitutive_nonlocal_v(s,1,g,ip,el)) * constitutive_nonlocal_Qeff0(s,myInstance) &
/ (kB * Temperature * (1.0_pReal + timeRatio)) &
* 0.5_pReal * wallFunc * (2.0_pReal - tauRel) &
/ ((1.0_pReal - tauRel) * (abs(tau(s)) - tauThreshold(s)))
dv_dtau(s) = abs(v(s,1)) * constitutive_nonlocal_Qeff0(s,myInstance) / (kB * Temperature * (1.0_pReal + timeRatio)) &
* 0.5_pReal * wallFunc * (2.0_pReal - tauRel) / ((1.0_pReal - tauRel) * (abs(tau(s)) - tauThreshold(s)))
endif
!*** If resolved stress exceeds threshold plus obstacle stress, the probability for thermal activation is 1.
!*** The tangent is zero, since no dependency of tau.
elseif (tauRel >= 1.0_pReal) then
constitutive_nonlocal_v(s,1:4,g,ip,el) = sign(constitutive_nonlocal_vs(myInstance), tau(s)) &
* constitutive_nonlocal_fattack(myInstance) &
/ (constitutive_nonlocal_vs(myInstance) * sqrt(rhoForest(s)) &
+ constitutive_nonlocal_fattack(myInstance))
v(s,1:4) = sign(constitutive_nonlocal_vs(myInstance), tau(s)) * constitutive_nonlocal_fattack(myInstance) &
/ (constitutive_nonlocal_vs(myInstance) * sqrt(rhoForest(s)) + constitutive_nonlocal_fattack(myInstance))
endif
enddo
endif
constitutive_nonlocal_v(1:ns,1:4,g,ip,el) = v
!$OMP FLUSH(constitutive_nonlocal_v)
!if (verboseDebugger .and. s) then
! !$OMP CRITICAL (write2out)
! write(6,*) '::: kinetics',g,ip,el
@ -1234,17 +1237,21 @@ myInstance = phase_constitutionInstance(material_phase(g,ip,el))
myStructure = constitutive_nonlocal_structure(myInstance)
ns = constitutive_nonlocal_totalNslip(myInstance)
!*** update dislocation velocity
call constitutive_nonlocal_kinetics(Tstar_v, Temperature, state(g,ip,el), g, ip, el, dv_dtau)
!*** shortcut to state variables
forall (t = 1:8) &
rhoSgl(1:ns,t) = state(g,ip,el)%p((t-1)*ns+1:t*ns)
forall (s = 1:ns, t = 5:8, rhoSgl(s,t) * constitutive_nonlocal_v(s,t-4,g,ip,el) < 0.0_pReal) & ! contribution of used rho for changing sign of v
rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t))
rhoForest = state(g,ip,el)%p(10*ns+1:11*ns)
tauThreshold = state(g,ip,el)%p(11*ns+1:12*ns)
call constitutive_nonlocal_kinetics(Tstar_v, Temperature, state(g,ip,el), g, ip, el, dv_dtau) ! update dislocation velocity
!*** Calculation of gdot and its tangent
@ -1252,23 +1259,21 @@ forall (t = 1:4) &
gdot(1:ns,t) = rhoSgl(1:ns,t) * constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) &
* constitutive_nonlocal_v(1:ns,t,g,ip,el)
gdotTotal = sum(gdot,2)
dgdotTotal_dtau = sum(rhoSgl,2) * constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) * dv_dtau
!*** Calculation of Lp and its tangent
do s = 1,ns
sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance)
sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance)
Lp = Lp + gdotTotal(s) * lattice_Sslip(1:3,1:3,sLattice,myStructure)
forall (i=1:3,j=1:3,k=1:3,l=1:3) &
dLp_dTstar3333(i,j,k,l) = dLp_dTstar3333(i,j,k,l) + dgdotTotal_dtau(s) * lattice_Sslip(i,j, sLattice,myStructure) &
* lattice_Sslip(k,l, sLattice,myStructure)
enddo
dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333)
!if (verboseDebugger .and. (debug_g==g .and. debug_i==i .and. debug_e==e)) then
! !$OMP CRITICAL (write2out)
! write(6,*) '::: LpandItsTangent',g,ip,el
@ -1448,6 +1453,7 @@ gdot = 0.0_pReal
dLower = 0.0_pReal
dUpper = 0.0_pReal
!*** shortcut to state variables
forall (t = 1:8) rhoSgl(1:ns,t) = state(g,ip,el)%p((t-1)*ns+1:t*ns)
@ -1456,6 +1462,7 @@ rhoForest = state(g,ip,el)%p(10*ns+1:11*ns)
tauThreshold = state(g,ip,el)%p(11*ns+1:12*ns)
Tdislocation_v = state(g,ip,el)%p(12*ns+1:12*ns+6)
!*** sanity check for timestep
if (timestep <= 0.0_pReal) then ! if illegal timestep...
@ -1464,6 +1471,7 @@ if (timestep <= 0.0_pReal) then
endif
!****************************************************************************
!*** Calculate shear rate
@ -1485,6 +1493,7 @@ if (verboseDebugger .and. (debug_g==g .and. debug_i==ip .and. debug_e==el)) then
endif
!****************************************************************************
!*** calculate limits for stable dipole height
@ -1519,6 +1528,7 @@ if (timestep > 0.0_pReal) then
endif
!****************************************************************************
!*** calculate dislocation multiplication
@ -1533,6 +1543,7 @@ where (rhoSgl(1:ns,1:2) > 0.0_pReal) &
/ constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance), 2, 2)
!****************************************************************************
!*** calculate dislocation fluxes (only for nonlocal constitution)
@ -1656,9 +1667,11 @@ if (.not. phase_localConstitution(material_phase(g,ip,el))) then
enddo ! neighbor loop
endif
if (numerics_integrationMode == 1_pInt) &
if (numerics_integrationMode == 1_pInt) then
constitutive_nonlocal_rhoDotFlux(1:ns,1:10,g,ip,el) = rhoDotFlux(1:ns,1:10) ! save flux calculation for output (if in central integration mode)
endif
!****************************************************************************
!*** calculate dipole formation and annihilation
@ -1740,9 +1753,7 @@ if (verboseDebugger .and. (debug_g==g .and. debug_i==ip .and. debug_e==el)) then
!$OMP END CRITICAL (write2out)
endif
!$OMP CRITICAL (copy2dotState)
dotState(g,ip,el)%p(1:10*ns) = dotState(g,ip,el)%p(1:10*ns) + reshape(rhoDot,(/10*ns/))
!$OMP END CRITICAL (copy2dotState)
dotState(g,ip,el)%p(1:10*ns) = dotState(g,ip,el)%p(1:10*ns) + reshape(rhoDot,(/10*ns/))
endsubroutine
@ -1804,8 +1815,10 @@ integer(pInt) Nneighbors, & !
s1, & ! slip system index (me)
s2 ! slip system index (my neighbor)
real(pReal), dimension(4) :: absoluteMisorientation ! absolute misorientation (without symmetry) between me and my neighbor
real(pReal), dimension(2,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(1,i,e)))) :: &
compatibility ! compatibility of one specific slip system to all neighbors slip systems's for edges and screws
real(pReal), dimension(2,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(1,i,e))),&
constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(1,i,e))),&
FE_NipNeighbors(mesh_element(2,e))) :: &
compatibility ! compatibility for current element and ip
real(pReal), dimension(3,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(1,i,e)))) :: &
slipNormal, &
slipDirection
@ -1820,20 +1833,19 @@ Nneighbors = FE_NipNeighbors(mesh_element(2,e))
my_phase = material_phase(1,i,e)
my_instance = phase_constitutionInstance(my_phase)
my_structure = constitutive_nonlocal_structure(my_instance)
ns = constitutive_nonlocal_totalNslip(my_instance)
ns = constitutive_nonlocal_totalNslip(my_instance)
slipNormal(1:3,1:ns) = lattice_sn(1:3, constitutive_nonlocal_slipSystemLattice(1:ns,my_instance), my_structure)
slipDirection(1:3,1:ns) = lattice_sd(1:3, constitutive_nonlocal_slipSystemLattice(1:ns,my_instance), my_structure)
!*** start out fully compatible
constitutive_nonlocal_compatibility(:,:,:,:,i,e) = 0.0_pReal
forall(s1 = 1:maxval(constitutive_nonlocal_totalNslip)) &
constitutive_nonlocal_compatibility(1:2,s1,s1,1:Nneighbors,i,e) = 1.0_pReal
compatibility = 0.0_pReal
forall(s1 = 1:ns) &
compatibility(1:2,s1,s1,1:Nneighbors) = 1.0_pReal
!*** Loop thrugh neighbors and check whether there is any compatibility.
!*** This is only the case for
do n = 1,Nneighbors
neighboring_e = mesh_ipNeighborhood(1,n,i,e)
@ -1845,21 +1857,21 @@ do n = 1,Nneighbors
if (neighboring_e <= 0 .or. neighboring_i <= 0) then
forall(s1 = 1:ns) &
constitutive_nonlocal_compatibility(1:2,s1,s1,n,i,e) = sqrt(constitutive_nonlocal_surfaceTransmissivity(my_instance))
compatibility(1:2,s1,s1,n) = sqrt(constitutive_nonlocal_surfaceTransmissivity(my_instance))
cycle
endif
!* PHASE BOUNDARY
!* If we encounter a different nonlocal "cpfem" phase at the neighbor,
!* we consider this to be a real "physical" phase boundary, so fully incompatible.
!* we consider this to be a real "physical" phase boundary, so completely incompatible.
!* If the neighboring "cpfem" phase has a local constitution,
!* we do not consider this to be a phase boundary, so fully compatible.
!* we do not consider this to be a phase boundary, so completely compatible.
neighboring_phase = material_phase(1,neighboring_i,neighboring_e)
if (neighboring_phase /= my_phase) then
if (.not. phase_localConstitution(neighboring_phase)) then
constitutive_nonlocal_compatibility(:,:,:,n,i,e) = 0.0_pReal
compatibility(1:2,1:ns,1:ns,n) = 0.0_pReal
endif
cycle
endif
@ -1879,35 +1891,32 @@ do n = 1,Nneighbors
0_pInt) ! no symmetry
do s1 = 1,ns ! my slip systems
do s2 = 1,ns ! my neighbor's slip systems
compatibility(1,s2) = math_mul3x3(slipNormal(1:3,s1), math_qRot(absoluteMisorientation, slipNormal(1:3,s2))) &
* abs(math_mul3x3(slipDirection(1:3,s1), math_qRot(absoluteMisorientation, slipDirection(1:3,s2))))
compatibility(2,s2) = abs(math_mul3x3(slipNormal(1:3,s1), math_qRot(absoluteMisorientation, slipNormal(1:3,s2)))) &
* abs(math_mul3x3(slipDirection(1:3,s1), math_qRot(absoluteMisorientation, slipDirection(1:3,s2))))
compatibility(1,s2,s1,n) = math_mul3x3(slipNormal(1:3,s1), math_qRot(absoluteMisorientation, slipNormal(1:3,s2))) &
* abs(math_mul3x3(slipDirection(1:3,s1), math_qRot(absoluteMisorientation, slipDirection(1:3,s2))))
compatibility(2,s2,s1,n) = abs(math_mul3x3(slipNormal(1:3,s1), math_qRot(absoluteMisorientation, slipNormal(1:3,s2)))) &
* abs(math_mul3x3(slipDirection(1:3,s1), math_qRot(absoluteMisorientation, slipDirection(1:3,s2))))
enddo
compatibilitySum = 0.0_pReal
belowThreshold = .true.
do while (compatibilitySum < 1.0_pReal .and. any(belowThreshold(1:ns)))
thresholdValue = maxval(compatibility(2,1:ns), belowThreshold(1:ns)) ! screws always positive
nThresholdValues = dble(count(compatibility(2,1:ns) == thresholdValue))
where (compatibility(2,1:ns) >= thresholdValue) &
thresholdValue = maxval(compatibility(2,1:ns,s1,n), belowThreshold(1:ns)) ! screws always positive
nThresholdValues = dble(count(compatibility(2,1:ns,s1,n) == thresholdValue))
where (compatibility(2,1:ns,s1,n) >= thresholdValue) &
belowThreshold(1:ns) = .false.
if (compatibilitySum + thresholdValue * nThresholdValues > 1.0_pReal) &
where (abs(compatibility(1:2,1:ns)) == thresholdValue) &
compatibility(1:2,1:ns) = sign((1.0_pReal - compatibilitySum) / nThresholdValues, compatibility(1:2,1:ns))
where (abs(compatibility(1:2,1:ns,s1,n)) == thresholdValue) &
compatibility(1:2,1:ns,s1,n) = sign((1.0_pReal - compatibilitySum) / nThresholdValues, compatibility(1:2,1:ns,s1,n))
compatibilitySum = compatibilitySum + nThresholdValues * thresholdValue
enddo
where (belowThreshold(1:ns)) compatibility(1,1:ns) = 0.0_pReal
where (belowThreshold(1:ns)) compatibility(2,1:ns) = 0.0_pReal
constitutive_nonlocal_compatibility(1:2,1:ns,s1,n,i,e) = compatibility(1:2,1:ns)
where (belowThreshold(1:ns)) compatibility(1,1:ns,s1,n) = 0.0_pReal
where (belowThreshold(1:ns)) compatibility(2,1:ns,s1,n) = 0.0_pReal
enddo ! my slip systems cycle
enddo ! neighbor cycle
constitutive_nonlocal_compatibility(1:2,1:ns,1:ns,1:Nneighbors,i,e) = compatibility
endsubroutine
@ -2412,4 +2421,5 @@ do o = 1,phase_Noutput(material_phase(g,ip,el))
enddo
endfunction
END MODULE

35
code/constitutive_phenopowerlaw.f90
View File

@ -150,16 +150,20 @@ subroutine constitutive_phenopowerlaw_init(file)
character(len=64) tag,formatting
character(len=1024) line
!$OMP CRITICAL (write2out)
write(6,*)
write(6,'(a20,a20,a12)') '<<<+- constitutive_',constitutive_phenopowerlaw_label,' init -+>>>'
write(6,*) '$Id$'
write(6,*)
!$OMP END CRITICAL (write2out)
maxNinstance = count(phase_constitution == constitutive_phenopowerlaw_label)
if (maxNinstance == 0) return
!$OMP CRITICAL (write2out)
write(6,'(a16,x,i5)') '# instances:',maxNinstance
write(6,*)
!$OMP END CRITICAL (write2out)
allocate(constitutive_phenopowerlaw_sizeDotState(maxNinstance)) ; constitutive_phenopowerlaw_sizeDotState = 0_pInt
allocate(constitutive_phenopowerlaw_sizeState(maxNinstance)) ; constitutive_phenopowerlaw_sizeState = 0_pInt
@ -656,11 +660,11 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temp
!* Calculation of Lp
tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,structID))
tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(1:6,index_myFamily+i,structID))
gdot_slip(j) = constitutive_phenopowerlaw_gdot0_slip(matID)*(abs(tau_slip(j))/state(ipc,ip,el)%p(j))**&
constitutive_phenopowerlaw_n_slip(matID)*sign(1.0_pReal,tau_slip(j))
Lp = Lp + (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F
gdot_slip(j)*lattice_Sslip(:,:,index_myFamily+i,structID)
gdot_slip(j)*lattice_Sslip(1:3,1:3,index_myFamily+i,structID)
!* Calculation of the tangent of Lp
@ -682,12 +686,12 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temp
!* Calculation of Lp
tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,structID))
tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,structID))
gdot_twin(j) = (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F
constitutive_phenopowerlaw_gdot0_twin(matID)*&
(abs(tau_twin(j))/state(ipc,ip,el)%p(nSlip+j))**&
constitutive_phenopowerlaw_n_twin(matID)*max(0.0_pReal,sign(1.0_pReal,tau_twin(j)))
Lp = Lp + gdot_twin(j)*lattice_Stwin(:,:,index_myFamily+i,structID)
Lp = Lp + gdot_twin(j)*lattice_Stwin(1:3,1:3,index_myFamily+i,structID)
!* Calculation of the tangent of Lp
@ -770,7 +774,6 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el
do i = 1,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family
j = j+1_pInt
h_slipslip(j) = c_slipslip*(1.0_pReal-state(ipc,ip,el)%p(j) / & ! system-dependent prefactor for slip--slip interaction
(constitutive_phenopowerlaw_tausat_slip(f,matID)+ssat_offset))** &
constitutive_phenopowerlaw_w0_slip(matID)
@ -778,7 +781,7 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el
!* Calculation of dot gamma
tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,structID))
tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(1:6,index_myFamily+i,structID))
gdot_slip(j) = constitutive_phenopowerlaw_gdot0_slip(matID)*(abs(tau_slip(j))/state(ipc,ip,el)%p(j))**&
constitutive_phenopowerlaw_n_slip(matID)*sign(1.0_pReal,tau_slip(j))
enddo
@ -794,7 +797,7 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el
!* Calculation of dot vol frac
tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,structID))
tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,structID))
gdot_twin(j) = (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F
constitutive_phenopowerlaw_gdot0_twin(matID)*&
(abs(tau_twin(j))/state(ipc,ip,el)%p(nSlip+j))**&
@ -808,9 +811,9 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el
do f = 1,lattice_maxNslipFamily ! loop over all slip families
do i = 1,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family
j = j+1_pInt
constitutive_phenopowerlaw_dotState(j) = & ! evolution of slip resistance j
h_slipslip(j) * dot_product(constitutive_phenopowerlaw_hardeningMatrix_slipslip(:,j,matID),abs(gdot_slip)) + & ! dot gamma_slip
h_sliptwin(j) * dot_product(constitutive_phenopowerlaw_hardeningMatrix_sliptwin(:,j,matID),gdot_twin) ! dot gamma_twin
constitutive_phenopowerlaw_dotState(j) = & ! evolution of slip resistance j
h_slipslip(j) * dot_product(constitutive_phenopowerlaw_hardeningMatrix_slipslip(1:nSlip,j,matID),abs(gdot_slip)) + & ! dot gamma_slip
h_sliptwin(j) * dot_product(constitutive_phenopowerlaw_hardeningMatrix_sliptwin(1:nTwin,j,matID),gdot_twin) ! dot gamma_twin
constitutive_phenopowerlaw_dotState(index_Gamma) = constitutive_phenopowerlaw_dotState(index_Gamma) + &
abs(gdot_slip(j))
enddo
@ -821,9 +824,9 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el
index_myFamily = sum(lattice_NtwinSystem(1:f-1,structID)) ! at which index starts my family
do i = 1,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family
j = j+1_pInt
constitutive_phenopowerlaw_dotState(j+nSlip) = & ! evolution of twin resistance j
h_twinslip(j) * dot_product(constitutive_phenopowerlaw_hardeningMatrix_twinslip(:,j,matID),abs(gdot_slip)) + & ! dot gamma_slip
h_twintwin(j) * dot_product(constitutive_phenopowerlaw_hardeningMatrix_twintwin(:,j,matID),gdot_twin) ! dot gamma_twin
constitutive_phenopowerlaw_dotState(j+nSlip) = & ! evolution of twin resistance j
h_twinslip(j) * dot_product(constitutive_phenopowerlaw_hardeningMatrix_twinslip(1:nSlip,j,matID),abs(gdot_slip)) + & ! dot gamma_slip
h_twintwin(j) * dot_product(constitutive_phenopowerlaw_hardeningMatrix_twintwin(1:nTwin,j,matID),gdot_twin) ! dot gamma_twin
constitutive_phenopowerlaw_dotState(index_F) = constitutive_phenopowerlaw_dotState(index_F) + &
gdot_twin(j)/lattice_shearTwin(index_myFamily+i,structID)
enddo
@ -929,7 +932,7 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,stat
index_myFamily = sum(lattice_NslipSystem(1:f-1,structID)) ! at which index starts my family
do i = 1,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family
j = j + 1_pInt
constitutive_phenopowerlaw_postResults(c+j) = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,structID))
constitutive_phenopowerlaw_postResults(c+j) = dot_product(Tstar_v,lattice_Sslip_v(1:6,index_myFamily+i,structID))
enddo; enddo
c = c + nSlip
@ -947,7 +950,7 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,stat
index_myFamily = sum(lattice_NtwinSystem(1:f-1,structID)) ! at which index starts my family
do i = 1,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family
j = j + 1_pInt
tau = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,structID))
tau = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,structID))
constitutive_phenopowerlaw_postResults(c+j) = (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F
constitutive_phenopowerlaw_gdot0_twin(matID)*&
(abs(tau)/state(ipc,ip,el)%p(j+nSlip))**&
@ -961,7 +964,7 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,stat
index_myFamily = sum(lattice_NtwinSystem(1:f-1,structID)) ! at which index starts my family
do i = 1,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family
j = j + 1_pInt
constitutive_phenopowerlaw_postResults(c+j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,structID))
constitutive_phenopowerlaw_postResults(c+j) = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,structID))
enddo; enddo
c = c + nTwin

1434
code/crystallite.f90
View File

File diff suppressed because it is too large Load Diff

40
code/debug.f90
View File

@ -62,10 +62,12 @@ subroutine debug_init()
character(len=64) tag
character(len=1024) line
write(6,*)
write(6,*) '<<<+- debug init -+>>>'
write(6,*) '$Id$'
write(6,*)
!$OMP CRITICAL (write2out)
write(6,*)
write(6,*) '<<<+- debug init -+>>>'
write(6,*) '$Id$'
write(6,*)
!$OMP END CRITICAL (write2out)
allocate(debug_StressLoopDistribution(nStress,2)) ; debug_StressLoopDistribution = 0_pInt
allocate(debug_LeapfrogBreakDistribution(nStress,2)) ; debug_LeapfrogBreakDistribution = 0_pInt
@ -77,8 +79,10 @@ subroutine debug_init()
! try to open the config file
if(IO_open_file(fileunit,debug_configFile)) then
write(6,*) ' ... using values from config file'
write(6,*)
!$OMP CRITICAL (write2out)
write(6,*) ' ... using values from config file'
write(6,*)
!$OMP END CRITICAL (write2out)
line = ''
! read variables from config file and overwrite parameters
@ -107,19 +111,25 @@ subroutine debug_init()
! no config file, so we use standard values
else
write(6,*) ' ... using standard values'
write(6,*)
!$OMP CRITICAL (write2out)
write(6,*) ' ... using standard values'
write(6,*)
!$OMP END CRITICAL (write2out)
endif
! writing parameters to output file
write(6,'(a24,x,l)') 'debug: ',debugger
write(6,'(a24,x,l)') 'verbose: ',verboseDebugger
write(6,'(a24,x,l)') 'selective: ',selectiveDebugger
!$OMP CRITICAL (write2out)
write(6,'(a24,x,l)') 'debug: ',debugger
write(6,'(a24,x,l)') 'verbose: ',verboseDebugger
write(6,'(a24,x,l)') 'selective: ',selectiveDebugger
!$OMP END CRITICAL (write2out)
if (selectiveDebugger) then
write(6,'(a24,x,i8)') ' element: ',debug_e
write(6,'(a24,x,i8)') ' ip: ',debug_i
write(6,'(a24,x,i8)') ' grain: ',debug_g
!$OMP CRITICAL (write2out)
write(6,'(a24,x,i8)') ' element: ',debug_e
write(6,'(a24,x,i8)') ' ip: ',debug_i
write(6,'(a24,x,i8)') ' grain: ',debug_g
!$OMP END CRITICAL (write2out)
else
debug_e = 0_pInt ! switch off selective debugging
debug_i = 0_pInt
@ -170,6 +180,7 @@ endsubroutine
call system_clock(count_rate=tickrate)
!$OMP CRITICAL (write2out)
write(6,*)
write(6,*) 'DEBUG Info'
write(6,*)
@ -265,6 +276,7 @@ endsubroutine
write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_MaterialpointLoopDistribution)
write(6,*)
!$OMP END CRITICAL (write2out)
endsubroutine

166
code/homogenization.f90
View File

@ -61,7 +61,7 @@ subroutine homogenization_init(Temperature)
use constitutive, only: constitutive_maxSizePostResults
use crystallite, only: crystallite_maxSizePostResults
use homogenization_isostrain
use homogenization_RGC ! RGC homogenization added <<<updated 31.07.2009>>>
use homogenization_RGC
real(pReal) Temperature
integer(pInt), parameter :: fileunit = 200
@ -73,7 +73,7 @@ subroutine homogenization_init(Temperature)
if(.not. IO_open_file(fileunit,material_configFile)) call IO_error (100) ! corrupt config file
call homogenization_isostrain_init(fileunit) ! parse all homogenizations of this type
call homogenization_RGC_init(fileunit) ! RGC homogenization added <<<updated 31.07.2009>>>
call homogenization_RGC_init(fileunit)
close(fileunit)
@ -131,8 +131,8 @@ subroutine homogenization_init(Temperature)
allocate(materialpoint_doneAndHappy(2,mesh_maxNips,mesh_NcpElems)); materialpoint_doneAndHappy = .true.
forall (i = 1:mesh_maxNips,e = 1:mesh_NcpElems)
materialpoint_F0(:,:,i,e) = math_I3
materialpoint_F(:,:,i,e) = math_I3
materialpoint_F0(1:3,1:3,i,e) = math_I3
materialpoint_F(1:3,1:3,i,e) = math_I3
end forall
do e = 1,mesh_NcpElems ! loop over elements
@ -149,7 +149,7 @@ subroutine homogenization_init(Temperature)
homogenization_sizeState(i,e) = homogenization_isostrain_sizeState(myInstance)
endif
homogenization_sizePostResults(i,e) = homogenization_isostrain_sizePostResults(myInstance)
!* RGC homogenization: added <<<updated 31.07.2009>>>
!* RGC homogenization
case (homogenization_RGC_label)
if (homogenization_RGC_sizeState(myInstance) > 0_pInt) then
allocate(homogenization_state0(i,e)%p(homogenization_RGC_sizeState(myInstance)))
@ -229,7 +229,8 @@ subroutine materialpoint_stressAndItsTangent(&
stepIncreaseHomog, &
nHomog, &
nMPstate
use math, only: math_det3x3
use math, only: math_det3x3, &
math_transpose3x3
use FEsolving, only: FEsolving_execElem, &
FEsolving_execIP, &
terminallyIll
@ -282,10 +283,10 @@ subroutine materialpoint_stressAndItsTangent(&
write (6,*)
write (6,*) 'Material Point start'
write (6,'(a,/,(f14.9,x))') 'Temp0 of 1 1',materialpoint_Temperature(1,1)
write (6,'(a,/,3(3(f14.9,x)/))') 'F0 of 1 1',materialpoint_F0(1:3,:,1,1)
write (6,'(a,/,3(3(f14.9,x)/))') 'F of 1 1',materialpoint_F(1:3,:,1,1)
write (6,'(a,/,3(3(f14.9,x)/))') 'Fp0 of 1 1 1',crystallite_Fp0(1:3,:,1,1,1)
write (6,'(a,/,3(3(f14.9,x)/))') 'Lp0 of 1 1 1',crystallite_Lp0(1:3,:,1,1,1)
write (6,'(a,/,3(3(f14.9,x)/))') 'F0 of 1 1',math_transpose3x3(materialpoint_F0(1:3,1:3,1,1))
write (6,'(a,/,3(3(f14.9,x)/))') 'F of 1 1',math_transpose3x3(materialpoint_F(1:3,1:3,1,1))
write (6,'(a,/,3(3(f14.9,x)/))') 'Fp0 of 1 1 1',math_transpose3x3(crystallite_Fp0(1:3,1:3,1,1,1))
write (6,'(a,/,3(3(f14.9,x)/))') 'Lp0 of 1 1 1',math_transpose3x3(crystallite_Lp0(1:3,1:3,1,1,1))
endif
@ -297,16 +298,16 @@ subroutine materialpoint_stressAndItsTangent(&
! initialize restoration points of grain...
forall (g = 1:myNgrains) constitutive_partionedState0(g,i,e)%p = constitutive_state0(g,i,e)%p ! ...microstructures
crystallite_partionedTemperature0(1:myNgrains,i,e) = materialpoint_Temperature(i,e) ! ...temperatures
crystallite_partionedFp0(:,:,1:myNgrains,i,e) = crystallite_Fp0(:,:,1:myNgrains,i,e) ! ...plastic def grads
crystallite_partionedLp0(:,:,1:myNgrains,i,e) = crystallite_Lp0(:,:,1:myNgrains,i,e) ! ...plastic velocity grads
crystallite_partioneddPdF0(:,:,:,:,1:myNgrains,i,e) = crystallite_dPdF0(:,:,:,:,1:myNgrains,i,e) ! ...stiffness
crystallite_partionedF0(:,:,1:myNgrains,i,e) = crystallite_F0(:,:,1:myNgrains,i,e) ! ...def grads
crystallite_partionedTstar0_v(:,1:myNgrains,i,e)= crystallite_Tstar0_v(:,1:myNgrains,i,e) ! ...2nd PK stress
crystallite_partionedFp0(1:3,1:3,1:myNgrains,i,e) = crystallite_Fp0(1:3,1:3,1:myNgrains,i,e) ! ...plastic def grads
crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e) = crystallite_Lp0(1:3,1:3,1:myNgrains,i,e) ! ...plastic velocity grads
crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,1:myNgrains,i,e) = crystallite_dPdF0(1:3,1:3,1:3,1:3,1:myNgrains,i,e) ! ...stiffness
crystallite_partionedF0(1:3,1:3,1:myNgrains,i,e) = crystallite_F0(1:3,1:3,1:myNgrains,i,e) ! ...def grads
crystallite_partionedTstar0_v(1:6,1:myNgrains,i,e) = crystallite_Tstar0_v(1:6,1:myNgrains,i,e) ! ...2nd PK stress
! initialize restoration points of ...
if (homogenization_sizeState(i,e) > 0_pInt) &
homogenization_subState0(i,e)%p = homogenization_state0(i,e)%p ! ...internal homogenization state
materialpoint_subF0(:,:,i,e) = materialpoint_F0(:,:,i,e) ! ...def grad
materialpoint_subF0(1:3,1:3,i,e) = materialpoint_F0(1:3,1:3,i,e) ! ...def grad
materialpoint_subFrac(i,e) = 0.0_pReal
materialpoint_subStep(i,e) = 1.0_pReal/subStepSizeHomog ! <<added to adopt flexibility in cutback size>>
@ -339,23 +340,26 @@ subroutine materialpoint_stressAndItsTangent(&
! calculate new subStep and new subFrac
materialpoint_subFrac(i,e) = materialpoint_subFrac(i,e) + materialpoint_subStep(i,e)
!$OMP FLUSH(materialpoint_subFrac)
materialpoint_subStep(i,e) = min(1.0_pReal-materialpoint_subFrac(i,e), &
stepIncreaseHomog*materialpoint_subStep(i,e)) ! <<introduce flexibility for step increase/acceleration>>
!$OMP FLUSH(materialpoint_subStep)
! still stepping needed
if (materialpoint_subStep(i,e) > subStepMinHomog) then
! wind forward grain starting point of...
crystallite_partionedTemperature0(1:myNgrains,i,e) = crystallite_Temperature(1:myNgrains,i,e) ! ...temperatures
crystallite_partionedF0(:,:,1:myNgrains,i,e) = crystallite_partionedF(:,:,1:myNgrains,i,e) ! ...def grads
crystallite_partionedFp0(:,:,1:myNgrains,i,e) = crystallite_Fp(:,:,1:myNgrains,i,e) ! ...plastic def grads
crystallite_partionedLp0(:,:,1:myNgrains,i,e) = crystallite_Lp(:,:,1:myNgrains,i,e) ! ...plastic velocity grads
crystallite_partioneddPdF0(:,:,:,:,1:myNgrains,i,e) = crystallite_dPdF(:,:,:,:,1:myNgrains,i,e)! ...stiffness
crystallite_partionedTstar0_v(:,1:myNgrains,i,e) = crystallite_Tstar_v(:,1:myNgrains,i,e) ! ...2nd PK stress
crystallite_partionedF0(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedF(1:3,1:3,1:myNgrains,i,e) ! ...def grads
crystallite_partionedFp0(1:3,1:3,1:myNgrains,i,e) = crystallite_Fp(1:3,1:3,1:myNgrains,i,e) ! ...plastic def grads
crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e) = crystallite_Lp(1:3,1:3,1:myNgrains,i,e) ! ...plastic velocity grads
crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,1:myNgrains,i,e) = crystallite_dPdF(1:3,1:3,1:3,1:3,1:myNgrains,i,e)! ...stiffness
crystallite_partionedTstar0_v(1:6,1:myNgrains,i,e) = crystallite_Tstar_v(1:6,1:myNgrains,i,e) ! ...2nd PK stress
forall (g = 1:myNgrains) constitutive_partionedState0(g,i,e)%p = constitutive_state(g,i,e)%p ! ...microstructures
if (homogenization_sizeState(i,e) > 0_pInt) &
homogenization_subState0(i,e)%p = homogenization_state(i,e)%p ! ...internal state of homog scheme
materialpoint_subF0(:,:,i,e) = materialpoint_subF(:,:,i,e) ! ...def grad
materialpoint_subF0(1:3,1:3,i,e) = materialpoint_subF(1:3,1:3,i,e) ! ...def grad
!$OMP FLUSH(materialpoint_subF0)
elseif (materialpoint_requested(i,e)) then ! this materialpoint just converged ! already at final time (??)
!$OMP CRITICAL (distributionHomog)
debug_MaterialpointLoopDistribution(min(nHomog+1,NiterationHomog)) = &
@ -373,7 +377,7 @@ subroutine materialpoint_stressAndItsTangent(&
!$OMP END CRITICAL (setTerminallyIll)
else ! cutback makes sense
materialpoint_subStep(i,e) = subStepSizeHomog * materialpoint_subStep(i,e) ! crystallite had severe trouble, so do a significant cutback
! <<modified to add more flexibility in cutback>>
!$OMP FLUSH(materialpoint_subStep)
if (verboseDebugger .and. (e == debug_e .and. i == debug_i)) then
!$OMP CRITICAL (write2out)
@ -385,10 +389,10 @@ subroutine materialpoint_stressAndItsTangent(&
! restore...
crystallite_Temperature(1:myNgrains,i,e) = crystallite_partionedTemperature0(1:myNgrains,i,e) ! ...temperatures
! ...initial def grad unchanged
crystallite_Fp(:,:,1:myNgrains,i,e) = crystallite_partionedFp0(:,:,1:myNgrains,i,e) ! ...plastic def grads
crystallite_Lp(:,:,1:myNgrains,i,e) = crystallite_partionedLp0(:,:,1:myNgrains,i,e) ! ...plastic velocity grads
crystallite_dPdF(:,:,:,:,1:myNgrains,i,e) = crystallite_partioneddPdF0(:,:,:,:,1:myNgrains,i,e) ! ...stiffness
crystallite_Tstar_v(:,1:myNgrains,i,e) = crystallite_partionedTstar0_v(:,1:myNgrains,i,e) ! ...2nd PK stress
crystallite_Fp(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedFp0(1:3,1:3,1:myNgrains,i,e) ! ...plastic def grads
crystallite_Lp(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e) ! ...plastic velocity grads
crystallite_dPdF(1:3,1:3,1:3,1:3,1:myNgrains,i,e) = crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,1:myNgrains,i,e) ! ...stiffness
crystallite_Tstar_v(1:6,1:myNgrains,i,e) = crystallite_partionedTstar0_v(1:6,1:myNgrains,i,e) ! ...2nd PK stress
forall (g = 1:myNgrains) constitutive_state(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructures
if (homogenization_sizeState(i,e) > 0_pInt) &
homogenization_state(i,e)%p = homogenization_subState0(i,e)%p ! ...internal state of homog scheme
@ -397,16 +401,16 @@ subroutine materialpoint_stressAndItsTangent(&
materialpoint_requested(i,e) = materialpoint_subStep(i,e) > subStepMinHomog
if (materialpoint_requested(i,e)) then
materialpoint_subF(:,:,i,e) = materialpoint_subF0(:,:,i,e) + &
materialpoint_subStep(i,e) * (materialpoint_F(:,:,i,e) - materialpoint_F0(:,:,i,e))
materialpoint_subdt(i,e) = materialpoint_subStep(i,e) * dt
materialpoint_doneAndHappy(:,i,e) = (/.false.,.true./)
materialpoint_subF(1:3,1:3,i,e) = materialpoint_subF0(1:3,1:3,i,e) + &
materialpoint_subStep(i,e) * (materialpoint_F(1:3,1:3,i,e) - materialpoint_F0(1:3,1:3,i,e))
materialpoint_subdt(i,e) = materialpoint_subStep(i,e) * dt
materialpoint_doneAndHappy(1:2,i,e) = (/.false.,.true./)
endif
enddo ! loop IPs
enddo ! loop elements
!$OMP END PARALLEL DO
!* Checks for cutback/substepping loops: added <<<updated 31.07.2009>>>
!* Checks for cutback/substepping loops
! write (6,'(a,/,8(L,x))') 'MP exceeds substep min',materialpoint_subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMinHomog
! write (6,'(a,/,8(L,x))') 'MP requested',materialpoint_requested(:,FEsolving_execELem(1):FEsolving_execElem(2))
! write (6,'(a,/,8(f6.4,x))') 'MP subFrac',materialpoint_subFrac(:,FEsolving_execELem(1):FEsolving_execElem(2))
@ -467,11 +471,13 @@ subroutine materialpoint_stressAndItsTangent(&
if ( materialpoint_requested(i,e) .and. &
.not. materialpoint_doneAndHappy(1,i,e)) then
if (.not. all(crystallite_converged(:,i,e))) then
materialpoint_doneAndHappy(:,i,e) = (/.true.,.false./)
materialpoint_doneAndHappy(1:2,i,e) = (/.true.,.false./)
materialpoint_converged(i,e) = .false.
else
materialpoint_doneAndHappy(:,i,e) = homogenization_updateState(i,e)
materialpoint_doneAndHappy(1:2,i,e) = homogenization_updateState(i,e)
materialpoint_converged(i,e) = all(homogenization_updateState(i,e)) ! converged if done and happy
endif
materialpoint_converged(i,e) = all(materialpoint_doneAndHappy(:,i,e)) ! converged if done and happy
!$OMP FLUSH(materialpoint_converged)
if (materialpoint_converged(i,e)) then
!$OMP CRITICAL (distributionMPState)
debug_MaterialpointStateLoopdistribution(NiterationMPstate) = &
@ -573,7 +579,7 @@ subroutine homogenization_partitionDeformation(&
use material, only: homogenization_type, homogenization_maxNgrains
use crystallite, only: crystallite_partionedF0,crystallite_partionedF
use homogenization_isostrain
use homogenization_RGC ! RGC homogenization added <<<updated 31.07.2009>>>
use homogenization_RGC
implicit none
@ -582,17 +588,17 @@ subroutine homogenization_partitionDeformation(&
select case(homogenization_type(mesh_element(3,el)))
case (homogenization_isostrain_label)
!* isostrain
call homogenization_isostrain_partitionDeformation(crystallite_partionedF(:,:,:,ip,el), &
crystallite_partionedF0(:,:,:,ip,el),&
materialpoint_subF(:,:,ip,el),&
call homogenization_isostrain_partitionDeformation(crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el), &
crystallite_partionedF0(1:3,1:3,1:homogenization_maxNgrains,ip,el),&
materialpoint_subF(1:3,1:3,ip,el),&
homogenization_state(ip,el), &
ip, &
el)
!* RGC homogenization added <<<updated 31.07.2009>>>
!* RGC homogenization
case (homogenization_RGC_label)
call homogenization_RGC_partitionDeformation(crystallite_partionedF(:,:,:,ip,el), &
crystallite_partionedF0(:,:,:,ip,el),&
materialpoint_subF(:,:,ip,el),&
call homogenization_RGC_partitionDeformation(crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el), &
crystallite_partionedF0(1:3,1:3,1:homogenization_maxNgrains,ip,el),&
materialpoint_subF(1:3,1:3,ip,el),&
homogenization_state(ip,el), &
ip, &
el)
@ -615,7 +621,7 @@ function homogenization_updateState(&
use crystallite, only: crystallite_P,crystallite_dPdF,crystallite_partionedF,crystallite_partionedF0 ! modified <<<updated 31.07.2009>>>
use homogenization_isostrain
use homogenization_RGC ! RGC homogenization added <<<updated 31.07.2009>>>
use homogenization_RGC
implicit none
integer(pInt), intent(in) :: ip,el
@ -624,23 +630,25 @@ function homogenization_updateState(&
select case(homogenization_type(mesh_element(3,el)))
!* isostrain
case (homogenization_isostrain_label)
homogenization_updateState = homogenization_isostrain_updateState( homogenization_state(ip,el), &
crystallite_P(:,:,:,ip,el), &
crystallite_dPdF(:,:,:,:,:,ip,el), &
ip, &
el)
!* RGC homogenization added <<<updated 31.07.2009>>>
homogenization_updateState = &
homogenization_isostrain_updateState( homogenization_state(ip,el), &
crystallite_P(1:3,1:3,1:homogenization_maxNgrains,ip,el), &
crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_maxNgrains,ip,el), &
ip, &
el)
!* RGC homogenization
case (homogenization_RGC_label)
homogenization_updateState = homogenization_RGC_updateState( homogenization_state(ip,el), &
homogenization_subState0(ip,el), &
crystallite_P(:,:,:,ip,el), &
crystallite_partionedF(:,:,:,ip,el), &
crystallite_partionedF0(:,:,:,ip,el),&
materialpoint_subF(:,:,ip,el),&
materialpoint_subdt(ip,el), &
crystallite_dPdF(:,:,:,:,:,ip,el), &
ip, &
el)
homogenization_updateState = &
homogenization_RGC_updateState( homogenization_state(ip,el), &
homogenization_subState0(ip,el), &
crystallite_P(1:3,1:3,1:homogenization_maxNgrains,ip,el), &
crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el), &
crystallite_partionedF0(1:3,1:3,1:homogenization_maxNgrains,ip,el),&
materialpoint_subF(1:3,1:3,ip,el),&
materialpoint_subdt(ip,el), &
crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_maxNgrains,ip,el), &
ip, &
el)
end select
return
@ -660,7 +668,7 @@ subroutine homogenization_averageStressAndItsTangent(&
use material, only: homogenization_type, homogenization_maxNgrains
use crystallite, only: crystallite_P,crystallite_dPdF
use homogenization_RGC ! RGC homogenization added <<<updated 31.07.2009>>>
use homogenization_RGC
use homogenization_isostrain
implicit none
@ -669,18 +677,18 @@ subroutine homogenization_averageStressAndItsTangent(&
select case(homogenization_type(mesh_element(3,el)))
!* isostrain
case (homogenization_isostrain_label)
call homogenization_isostrain_averageStressAndItsTangent( materialpoint_P(:,:,ip,el), &
materialpoint_dPdF(:,:,:,:,ip,el),&
crystallite_P(:,:,:,ip,el), &
crystallite_dPdF(:,:,:,:,:,ip,el), &
ip, &
el)
!* RGC homogenization added <<<updated 31.07.2009>>>
call homogenization_isostrain_averageStressAndItsTangent(materialpoint_P(1:3,1:3,ip,el), &
materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),&
crystallite_P(1:3,1:3,1:homogenization_maxNgrains,ip,el), &
crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_maxNgrains,ip,el), &
ip, &
el)
!* RGC homogenization
case (homogenization_RGC_label)
call homogenization_RGC_averageStressAndItsTangent( materialpoint_P(:,:,ip,el), &
materialpoint_dPdF(:,:,:,:,ip,el),&
crystallite_P(:,:,:,ip,el), &
crystallite_dPdF(:,:,:,:,:,ip,el), &
call homogenization_RGC_averageStressAndItsTangent( materialpoint_P(1:3,1:3,ip,el), &
materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),&
crystallite_P(1:3,1:3,1:homogenization_maxNgrains,ip,el), &
crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_maxNgrains,ip,el), &
ip, &
el)
end select
@ -703,7 +711,7 @@ subroutine homogenization_averageTemperature(&
use crystallite, only: crystallite_Temperature
use homogenization_isostrain
use homogenization_RGC ! RGC homogenization added <<<updated 31.07.2009>>>
use homogenization_RGC
implicit none
integer(pInt), intent(in) :: ip,el
@ -711,10 +719,12 @@ subroutine homogenization_averageTemperature(&
select case(homogenization_type(mesh_element(3,el)))
!* isostrain
case (homogenization_isostrain_label)
materialpoint_Temperature(ip,el) = homogenization_isostrain_averageTemperature(crystallite_Temperature(:,ip,el), ip, el)
!* RGC homogenization added <<<updated 31.07.2009>>>
materialpoint_Temperature(ip,el) = &
homogenization_isostrain_averageTemperature(crystallite_Temperature(1:homogenization_maxNgrains,ip,el), ip, el)
!* RGC homogenization
case (homogenization_RGC_label)
materialpoint_Temperature(ip,el) = homogenization_RGC_averageTemperature(crystallite_Temperature(:,ip,el), ip, el)
materialpoint_Temperature(ip,el) = &
homogenization_RGC_averageTemperature(crystallite_Temperature(1:homogenization_maxNgrains,ip,el), ip, el)
end select
return
@ -734,7 +744,7 @@ function homogenization_postResults(&
use mesh, only: mesh_element
use material, only: homogenization_type
use homogenization_isostrain
use homogenization_RGC ! RGC homogenization added <<<updated 31.07.2009>>>
use homogenization_RGC
implicit none
!* Definition of variables
@ -746,7 +756,7 @@ function homogenization_postResults(&
!* isostrain
case (homogenization_isostrain_label)
homogenization_postResults = homogenization_isostrain_postResults(homogenization_state(ip,el),ip,el)
!* RGC homogenization added <<<updated 31.07.2009>>>
!* RGC homogenization
case (homogenization_RGC_label)
homogenization_postResults = homogenization_RGC_postResults(homogenization_state(ip,el),ip,el)
end select

40
code/homogenization_RGC.f90
View File

@ -59,10 +59,12 @@ subroutine homogenization_RGC_init(&
character(len=64) tag
character(len=1024) line
!$OMP CRITICAL (write2out)
write(6,*)
write(6,'(a21,a20,a12)') '<<<+- homogenization',homogenization_RGC_label,' init -+>>>'
write(6,*) '$Id$'
write(6,*)
!$OMP END CRITICAL (write2out)
maxNinstance = count(homogenization_type == homogenization_RGC_label)
if (maxNinstance == 0) return
@ -148,6 +150,7 @@ subroutine homogenization_RGC_init(&
enddo
100 do i = 1,maxNinstance ! sanity checks
!$OMP CRITICAL (write2out)
write(6,'(a15,x,i4)') 'instance: ', i
write(6,*)
write(6,'(a25,3(x,i8))') 'cluster size: ',(homogenization_RGC_Ngrains(j,i),j=1,3)
@ -155,6 +158,7 @@ subroutine homogenization_RGC_init(&
write(6,'(a25,x,e10.3)') 'over-proportionality: ', homogenization_RGC_ciAlpha(i)
write(6,'(a25,3(x,e10.3))') 'grain size: ',(homogenization_RGC_dAlpha(j,i),j=1,3)
write(6,'(a25,3(x,e10.3))') 'cluster orientation: ',(homogenization_RGC_angles(j,i),j=1,3)
!$OMP END CRITICAL (write2out)
enddo
do i = 1,maxNinstance
@ -256,6 +260,7 @@ subroutine homogenization_RGC_partitionDeformation(&
!* Debugging the overall deformation gradient
if (RGCdebug) then
!$OMP CRITICAL (write2out)
write(6,'(x,a,i3,a,i3,a)')'========== Increment: ',theInc,' Cycle: ',cycleCounter,' =========='
write(6,'(x,a32)')'Overall deformation gradient: '
do i = 1,3
@ -263,6 +268,7 @@ subroutine homogenization_RGC_partitionDeformation(&
enddo
write(6,*)' '
call flush(6)
!$OMP END CRITICAL (write2out)
endif
!* Compute the deformation gradient of individual grains due to relaxations
@ -281,12 +287,14 @@ subroutine homogenization_RGC_partitionDeformation(&
!* Debugging the grain deformation gradients
if (RGCdebug) then
!$OMP CRITICAL (write2out)
write(6,'(x,a32,x,i3)')'Deformation gradient of grain: ',iGrain
do i = 1,3
write(6,'(x,3(e14.8,x))')(F(i,j,iGrain), j = 1,3)
enddo
write(6,*)' '
call flush(6)
!$OMP END CRITICAL (write2out)
endif
enddo
@ -371,11 +379,13 @@ function homogenization_RGC_updateState(&
!* Debugging the obtained state
if (RGCdebug) then
!$OMP CRITICAL (write2out)
write(6,'(x,a30)')'Obtained state: '
do i = 1,3*nIntFaceTot
write(6,'(x,2(e14.8,x))')state%p(i)
enddo
write(6,*)' '
!$OMP END CRITICAL (write2out)
endif
!* Computing interface mismatch and stress penalty tensor for all interfaces of all grains
@ -386,6 +396,7 @@ function homogenization_RGC_updateState(&
!* Debugging the mismatch, stress and penalties of grains
if (RGCdebug) then
!$OMP CRITICAL (write2out)
do iGrain = 1,nGrain
write(6,'(x,a30,x,i3,x,a4,3(x,e14.8))')'Mismatch magnitude of grain(',iGrain,') :',NN(1,iGrain),NN(2,iGrain),NN(3,iGrain)
write(6,*)' '
@ -397,6 +408,7 @@ function homogenization_RGC_updateState(&
enddo
write(6,*)' '
enddo
!$OMP END CRITICAL (write2out)
endif
!* End of initialization
@ -433,9 +445,11 @@ function homogenization_RGC_updateState(&
!* Debugging the residual stress
if (RGCdebug) then
!$OMP CRITICAL (write2out)
write(6,'(x,a30,x,i3)')'Traction at interface: ',iNum
write(6,'(x,3(e14.8,x))')(tract(iNum,j), j = 1,3)
write(6,*)' '
!$OMP END CRITICAL (write2out)
endif
enddo
!* End of residual stress calculation
@ -449,6 +463,7 @@ function homogenization_RGC_updateState(&
!* Debugging the convergent criteria
if (RGCcheck) then
!$OMP CRITICAL (write2out)
write(6,'(x,a)')' '
write(6,'(x,a,x,i2,x,i4)')'RGC residual check ...',ip,el
write(6,'(x,a15,x,e14.8,x,a7,i3,x,a12,i2,i2)')'Max stress: ',stresMax, &
@ -456,6 +471,7 @@ function homogenization_RGC_updateState(&
write(6,'(x,a15,x,e14.8,x,a7,i3,x,a12,i2)')'Max residual: ',residMax, &
'@ iface',residLoc(1),'in direction',residLoc(2)
call flush(6)
!$OMP END CRITICAL (write2out)
endif
homogenization_RGC_updateState = .false.
@ -464,9 +480,11 @@ function homogenization_RGC_updateState(&
homogenization_RGC_updateState = .true.
if (RGCcheck) then
!$OMP CRITICAL (write2out)
write(6,'(x,a55)')'... done and happy'
write(6,*)' '
call flush(6)
!$OMP END CRITICAL (write2out)
endif
! write(6,'(x,a,x,i3,x,a6,x,i3,x,a12)')'RGC_updateState: ip',ip,'| el',el,'converged :)'
@ -492,6 +510,7 @@ function homogenization_RGC_updateState(&
state%p(3*nIntFaceTot+8) = maxval(abs(drelax))/dt ! the maximum rate of relaxation vectors
if (RGCcheck) then
!$OMP CRITICAL (write2out)
write(6,'(x,a30,x,e14.8)')'Constitutive work: ',constitutiveWork
write(6,'(x,a30,3(x,e14.8))')'Magnitude mismatch: ',sum(NN(1,:))/dble(nGrain), &
sum(NN(2,:))/dble(nGrain), &
@ -503,6 +522,7 @@ function homogenization_RGC_updateState(&
write(6,'(x,a30,x,e14.8)')'Average relaxation rate: ',sum(abs(drelax))/dt/dble(3*nIntFaceTot)
write(6,*)''
call flush(6)
!$OMP END CRITICAL (write2out)
endif
deallocate(tract,resid,relax,drelax)
@ -514,9 +534,11 @@ function homogenization_RGC_updateState(&
homogenization_RGC_updateState = (/.true.,.false./) ! with direct cut-back
if (RGCcheck) then
!$OMP CRITICAL (write2out)
write(6,'(x,a55)')'... broken'
write(6,*)' '
call flush(6)
!$OMP END CRITICAL (write2out)
endif
deallocate(tract,resid,relax,drelax)
@ -526,9 +548,11 @@ function homogenization_RGC_updateState(&
else
if (RGCcheck) then
!$OMP CRITICAL (write2out)
write(6,'(x,a55)')'... not yet done'
write(6,*)' '
call flush(6)
!$OMP END CRITICAL (write2out)
endif
endif
@ -578,12 +602,14 @@ function homogenization_RGC_updateState(&
!* Debugging the global Jacobian matrix of stress tangent
if (RGCdebugJacobi) then
!$OMP CRITICAL (write2out)
write(6,'(x,a30)')'Jacobian matrix of stress'
do i = 1,3*nIntFaceTot
write(6,'(x,100(e10.4,x))')(smatrix(i,j), j = 1,3*nIntFaceTot)
enddo
write(6,*)' '
call flush(6)
!$OMP END CRITICAL (write2out)
endif
!* ... of the stress penalty tangent (mismatch penalty and volume penalty,
@ -632,12 +658,14 @@ function homogenization_RGC_updateState(&
!* Debugging the global Jacobian matrix of penalty tangent
if (RGCdebugJacobi) then
!$OMP CRITICAL (write2out)
write(6,'(x,a30)')'Jacobian matrix of penalty'
do i = 1,3*nIntFaceTot
write(6,'(x,100(e10.4,x))')(pmatrix(i,j), j = 1,3*nIntFaceTot)
enddo
write(6,*)' '
call flush(6)
!$OMP END CRITICAL (write2out)
endif
!* ... of the numerical viscosity traction "rmatrix"
@ -650,24 +678,28 @@ function homogenization_RGC_updateState(&
!* Debugging the global Jacobian matrix of numerical viscosity tangent
if (RGCdebugJacobi) then
!$OMP CRITICAL (write2out)
write(6,'(x,a30)')'Jacobian matrix of penalty'
do i = 1,3*nIntFaceTot
write(6,'(x,100(e10.4,x))')(rmatrix(i,j), j = 1,3*nIntFaceTot)
enddo
write(6,*)' '
call flush(6)
!$OMP END CRITICAL (write2out)
endif
!* The overall Jacobian matrix summarizing contributions of smatrix, pmatrix, rmatrix
allocate(jmatrix(3*nIntFaceTot,3*nIntFaceTot)); jmatrix = smatrix + pmatrix + rmatrix
if (RGCdebugJacobi) then
!$OMP CRITICAL (write2out)
write(6,'(x,a30)')'Jacobian matrix (total)'
do i = 1,3*nIntFaceTot
write(6,'(x,100(e10.4,x))')(jmatrix(i,j), j = 1,3*nIntFaceTot)
enddo
write(6,*)' '
call flush(6)
!$OMP END CRITICAL (write2out)
endif
!*** End of construction and assembly of Jacobian matrix
@ -679,12 +711,14 @@ function homogenization_RGC_updateState(&
!* Debugging the inverse Jacobian matrix
if (RGCdebugJacobi) then
!$OMP CRITICAL (write2out)
write(6,'(x,a30)')'Jacobian inverse'
do i = 1,3*nIntFaceTot
write(6,'(x,100(e10.4,x))')(jnverse(i,j), j = 1,3*nIntFaceTot)
enddo
write(6,*)' '
call flush(6)
!$OMP END CRITICAL (write2out)
endif
!* Calculate the state update (global relaxation vectors) for the next Newton-Raphson iteration
@ -698,19 +732,23 @@ function homogenization_RGC_updateState(&
state%p(1:3*nIntFaceTot) = relax
if (any(abs(drelax(:)) > maxdRelax_RGC)) then ! Forcing cutback when the incremental change of relaxation vector becomes too large
homogenization_RGC_updateState = (/.true.,.false./)
!$OMP CRITICAL (write2out)
write(6,'(x,a,x,i3,x,a,x,i3,x,a)')'RGC_updateState: ip',ip,'| el',el,'enforces cutback'
write(6,'(x,a,x,e14.8)')'due to large relaxation change =',maxval(abs(drelax))
call flush(6)
!$OMP END CRITICAL (write2out)
endif
!* Debugging the return state
if (RGCdebugJacobi) then
!$OMP CRITICAL (write2out)
write(6,'(x,a30)')'Returned state: '
do i = 1,3*nIntFaceTot
write(6,'(x,2(e14.8,x))')state%p(i)
enddo
write(6,*)' '
call flush(6)
!$OMP END CRITICAL (write2out)
endif
deallocate(tract,resid,jmatrix,jnverse,relax,drelax,pmatrix,smatrix,p_relax,p_resid)
@ -757,6 +795,7 @@ subroutine homogenization_RGC_averageStressAndItsTangent(&
!* Debugging the grain tangent
if (RGCdebug) then
!$OMP CRITICAL (write2out)
do iGrain = 1,Ngrains
dPdF99 = math_Plain3333to99(dPdF(:,:,:,:,iGrain))
write(6,'(x,a30,x,i3)')'Stress tangent of grain: ',iGrain
@ -766,6 +805,7 @@ subroutine homogenization_RGC_averageStressAndItsTangent(&
write(6,*)' '
enddo
call flush(6)
!$OMP END CRITICAL (write2out)
endif
!* Computing the average first Piola-Kirchhoff stress P and the average tangent dPdF

18
code/homogenization_isostrain.f90
View File

@ -1,13 +1,13 @@
!* $Id$
!*****************************************************
!* Module: HOMOGENIZATION_ISOSTRAIN *
!* Module: HOMOGENIZATION_ISOSTRAIN *
!*****************************************************
!* contains: *
!*****************************************************
! [isostrain]
! type isostrain
! Ngrains 6
! [isostrain]
! type isostrain
! Ngrains 6
! (output) Ngrains
MODULE homogenization_isostrain
@ -115,9 +115,9 @@ subroutine homogenization_isostrain_init(&
end select
if (mySize > 0_pInt) then ! any meaningful output found
homogenization_isostrain_sizePostResult(j,i) = mySize
homogenization_isostrain_sizePostResults(i) = &
homogenization_isostrain_sizePostResults(i) + mySize
homogenization_isostrain_sizePostResult(j,i) = mySize
homogenization_isostrain_sizePostResults(i) = &
homogenization_isostrain_sizePostResults(i) + mySize
endif
enddo
enddo
@ -137,7 +137,7 @@ function homogenization_isostrain_stateInit(myInstance)
!* Definition of variables
integer(pInt), intent(in) :: myInstance
real(pReal), dimension(homogenization_isostrain_sizeState(myInstance)) :: &
homogenization_isostrain_stateInit ! modified <<<updated 31.07.2009>>>
homogenization_isostrain_stateInit
homogenization_isostrain_stateInit = 0.0_pReal
@ -173,7 +173,7 @@ subroutine homogenization_isostrain_partitionDeformation(&
! homID = homogenization_typeInstance(mesh_element(3,el))
forall (i = 1:homogenization_Ngrains(mesh_element(3,el))) &
F(:,:,i) = avgF
F(1:3,1:3,i) = avgF
return

4
code/lattice.f90
View File

@ -694,10 +694,12 @@ subroutine lattice_init()
integer(pInt), parameter :: fileunit = 200
integer(pInt) i,Nsections
!$OMP CRITICAL (write2out)
write(6,*)
write(6,*) '<<<+- lattice init -+>>>'
write(6,*) '$Id$'
write(6,*)
!$OMP END CRITICAL (write2out)
if(.not. IO_open_file(fileunit,material_configFile)) call IO_error(100) ! cannot open config file
Nsections = IO_countSections(fileunit,material_partPhase)
@ -705,9 +707,11 @@ subroutine lattice_init()
! lattice_Nstructure = Nsections + 2_pInt ! most conservative assumption
close(fileunit)
!$OMP CRITICAL (write2out)
write(6,'(a16,x,i5)') '# phases:',Nsections
write(6,'(a16,x,i5)') '# structures:',lattice_Nstructure
write(6,*)
!$OMP END CRITICAL (write2out)
allocate(lattice_Sslip(3,3,lattice_maxNslip,lattice_Nstructure)); lattice_Sslip = 0.0_pReal
allocate(lattice_Sslip_v(6,lattice_maxNslip,lattice_Nstructure)); lattice_Sslip_v = 0.0_pReal

13
code/material.f90
View File

@ -93,11 +93,13 @@ subroutine material_init()
integer(pInt), parameter :: fileunit = 200
integer(pInt) i,j
!$OMP CRITICAL (write2out)
write(6,*)
write(6,*) '<<<+- material init -+>>>'
write(6,*) '$Id$'
write(6,*)
!$OMP END CRITICAL (write2out)
if(.not. IO_open_file(fileunit,material_configFile)) call IO_error(100) ! cannot open config file
call material_parseHomogenization(fileunit,material_partHomogenization)
call material_parseMicrostructure(fileunit,material_partMicrostructure)
@ -106,6 +108,7 @@ subroutine material_init()
call material_parsePhase(fileunit,material_partPhase)
close(fileunit)
!$OMP CRITICAL (write2out)
do i = 1,material_Nmicrostructure
if (microstructure_crystallite(i) < 1 .or. &
microstructure_crystallite(i) > material_Ncrystallite) call IO_error(150,i)
@ -141,7 +144,8 @@ subroutine material_init()
write (6,*)
endif
enddo
!$OMP END CRITICAL (write2out)
call material_populateGrains()
endsubroutine
@ -577,17 +581,21 @@ subroutine material_populateGrains()
allocate(phaseOfGrain(maxval(Ngrains))) ! reserve memory for maximum case
allocate(orientationOfGrain(3,maxval(Ngrains))) ! reserve memory for maximum case
!$OMP CRITICAL (write2out)
write (6,*)
write (6,*) 'MATERIAL grain population'
write (6,*)
write (6,'(a32,x,a32,x,a6)') 'homogenization_name','microstructure_name','grain#'
!$OMP END CRITICAL (write2out)
do homog = 1,material_Nhomogenization ! loop over homogenizations
dGrains = homogenization_Ngrains(homog) ! grain number per material point
do micro = 1,material_Nmicrostructure ! all pairs of homog and micro
if (Ngrains(homog,micro) > 0) then ! an active pair of homog and micro
myNgrains = Ngrains(homog,micro) ! assign short name for total number of grains to populate
!$OMP CRITICAL (write2out)
write (6,*)
write (6,'(a32,x,a32,x,i6)') homogenization_name(homog),microstructure_name(micro),myNgrains
!$OMP END CRITICAL (write2out)
! ---------------------------------------------------------------------------- calculate volume of each grain
volumeOfGrain = 0.0_pReal
@ -729,7 +737,6 @@ subroutine material_populateGrains()
endif ! active homog,micro pair
enddo
enddo
deallocate(volumeOfGrain)
deallocate(phaseOfGrain)

9
code/math.f90
View File

@ -130,10 +130,12 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = &
integer(pInt), dimension(1) :: randInit
!$OMP CRITICAL (write2out)
write(6,*)
write(6,*) '<<<+- math init -+>>>'
write(6,*) '$Id$'
write(6,*)
!$OMP END CRITICAL (write2out)
if (fixedSeed > 0_pInt) then
randInit = fixedSeed
@ -143,8 +145,10 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = &
endif
call random_seed(get=randInit)
!$OMP CRITICAL (write2out)
write(6,*) 'random seed: ',randInit(1)
write(6,*)
!$OMP END CRITICAL (write2out)
call halton_seed_set(randInit(1))
call halton_ndim_set(3)
@ -2576,7 +2580,7 @@ endfunction
r(1:ndim) = 0.0_pReal
if ( any ( base(1:ndim) <= 1 ) ) then
!$OMP CRITICAL (write2out)
!$OMP CRITICAL (write2out)
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'I_TO_HALTON - Fatal error!'
write ( *, '(a)' ) ' An input base BASE is <= 1!'
@ -2584,7 +2588,7 @@ endfunction
write ( *, '(i6,i6)' ) i, base(i)
enddo
call flush(6)
!$OMP END CRITICAL (write2out)
!$OMP END CRITICAL (write2out)
stop
end if
@ -2855,7 +2859,6 @@ endfunction
else
prime = 0
!$OMP CRITICAL (write2out)
write ( 6, '(a)' ) ' '
write ( 6, '(a)' ) 'PRIME - Fatal error!'
write ( 6, '(a,i6)' ) ' Illegal prime index N = ', n

4
code/mesh.f90
View File

@ -239,10 +239,12 @@
integer(pInt), parameter :: fileUnit = 222
integer(pInt) e,element,ip
!$OMP CRITICAL (write2out)
write(6,*)
write(6,*) '<<<+- mesh init -+>>>'
write(6,*) '$Id$'
write(6,*)
!$OMP END CRITICAL (write2out)
call mesh_build_FEdata() ! --- get properties of the different types of elements
@ -304,7 +306,9 @@
forall (e = 1:mesh_NcpElems) FEsolving_execIP(2,e) = FE_Nips(mesh_element(2,e))
allocate(calcMode(mesh_maxNips,mesh_NcpElems))
!$OMP CRITICAL (write2out)
write(6,*) '<<<+- mesh init done -+>>>'
!$OMP END CRITICAL (write2out)
calcMode = .false. ! pretend to have collected what first call is asking (F = I)
calcMode(ip,mesh_FEasCP('elem',element)) = .true. ! first ip,el needs to be already pingponged to "calc"
lastMode = .true. ! and its mode is already known...

15
code/mpie_cpfem_marc.f90
View File

@ -247,14 +247,15 @@ subroutine hypela2(&
d_max, d_min
! OpenMP variable
!$ integer(pInt) defaultNumThreadsInt ! default value set by Marc
!$ integer(pInt) defaultNumThreadsInt ! default value set by Marc
!$ defaultNumThreadsInt = omp_get_num_threads() ! remember number of threads 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)
cp_en = mesh_FEasCP('elem',n(1))
!$ call omp_set_num_threads(mpieNumThreadsInt) ! set number of threads for parallel execution set by MPIE_NUM_THREADS
!$ call omp_set_num_threads(mpieNumThreadsInt) ! set number of threads for parallel execution set by MPIE_NUM_THREADS
if (lovl == 4) then ! Marc requires stiffness in separate call
if ( timinc < theDelta .and. theInc == inc ) then ! first after cutback
@ -263,10 +264,7 @@ subroutine hypela2(&
computationMode = 6 ! --> just return known tangent
endif
else ! stress requested (lovl == 6)
cp_en = mesh_FEasCP('elem',n(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 <<
@ -294,16 +292,13 @@ subroutine hypela2(&
write (6,'(i6,x,i2,x,a)') n(1),nn,'<< hypela2 >> new increment..!'; call flush(6)
!$OMP END CRITICAL (write2out)
endif
else if ( timinc < theDelta ) then ! >> cutBack <<
terminallyIll = .false.
cycleCounter = -1 ! first calc step increments this to cycle = 0
calcMode = .true. ! pretend last step was calculation
!$OMP CRITICAL (write2out)
write(6,'(i6,x,i2,x,a)') n(1),nn,'<< hypela2 >> cutback detected..!'; call flush(6)
!$OMP END CRITICAL (write2out)
endif ! convergence treatment end
calcMode(nn,cp_en) = .not. calcMode(nn,cp_en) ! ping pong (calc <--> collect)
@ -328,6 +323,7 @@ subroutine hypela2(&
if ( lastMode /= calcMode(nn,cp_en) .and. &
.not. terminallyIll ) then
call debug_info() ! first after ping pong reports (meaningful) debugging
!$OMP CRITICAL (write2out)
write(6,*)
write(6,*) 'EXTREME VALUES OF RETURNED VARIABLES (from previous cycle)'
write(6,*)
@ -337,6 +333,7 @@ subroutine hypela2(&
write(6,'(a14,x,e12.3,x,i6,x,i4)') 'jacobian min :', d_min, d_min_e, d_min_i
write(6,'(a14,x,e12.3,x,i6,x,i4)') ' max :', d_max, d_max_e, d_max_i
write(6,*)
!$OMP END CRITICAL (write2out)
endif
if ( lastIncConverged ) then
lastIncConverged = .false. ! reset flag

122
code/numerics.f90
View File

@ -96,10 +96,12 @@ subroutine numerics_init()
! OpenMP variable
!$ character(len=4) mpieNumThreadsString !environment variable MPIE_NUMTHREADS
write(6,*)
write(6,*) '<<<+- numerics init -+>>>'
write(6,*) '$Id$'
write(6,*)
!$OMP CRITICAL (write2out)
write(6,*)
write(6,*) '<<<+- numerics init -+>>>'
write(6,*) '$Id$'
write(6,*)
!$OMP END CRITICAL (write2out)
! initialize all parameters with standard values
relevantStrain = 1.0e-7_pReal
@ -161,8 +163,10 @@ subroutine numerics_init()
! try to open the config file
if(IO_open_file(fileunit,numerics_configFile)) then
write(6,*) ' ... using values from config file'
write(6,*)
!$OMP CRITICAL (write2out)
write(6,*) ' ... using values from config file'
write(6,*)
!$OMP END CRITICAL (write2out)
line = ''
! read variables from config file and overwrite parameters
@ -269,64 +273,68 @@ subroutine numerics_init()
! no config file, so we use standard values
else
write(6,*) ' ... using standard values'
write(6,*)
!$OMP CRITICAL (write2out)
write(6,*) ' ... using standard values'
write(6,*)
!$OMP END CRITICAL (write2out)
endif
! writing parameters to output file
write(6,'(a24,x,e8.1)') 'relevantStrain: ',relevantStrain
write(6,'(a24,x,e8.1)') 'defgradTolerance: ',defgradTolerance
write(6,'(a24,x,i8)') 'iJacoStiffness: ',iJacoStiffness
write(6,'(a24,x,i8)') 'iJacoLpresiduum: ',iJacoLpresiduum
write(6,'(a24,x,e8.1)') 'pert_Fg: ',pert_Fg
write(6,'(a24,x,i8)') 'pert_method: ',pert_method
write(6,'(a24,x,i8)') 'nCryst: ',nCryst
write(6,'(a24,x,e8.1)') 'subStepMinCryst: ',subStepMinCryst
write(6,'(a24,x,e8.1)') 'subStepSizeCryst: ',subStepSizeCryst
write(6,'(a24,x,e8.1)') 'stepIncreaseCryst: ',stepIncreaseCryst
write(6,'(a24,x,i8)') 'nState: ',nState
write(6,'(a24,x,i8)') 'nStress: ',nStress
write(6,'(a24,x,e8.1)') 'rTol_crystalliteState: ',rTol_crystalliteState
write(6,'(a24,x,e8.1)') 'rTol_crystalliteTemp: ',rTol_crystalliteTemperature
write(6,'(a24,x,e8.1)') 'rTol_crystalliteStress: ',rTol_crystalliteStress
write(6,'(a24,x,e8.1)') 'aTol_crystalliteStress: ',aTol_crystalliteStress
write(6,'(a24,2(x,i8))')'integrator: ',numerics_integrator
write(6,*)
write(6,'(a24,x,i8)') 'nHomog: ',nHomog
write(6,'(a24,x,e8.1)') 'subStepMinHomog: ',subStepMinHomog
write(6,'(a24,x,e8.1)') 'subStepSizeHomog: ',subStepSizeHomog
write(6,'(a24,x,e8.1)') 'stepIncreaseHomog: ',stepIncreaseHomog
write(6,'(a24,x,i8)') 'nMPstate: ',nMPstate
write(6,*)
!$OMP CRITICAL (write2out)
write(6,'(a24,x,e8.1)') 'relevantStrain: ',relevantStrain
write(6,'(a24,x,e8.1)') 'defgradTolerance: ',defgradTolerance
write(6,'(a24,x,i8)') 'iJacoStiffness: ',iJacoStiffness
write(6,'(a24,x,i8)') 'iJacoLpresiduum: ',iJacoLpresiduum
write(6,'(a24,x,e8.1)') 'pert_Fg: ',pert_Fg
write(6,'(a24,x,i8)') 'pert_method: ',pert_method
write(6,'(a24,x,i8)') 'nCryst: ',nCryst
write(6,'(a24,x,e8.1)') 'subStepMinCryst: ',subStepMinCryst
write(6,'(a24,x,e8.1)') 'subStepSizeCryst: ',subStepSizeCryst
write(6,'(a24,x,e8.1)') 'stepIncreaseCryst: ',stepIncreaseCryst
write(6,'(a24,x,i8)') 'nState: ',nState
write(6,'(a24,x,i8)') 'nStress: ',nStress
write(6,'(a24,x,e8.1)') 'rTol_crystalliteState: ',rTol_crystalliteState
write(6,'(a24,x,e8.1)') 'rTol_crystalliteTemp: ',rTol_crystalliteTemperature
write(6,'(a24,x,e8.1)') 'rTol_crystalliteStress: ',rTol_crystalliteStress
write(6,'(a24,x,e8.1)') 'aTol_crystalliteStress: ',aTol_crystalliteStress
write(6,'(a24,2(x,i8))')'integrator: ',numerics_integrator
write(6,*)
write(6,'(a24,x,i8)') 'nHomog: ',nHomog
write(6,'(a24,x,e8.1)') 'subStepMinHomog: ',subStepMinHomog
write(6,'(a24,x,e8.1)') 'subStepSizeHomog: ',subStepSizeHomog
write(6,'(a24,x,e8.1)') 'stepIncreaseHomog: ',stepIncreaseHomog
write(6,'(a24,x,i8)') 'nMPstate: ',nMPstate
write(6,*)
!* RGC parameters
write(6,'(a24,x,e8.1)') 'aTol_RGC: ',absTol_RGC
write(6,'(a24,x,e8.1)') 'rTol_RGC: ',relTol_RGC
write(6,'(a24,x,e8.1)') 'aMax_RGC: ',absMax_RGC
write(6,'(a24,x,e8.1)') 'rMax_RGC: ',relMax_RGC
write(6,'(a24,x,e8.1)') 'perturbPenalty_RGC: ',pPert_RGC
write(6,'(a24,x,e8.1)') 'relevantMismatch_RGC: ',xSmoo_RGC
write(6,'(a24,x,e8.1)') 'viscosityrate_RGC: ',viscPower_RGC
write(6,'(a24,x,e8.1)') 'viscositymodulus_RGC: ',viscModus_RGC
write(6,'(a24,x,e8.1)') 'maxrelaxation_RGC: ',maxdRelax_RGC
write(6,'(a24,x,e8.1)') 'maxVolDiscrepancy_RGC: ',maxVolDiscr_RGC
write(6,'(a24,x,e8.1)') 'volDiscrepancyMod_RGC: ',volDiscrMod_RGC
write(6,'(a24,x,e8.1)') 'discrepancyPower_RGC: ',volDiscrPow_RGC
write(6,*)
write(6,'(a24,x,e8.1)') 'aTol_RGC: ',absTol_RGC
write(6,'(a24,x,e8.1)') 'rTol_RGC: ',relTol_RGC
write(6,'(a24,x,e8.1)') 'aMax_RGC: ',absMax_RGC
write(6,'(a24,x,e8.1)') 'rMax_RGC: ',relMax_RGC
write(6,'(a24,x,e8.1)') 'perturbPenalty_RGC: ',pPert_RGC
write(6,'(a24,x,e8.1)') 'relevantMismatch_RGC: ',xSmoo_RGC
write(6,'(a24,x,e8.1)') 'viscosityrate_RGC: ',viscPower_RGC
write(6,'(a24,x,e8.1)') 'viscositymodulus_RGC: ',viscModus_RGC
write(6,'(a24,x,e8.1)') 'maxrelaxation_RGC: ',maxdRelax_RGC
write(6,'(a24,x,e8.1)') 'maxVolDiscrepancy_RGC: ',maxVolDiscr_RGC
write(6,'(a24,x,e8.1)') 'volDiscrepancyMod_RGC: ',volDiscrMod_RGC
write(6,'(a24,x,e8.1)') 'discrepancyPower_RGC: ',volDiscrPow_RGC
write(6,*)
!* spectral parameters
write(6,'(a24,x,e8.1)') 'err_div_tol: ',err_div_tol
write(6,'(a24,x,e8.1)') 'err_defgrad_tol: ',err_defgrad_tol
write(6,'(a24,x,e8.1)') 'err_stress_tolrel: ',err_stress_tolrel
write(6,'(a24,x,i8)') 'itmax: ',itmax
write(6,'(a24,x,L8)') 'memory_efficient: ',memory_efficient
write(6,*)
write(6,'(a24,x,e8.1)') 'err_div_tol: ',err_div_tol
write(6,'(a24,x,e8.1)') 'err_defgrad_tol: ',err_defgrad_tol
write(6,'(a24,x,e8.1)') 'err_stress_tolrel: ',err_stress_tolrel
write(6,'(a24,x,i8)') 'itmax: ',itmax
write(6,'(a24,x,L8)') 'memory_efficient: ',memory_efficient
write(6,*)
!* Random seeding parameters
write(6,'(a24,x,i8)') 'fixed_seed: ',fixedSeed
write(6,*)
write(6,'(a24,x,i8)') 'fixed_seed: ',fixedSeed
write(6,*)
!$OMP END CRITICAL (write2out)
!* openMP parameter
!$ write(6,'(a24,x,i8)') 'number of threads: ',OMP_get_max_threads()
@ -379,7 +387,11 @@ subroutine numerics_init()
if (err_stress_tolrel <= 0.0_pReal) call IO_error(49)
if (itmax <= 1.0_pInt) call IO_error(49)
if (fixedSeed <= 0_pInt) write(6,'(a)') 'Random is random!'
if (fixedSeed <= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(a)') 'Random is random!'
!$OMP END CRITICAL (write2out)
endif
endsubroutine
END MODULE numerics

35
code/prec.f90
View File

@ -3,30 +3,31 @@
MODULE prec
!##############################################################
implicit none
implicit none
! *** Precision of real and integer variables ***
integer, parameter :: pReal = selected_real_kind(15,300) ! 15 significant digits, up to 1e+-300
integer, parameter :: pInt = selected_int_kind(9) ! up to +- 1e9
integer, parameter :: pLongInt = 8 ! should be 64bit
real(pReal), parameter :: tol_math_check = 1.0e-8_pReal
real(pReal), parameter :: tol_gravityNodePos = 1.0e-100_pReal
integer, parameter :: pReal = selected_real_kind(15,300) ! 15 significant digits, up to 1e+-300
integer, parameter :: pInt = selected_int_kind(9) ! up to +- 1e9
integer, parameter :: pLongInt = 8 ! should be 64bit
real(pReal), parameter :: tol_math_check = 1.0e-8_pReal
real(pReal), parameter :: tol_gravityNodePos = 1.0e-100_pReal
type :: p_vec
real(pReal), dimension(:), pointer :: p
end type p_vec
type :: p_vec
real(pReal), dimension(:), pointer :: p
end type p_vec
CONTAINS
subroutine prec_init
implicit none
implicit none
!$OMP CRITICAL (write2out)
write(6,*)
write(6,*) '<<<+- prec init -+>>>'
write(6,*) '$Id$'
write(6,*)
!$OMP END CRITICAL (write2out)
write(6,*)
write(6,*) '<<<+- prec init -+>>>'
write(6,*) '$Id$'
write(6,*)
return
end subroutine
END MODULE prec
END MODULE prec

33
code/prec_single.f90
View File

@ -3,28 +3,29 @@
MODULE prec
!##############################################################
implicit none
implicit none
! *** Precision of real and integer variables ***
integer, parameter :: pReal = selected_real_kind(6,37) ! 6 significant digits, up to 1e+-37
integer, parameter :: pInt = selected_int_kind(9) ! up to +- 1e9
integer, parameter :: pLongInt = 8 ! should be 64bit
real(pReal), parameter :: tol_math_check = 1.0e-5_pReal
real(pReal), parameter :: tol_gravityNodePos = 1.0e-36_pReal
integer, parameter :: pReal = selected_real_kind(6,37) ! 6 significant digits, up to 1e+-37
integer, parameter :: pInt = selected_int_kind(9) ! up to +- 1e9
integer, parameter :: pLongInt = 8 ! should be 64bit
real(pReal), parameter :: tol_math_check = 1.0e-5_pReal
real(pReal), parameter :: tol_gravityNodePos = 1.0e-36_pReal
type :: p_vec
real(pReal), dimension(:), pointer :: p
end type p_vec
type :: p_vec
real(pReal), dimension(:), pointer :: p
end type p_vec
CONTAINS
subroutine prec_init
write(6,*)
write(6,*) '<<<+- prec init -+>>>'
write(6,*) '$Id: prec.f90 407 2009-08-31 15:09:15Z MPIE\f.roters $'
write(6,*)
return
!$OMP CRITICAL (write2out)
write(6,*)
write(6,*) '<<<+- prec init -+>>>'
write(6,*) '$Id: prec.f90 407 2009-08-31 15:09:15Z MPIE\f.roters $'
write(6,*)
!$OMP END CRITICAL (write2out)
end subroutine
END MODULE prec
END MODULE prec