From 235266b16927c75256fd9f918c122a0e9a93f890 Mon Sep 17 00:00:00 2001 From: Christoph Kords Date: Thu, 17 Mar 2011 10:46:17 +0000 Subject: [PATCH] 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) !$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. --- code/CPFEM.f90 | 72 +- code/IO.f90 | 7 +- code/constitutive.f90 | 336 ++++--- code/constitutive_j2.f90 | 10 +- code/constitutive_nonlocal.f90 | 118 ++- code/constitutive_phenopowerlaw.f90 | 35 +- code/crystallite.f90 | 1434 +++++++++++++++------------ code/debug.f90 | 40 +- code/homogenization.f90 | 166 ++-- code/homogenization_RGC.f90 | 40 + code/homogenization_isostrain.f90 | 18 +- code/lattice.f90 | 4 + code/material.f90 | 13 +- code/math.f90 | 9 +- code/mesh.f90 | 4 + code/mpie_cpfem_marc.f90 | 15 +- code/numerics.f90 | 122 ++- code/prec.f90 | 35 +- code/prec_single.f90 | 33 +- 19 files changed, 1427 insertions(+), 1084 deletions(-) diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index 8f7a94476..b13e5c6b9 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -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 diff --git a/code/IO.f90 b/code/IO.f90 index 0e768bf02..51c7b2988 100644 --- a/code/IO.f90 +++ b/code/IO.f90 @@ -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) diff --git a/code/constitutive.f90 b/code/constitutive.f90 index cb07bdf07..106a6ecbb 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -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 \ No newline at end of file diff --git a/code/constitutive_j2.f90 b/code/constitutive_j2.f90 index 34b33a6ba..86adfca64 100644 --- a/code/constitutive_j2.f90 +++ b/code/constitutive_j2.f90 @@ -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 diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index ee8fd504b..60782dbd0 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -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 diff --git a/code/constitutive_phenopowerlaw.f90 b/code/constitutive_phenopowerlaw.f90 index 13a12212d..c5b514312 100644 --- a/code/constitutive_phenopowerlaw.f90 +++ b/code/constitutive_phenopowerlaw.f90 @@ -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 diff --git a/code/crystallite.f90 b/code/crystallite.f90 index e29cf445a..0c22a1a15 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -72,9 +72,7 @@ logical, dimension (:,:,:), allocatable :: & crystallite_localConstitution, & ! indicates this grain to have purely local constitutive law crystallite_requested, & ! flag to request crystallite calculation crystallite_todo, & ! flag to indicate need for further computation - crystallite_converged, & ! convergence flag - crystallite_stateConverged, & ! flag indicating convergence of state - crystallite_temperatureConverged ! flag indicating convergence of temperature + crystallite_converged ! convergence flag CONTAINS @@ -153,10 +151,12 @@ character(len=64) tag character(len=1024) line -write(6,*) -write(6,*) '<<<+- crystallite init -+>>>' -write(6,*) '$Id$' -write(6,*) +!$OMP CRITICAL (write2out) + write(6,*) + write(6,*) '<<<+- crystallite init -+>>>' + write(6,*) '$Id$' + write(6,*) +!$OMP END CRITICAL (write2out) gMax = homogenization_maxNgrains @@ -206,8 +206,6 @@ allocate(crystallite_localConstitution(gMax,iMax,eMax)); crystallite_loc allocate(crystallite_requested(gMax,iMax,eMax)); crystallite_requested = .false. allocate(crystallite_todo(gMax,iMax,eMax)); crystallite_todo = .false. allocate(crystallite_converged(gMax,iMax,eMax)); crystallite_converged = .true. -allocate(crystallite_stateConverged(gMax,iMax,eMax)); crystallite_stateConverged = .false. -allocate(crystallite_temperatureConverged(gMax,iMax,eMax)); crystallite_temperatureConverged = .false. allocate(crystallite_output(maxval(crystallite_Noutput), & material_Ncrystallite)) ; crystallite_output = '' @@ -290,26 +288,36 @@ enddo close(file) -!$OMP PARALLEL PRIVATE(myNgrains,myPhase,myStructure) +!$OMP PARALLEL PRIVATE(myNgrains,myPhase,myMat,myStructure) !$OMP DO do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over all cp elements myNgrains = homogenization_Ngrains(mesh_element(3,e)) ! look up homogenization-->grainCount do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element do g = 1,myNgrains - crystallite_partionedTemperature0(g,i,e) = Temperature ! isothermal assumption - crystallite_Fp0(:,:,g,i,e) = math_EulerToR(material_EulerAngles(1:3,g,i,e)) ! plastic def gradient reflects init orientation - crystallite_Fe(:,:,g,i,e) = math_transpose3x3(crystallite_Fp0(1:3,1:3,g,i,e)) - crystallite_F0(:,:,g,i,e) = math_I3 - crystallite_partionedFp0(:,:,g,i,e) = crystallite_Fp0(:,:,g,i,e) - crystallite_partionedF0(:,:,g,i,e) = crystallite_F0(:,:,g,i,e) - crystallite_partionedF(:,:,g,i,e) = crystallite_F0(:,:,g,i,e) - crystallite_requested(g,i,e) = .true. + crystallite_Fp0(1:3,1:3,g,i,e) = math_EulerToR(material_EulerAngles(1:3,g,i,e)) ! plastic def gradient reflects init orientation + crystallite_F0(1:3,1:3,g,i,e) = math_I3 crystallite_localConstitution(g,i,e) = phase_localConstitution(material_phase(g,i,e)) + !$OMP FLUSH(crystallite_Fp0) + crystallite_Fe(1:3,1:3,g,i,e) = math_transpose3x3(crystallite_Fp0(1:3,1:3,g,i,e)) enddo enddo enddo !$OMP ENDDO +!$OMP DO + do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over all cp elements + myNgrains = homogenization_Ngrains(mesh_element(3,e)) ! look up homogenization-->grainCount + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element + do g = 1,myNgrains + enddo + enddo + enddo +!$OMP ENDDO +crystallite_partionedTemperature0 = Temperature ! isothermal assumption +crystallite_partionedFp0 = crystallite_Fp0 +crystallite_partionedF0 = crystallite_F0 +crystallite_partionedF = crystallite_F0 +crystallite_requested = .true. ! Initialize crystallite_symmetryID(g,i,e) @@ -350,6 +358,7 @@ call crystallite_stressAndItsTangent(.true.) ! request elastic crystallite_fallbackdPdF = crystallite_dPdF ! use initial elastic stiffness as fallback ! *** Output to MARC output file *** +!$OMP CRITICAL (write2out) write(6,'(a35,x,7(i5,x))') 'crystallite_Temperature: ', shape(crystallite_Temperature) write(6,'(a35,x,7(i5,x))') 'crystallite_dotTemperature: ', shape(crystallite_dotTemperature) write(6,'(a35,x,7(i5,x))') 'crystallite_Fe: ', shape(crystallite_Fe) @@ -391,19 +400,16 @@ write(6,'(a35,x,7(i5,x))') 'crystallite_localConstitution: ', shape(crystall write(6,'(a35,x,7(i5,x))') 'crystallite_requested: ', shape(crystallite_requested) write(6,'(a35,x,7(i5,x))') 'crystallite_todo: ', shape(crystallite_todo) write(6,'(a35,x,7(i5,x))') 'crystallite_converged: ', shape(crystallite_converged) -write(6,'(a35,x,7(i5,x))') 'crystallite_stateConverged: ', shape(crystallite_stateConverged) -write(6,'(a35,x,7(i5,x))') 'crystallite_temperatureConverged: ', shape(crystallite_temperatureConverged) write(6,'(a35,x,7(i5,x))') 'crystallite_sizePostResults: ', shape(crystallite_sizePostResults) write(6,'(a35,x,7(i5,x))') 'crystallite_sizePostResult: ', shape(crystallite_sizePostResult) write(6,*) write(6,*) 'Number of nonlocal grains: ',count(.not. crystallite_localConstitution) call flush(6) +!$OMP END CRITICAL (write2out) call debug_info() call debug_reset() -return - endsubroutine @@ -426,7 +432,6 @@ use numerics, only: subStepMinCryst, & numerics_integrator, & numerics_integrationMode use debug, only: debugger, & - selectiveDebugger, & verboseDebugger, & debug_e, & debug_i, & @@ -439,6 +444,7 @@ use math, only: math_inv3x3, & math_mul66x6, & math_Mandel6to33, & math_Mandel33to6, & + math_transpose3x3, & math_I3 use FEsolving, only: FEsolving_execElem, & FEsolving_execIP, & @@ -517,12 +523,12 @@ crystallite_subStep = 0.0_pReal do g = 1,myNgrains if (crystallite_requested(g,i,e)) then ! initialize restoration point of ... crystallite_subTemperature0(g,i,e) = crystallite_partionedTemperature0(g,i,e) ! ...temperature - constitutive_subState0(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructure - crystallite_subFp0(:,:,g,i,e) = crystallite_partionedFp0(:,:,g,i,e) ! ...plastic def grad - crystallite_subLp0(:,:,g,i,e) = crystallite_partionedLp0(:,:,g,i,e) ! ...plastic velocity grad - crystallite_dPdF0(:,:,:,:,g,i,e) = crystallite_partioneddPdF0(:,:,:,:,g,i,e) ! ...stiffness - crystallite_subF0(:,:,g,i,e) = crystallite_partionedF0(:,:,g,i,e) ! ...def grad - crystallite_subTstar0_v(:,g,i,e) = crystallite_partionedTstar0_v(:,g,i,e) !...2nd PK stress + constitutive_subState0(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructure + crystallite_subFp0(1:3,1:3,g,i,e) = crystallite_partionedFp0(1:3,1:3,g,i,e) ! ...plastic def grad + crystallite_subLp0(1:3,1:3,g,i,e) = crystallite_partionedLp0(1:3,1:3,g,i,e) ! ...plastic velocity grad + crystallite_dPdF0(1:3,1:3,1:3,1:3,g,i,e) = crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,g,i,e) ! ...stiffness + crystallite_subF0(1:3,1:3,g,i,e) = crystallite_partionedF0(1:3,1:3,g,i,e) ! ...def grad + crystallite_subTstar0_v(1:6,g,i,e) = crystallite_partionedTstar0_v(1:6,g,i,e) !...2nd PK stress crystallite_subFrac(g,i,e) = 0.0_pReal crystallite_subStep(g,i,e) = 1.0_pReal/subStepSizeCryst @@ -560,15 +566,18 @@ do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2 endif crystallite_subFrac(g,i,e) = crystallite_subFrac(g,i,e) + crystallite_subStep(g,i,e) formerSubStep = crystallite_subStep(g,i,e) + !$OMP FLUSH(crystallite_subFrac) crystallite_subStep(g,i,e) = min( 1.0_pReal - crystallite_subFrac(g,i,e), & stepIncreaseCryst * crystallite_subStep(g,i,e) ) + !$OMP FLUSH(crystallite_subStep) if (crystallite_subStep(g,i,e) > subStepMinCryst) then crystallite_subTemperature0(g,i,e) = crystallite_Temperature(g,i,e) ! wind forward... - crystallite_subF0(:,:,g,i,e) = crystallite_subF(:,:,g,i,e) ! ...def grad - crystallite_subFp0(:,:,g,i,e) = crystallite_Fp(:,:,g,i,e) ! ...plastic def grad - crystallite_subLp0(:,:,g,i,e) = crystallite_Lp(:,:,g,i,e) ! ...plastic velocity gradient + crystallite_subF0(1:3,1:3,g,i,e) = crystallite_subF(1:3,1:3,g,i,e) ! ...def grad + crystallite_subFp0(1:3,1:3,g,i,e) = crystallite_Fp(1:3,1:3,g,i,e) ! ...plastic def grad + crystallite_subLp0(1:3,1:3,g,i,e) = crystallite_Lp(1:3,1:3,g,i,e) ! ...plastic velocity gradient constitutive_subState0(g,i,e)%p = constitutive_state(g,i,e)%p ! ...microstructure - crystallite_subTstar0_v(:,g,i,e) = crystallite_Tstar_v(:,g,i,e) ! ...2nd PK stress + crystallite_subTstar0_v(1:6,g,i,e) = crystallite_Tstar_v(1:6,g,i,e) ! ...2nd PK stress + !$OMP FLUSH(crystallite_subF0) elseif (formerSubStep > subStepMinCryst) then ! this crystallite just converged !$OMP CRITICAL (distributionCrystallite) debug_CrystalliteLoopDistribution(min(nCryst+1,NiterationCrystallite)) = & @@ -581,12 +590,13 @@ do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2 else crystallite_subStep(g,i,e) = subStepSizeCryst * crystallite_subStep(g,i,e) ! cut step in half and restore... crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) ! ...temperature - crystallite_Fp(:,:,g,i,e) = crystallite_subFp0(:,:,g,i,e) ! ...plastic def grad - crystallite_invFp(:,:,g,i,e) = math_inv3x3(crystallite_Fp(1:3,1:3,g,i,e)) - crystallite_Lp(:,:,g,i,e) = crystallite_subLp0(:,:,g,i,e) ! ...plastic velocity grad - constitutive_state(g,i,e)%p = constitutive_subState0(g,i,e)%p ! ...microstructure - crystallite_Tstar_v(:,g,i,e) = crystallite_subTstar0_v(:,g,i,e) ! ...2nd PK stress + crystallite_Fp(1:3,1:3,g,i,e) = crystallite_subFp0(1:3,1:3,g,i,e) ! ...plastic def grad + crystallite_invFp(1:3,1:3,g,i,e) = math_inv3x3(crystallite_Fp(1:3,1:3,g,i,e)) + crystallite_Lp(1:3,1:3,g,i,e) = crystallite_subLp0(1:3,1:3,g,i,e) ! ...plastic velocity grad + constitutive_state(g,i,e)%p = constitutive_subState0(g,i,e)%p ! ...microstructure + crystallite_Tstar_v(1:6,g,i,e) = crystallite_subTstar0_v(1:6,g,i,e) ! ...2nd PK stress ! cant restore dotState here, since not yet calculated in first cutback after initialization + !$OMP FLUSH(crystallite_subStep,crystallite_invFp) if (debugger .and. (e == debug_e .and. i == debug_i .and. g == debug_g)) then !$OMP CRITICAL (write2out) write(6,'(a78,f10.8)') 'cutback step in crystallite_stressAndItsTangent with new crystallite_subStep: ',& @@ -598,13 +608,15 @@ do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2 ! --- prepare for integration --- + !$OMP FLUSH(crystallite_subStep) crystallite_todo(g,i,e) = crystallite_subStep(g,i,e) > subStepMinCryst ! still on track or already done (beyond repair) if (crystallite_todo(g,i,e)) then - crystallite_subF(:,:,g,i,e) = crystallite_subF0(:,:,g,i,e) + & - crystallite_subStep(g,i,e) * & - (crystallite_partionedF(:,:,g,i,e) - crystallite_partionedF0(:,:,g,i,e)) - crystallite_Fe(:,:,g,i,e) = math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), crystallite_invFp(1:3,1:3,g,i,e)) - crystallite_subdt(g,i,e) = crystallite_subStep(g,i,e) * crystallite_dt(g,i,e) + crystallite_subF(1:3,1:3,g,i,e) = crystallite_subF0(1:3,1:3,g,i,e) & + + crystallite_subStep(g,i,e) & + * (crystallite_partionedF(1:3,1:3,g,i,e) - crystallite_partionedF0(1:3,1:3,g,i,e)) + !$OMP FLUSH(crystallite_subF) + crystallite_Fe(1:3,1:3,g,i,e) = math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), crystallite_invFp(1:3,1:3,g,i,e)) + crystallite_subdt(g,i,e) = crystallite_subStep(g,i,e) * crystallite_dt(g,i,e) crystallite_converged(g,i,e) = .false. ! start out non-converged endif @@ -647,16 +659,17 @@ enddo Fe_guess = math_mul33x33(crystallite_partionedF(1:3,1:3,g,i,e), invFp) Tstar = math_Mandel6to33( math_mul66x6( 0.5_pReal*constitutive_homogenizedC(g,i,e), & math_Mandel33to6( math_mul33x33(transpose(Fe_guess),Fe_guess) - math_I3 ) ) ) - crystallite_P(:,:,g,i,e) = math_mul33x33(Fe_guess,math_mul33x33(Tstar,transpose(invFp))) + crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(Fe_guess,math_mul33x33(Tstar,transpose(invFp))) + !$OMP FLUSH(crystallite_P) endif if (debugger .and. (e == debug_e .and. i == debug_i .and. g == debug_g)) then !$OMP CRITICAL (write2out) write (6,*) '#############' write (6,*) 'central solution of cryst_StressAndTangent' write (6,*) '#############' - write (6,'(a8,3(x,i4),/,3(3(f12.4,x)/))') ' P of', g, i, e, math_transpose3x3(crystallite_P(:,:,g,i,e))/1e6 - write (6,'(a8,3(x,i4),/,3(3(f14.9,x)/))') ' Fp of', g, i, e, math_transpose3x3(crystallite_Fp(:,:,g,i,e)) - write (6,'(a8,3(x,i4),/,3(3(f14.9,x)/))') ' Lp of', g, i, e, math_transpose3x3(crystallite_Lp(:,:,g,i,e)) + write (6,'(a8,3(x,i4),/,3(3(f12.4,x)/))') ' P of', g, i, e, math_transpose3x3(crystallite_P(1:3,1:3,g,i,e)) / 1e6 + write (6,'(a8,3(x,i4),/,3(3(f14.9,x)/))') ' Fp of', g, i, e, math_transpose3x3(crystallite_Fp(1:3,1:3,g,i,e)) + write (6,'(a8,3(x,i4),/,3(3(f14.9,x)/))') ' Lp of', g, i, e, math_transpose3x3(crystallite_Lp(1:3,1:3,g,i,e)) !$OMP END CRITICAL (write2out) endif enddo @@ -735,9 +748,9 @@ if(updateJaco) then if (crystallite_requested(g,i,e) .and. crystallite_converged(g,i,e)) then ! converged state warrants stiffness update select case(perturbation) case(1) - dPdF_perturbation1(:,:,k,l,g,i,e) = (crystallite_P(:,:,g,i,e) - P_backup(:,:,g,i,e)) / myPert ! tangent dP_ij/dFg_kl + dPdF_perturbation1(1:3,1:3,k,l,g,i,e) = (crystallite_P(1:3,1:3,g,i,e) - P_backup(1:3,1:3,g,i,e)) / myPert ! tangent dP_ij/dFg_kl case(2) - dPdF_perturbation2(:,:,k,l,g,i,e) = (crystallite_P(:,:,g,i,e) - P_backup(:,:,g,i,e)) / myPert ! tangent dP_ij/dFg_kl + dPdF_perturbation2(1:3,1:3,k,l,g,i,e) = (crystallite_P(1:3,1:3,g,i,e) - P_backup(1:3,1:3,g,i,e)) / myPert ! tangent dP_ij/dFg_kl end select endif enddo; enddo; enddo @@ -782,19 +795,20 @@ if(updateJaco) then if (crystallite_requested(g,i,e) .and. crystallite_converged(g,i,e)) then ! central solution converged select case(pert_method) case(1) - crystallite_dPdF(:,:,:,:,g,i,e) = dPdF_perturbation1(:,:,:,:,g,i,e) + crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = dPdF_perturbation1(1:3,1:3,1:3,1:3,g,i,e) case(2) - crystallite_dPdF(:,:,:,:,g,i,e) = dPdF_perturbation2(:,:,:,:,g,i,e) + crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = dPdF_perturbation2(1:3,1:3,1:3,1:3,g,i,e) case(3) - crystallite_dPdF(:,:,:,:,g,i,e) = 0.5_pReal* (dPdF_perturbation1(:,:,:,:,g,i,e) + dPdF_perturbation2(:,:,:,:,g,i,e)) + crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = 0.5_pReal* ( dPdF_perturbation1(1:3,1:3,1:3,1:3,g,i,e) & + + dPdF_perturbation2(1:3,1:3,1:3,1:3,g,i,e)) end select elseif (crystallite_requested(g,i,e) .and. .not. crystallite_converged(g,i,e)) then ! central solution did not converge - crystallite_dPdF(:,:,:,:,g,i,e) = crystallite_fallbackdPdF(:,:,:,:,g,i,e) ! use (elastic) fallback + crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = crystallite_fallbackdPdF(1:3,1:3,1:3,1:3,g,i,e) ! use (elastic) fallback endif enddo; enddo; enddo !OMP END PARALLEL DO -endif ! jacobian calculation +endif ! jacobian calculation endsubroutine @@ -811,7 +825,6 @@ use prec, only: pInt, & pReal use numerics, only: numerics_integrationMode use debug, only: debugger, & - selectiveDebugger, & verboseDebugger, & debug_e, & debug_i, & @@ -861,14 +874,14 @@ logical singleRun ! flag i if (present(ee) .and. present(ii) .and. present(gg)) then eIter = ee - iIter(:,ee) = ii - gIter(:,ee) = gg + iIter(1:2,ee) = ii + gIter(1:2,ee) = gg singleRun = .true. else eIter = FEsolving_execElem(1:2) do e = eIter(1),eIter(2) - iIter(:,e) = FEsolving_execIP(1:2,e) - gIter(:,e) = (/1,homogenization_Ngrains(mesh_element(3,e))/) + iIter(1:2,e) = FEsolving_execIP(1:2,e) + gIter(1:2,e) = (/1,homogenization_Ngrains(mesh_element(3,e))/) enddo singleRun = .false. endif @@ -879,9 +892,9 @@ endif !$OMP PARALLEL PRIVATE(mySizeDotState) !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), & + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states endif constitutive_dotState(g,i,e)%p = 0.0_pReal ! reset dotState to zero @@ -893,15 +906,21 @@ endif RK4dotTemperature = 0.0_pReal ! initialize Runge-Kutta dotTemperature !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains constitutive_RK4dotState(g,i,e)%p = 0.0_pReal ! initialize Runge-Kutta dotState if (crystallite_todo(g,i,e)) then - call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_Fe, crystallite_Fp, & + call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_orientation, g,i,e) - crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), & + crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Temperature(g,i,e),g,i,e) - if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState - .or. crystallite_dotTemperature(g,i,e)/=crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature + endif + enddo; enddo; enddo +!$OMP ENDDO +!$OMP DO + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + if ( any(constitutive_dotState(g,i,e)%p /= constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState + .or. crystallite_dotTemperature(g,i,e) /= crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped @@ -922,14 +941,23 @@ do n = 1,4 ! --- state update --- !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then mySizeDotState = constitutive_sizeDotState(g,i,e) - constitutive_RK4dotState(g,i,e)%p = constitutive_RK4dotState(g,i,e)%p + weight(n)*constitutive_dotState(g,i,e)%p - RK4dotTemperature(g,i,e) = RK4dotTemperature(g,i,e) + weight(n)*crystallite_dotTemperature(g,i,e) - if (n == 4) then - constitutive_dotState(g,i,e)%p = constitutive_RK4dotState(g,i,e)%p / 6.0_pReal ! use weighted RKdotState for final integration + if (n < 4) then + constitutive_RK4dotState(g,i,e)%p = constitutive_RK4dotState(g,i,e)%p + weight(n)*constitutive_dotState(g,i,e)%p + RK4dotTemperature(g,i,e) = RK4dotTemperature(g,i,e) + weight(n)*crystallite_dotTemperature(g,i,e) + elseif (n == 4) then + constitutive_dotState(g,i,e)%p = (constitutive_RK4dotState(g,i,e)%p + weight(n)*constitutive_dotState(g,i,e)%p) /6.0_pReal ! use weighted RKdotState for final integration + crystallite_dotTemperature(g,i,e) = (RK4dotTemperature(g,i,e) + weight(n)*crystallite_dotTemperature(g,i,e)) / 6.0_pReal endif + endif + enddo; enddo; enddo + !$OMP ENDDO + !$OMP DO + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + mySizeDotState = constitutive_sizeDotState(g,i,e) constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_subState0(g,i,e)%p(1:mySizeDotState) & + constitutive_dotState(g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e) * timeStepFraction(n) crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) & @@ -942,9 +970,9 @@ do n = 1,4 ! --- update dependent states --- !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), & + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states endif enddo; enddo; enddo @@ -954,11 +982,11 @@ do n = 1,4 ! --- stress integration --- !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then if (crystallite_integrateStress(g,i,e,timeStepFraction(n))) then ! fraction of original times step if (n == 4) then ! final integration step - if (verboseDebugger .and. selectiveDebugger .and. e == debug_e .and. i == debug_i .and. g == debug_g) then + if (verboseDebugger .and. e == debug_e .and. i == debug_i .and. g == debug_g) then mySizeDotState = constitutive_sizeDotState(g,i,e) !$OMP CRITICAL (write2out) write(6,*) '::: updateState',g,i,e @@ -994,15 +1022,21 @@ do n = 1,4 if (n < 4) then !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_Fe, crystallite_Fp, & + call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & crystallite_Temperature(g,i,e), timeStepFraction(n)*crystallite_subdt(g,i,e), & ! fraction of original timestep crystallite_orientation, g,i,e) - crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), & + crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Temperature(g,i,e),g,i,e) - if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState - .or. crystallite_dotTemperature(g,i,e)/=crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature + endif + enddo; enddo; enddo + !$OMP ENDDO + !$OMP DO + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + if ( any(constitutive_dotState(g,i,e)%p /= constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState + .or. crystallite_dotTemperature(g,i,e) /= crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped @@ -1045,7 +1079,6 @@ subroutine crystallite_integrateStateRKCK45(gg,ii,ee) use prec, only: pInt, & pReal use debug, only: debugger, & - selectiveDebugger, & verboseDebugger, & debug_e, & debug_i, & @@ -1155,14 +1188,14 @@ c(5) = 0.875_pReal if (present(ee) .and. present(ii) .and. present(gg)) then eIter = ee - iIter(:,ee) = ii - gIter(:,ee) = gg + iIter(1:2,ee) = ii + gIter(1:2,ee) = gg singleRun = .true. else eIter = FEsolving_execElem(1:2) do e = eIter(1),eIter(2) - iIter(:,e) = FEsolving_execIP(1:2,e) - gIter(:,e) = (/1,homogenization_Ngrains(mesh_element(3,e))/) + iIter(1:2,e) = FEsolving_execIP(1:2,e) + gIter(1:2,e) = (/1,homogenization_Ngrains(mesh_element(3,e))/) enddo singleRun = .false. endif @@ -1173,9 +1206,9 @@ endif !$OMP PARALLEL PRIVATE(mySizeDotState) !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), & + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states endif constitutive_dotState(g,i,e)%p = 0.0_pReal ! reset dotState to zero @@ -1192,14 +1225,20 @@ if (verboseDebugger) then !$OMP END SINGLE endif !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_Fe, crystallite_Fp, & + call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_orientation, g,i,e) - crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), & + crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Temperature(g,i,e),g,i,e) - if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState - .or. crystallite_dotTemperature(g,i,e)/=crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature + endif + enddo; enddo; enddo +!$OMP ENDDO +!$OMP DO + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + if ( any(constitutive_dotState(g,i,e)%p /= constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState + .or. crystallite_dotTemperature(g,i,e) /= crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped @@ -1220,17 +1259,60 @@ do n = 1,5 ! --- state update --- !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then mySizeDotState = constitutive_sizeDotState(g,i,e) constitutive_RKCK45dotState(n,g,i,e)%p = constitutive_dotState(g,i,e)%p ! store Runge-Kutta dotState RKCK45dotTemperature(n,g,i,e) = crystallite_dotTemperature(g,i,e) ! store Runge-Kutta dotTemperature - constitutive_dotState(g,i,e)%p = 0.0_pReal ! reset dotState - crystallite_dotTemperature(g,i,e) = 0.0_pReal ! reset dotTemperature - do j = 1,n ! rates as combination of Runge-Kutta rates - constitutive_dotState(g,i,e)%p = constitutive_dotState(g,i,e)%p + a(j,n) * constitutive_RKCK45dotState(j,g,i,e)%p - crystallite_dotTemperature(g,i,e) = crystallite_dotTemperature(g,i,e) + a(j,n) * RKCK45dotTemperature(j,g,i,e) - enddo + endif + enddo; enddo; enddo + !$OMP ENDDO + !$OMP DO + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + if (n == 1) then ! NEED TO DO THE ADDITION IN THIS LENGTHY WAY BECAUSE OF PARALLELIZATION (CAN'T USE A REDUCTION CLAUSE ON A POINTER OR USER DEFINED TYPE) + constitutive_dotState(g,i,e)%p = a(1,1) * constitutive_RKCK45dotState(1,g,i,e)%p + crystallite_dotTemperature(g,i,e) = a(1,1) * RKCK45dotTemperature(1,g,i,e) + elseif (n == 2) then + constitutive_dotState(g,i,e)%p = a(1,2) * constitutive_RKCK45dotState(1,g,i,e)%p & + + a(2,2) * constitutive_RKCK45dotState(2,g,i,e)%p + crystallite_dotTemperature(g,i,e) = a(1,2) * RKCK45dotTemperature(1,g,i,e) & + + a(2,2) * RKCK45dotTemperature(2,g,i,e) + elseif (n == 3) then + constitutive_dotState(g,i,e)%p = a(1,3) * constitutive_RKCK45dotState(1,g,i,e)%p & + + a(2,3) * constitutive_RKCK45dotState(2,g,i,e)%p & + + a(3,3) * constitutive_RKCK45dotState(3,g,i,e)%p + crystallite_dotTemperature(g,i,e) = a(1,3) * RKCK45dotTemperature(1,g,i,e) & + + a(2,3) * RKCK45dotTemperature(2,g,i,e) & + + a(3,3) * RKCK45dotTemperature(3,g,i,e) + elseif (n == 4) then + constitutive_dotState(g,i,e)%p = a(1,4) * constitutive_RKCK45dotState(1,g,i,e)%p & + + a(2,4) * constitutive_RKCK45dotState(2,g,i,e)%p & + + a(3,4) * constitutive_RKCK45dotState(3,g,i,e)%p & + + a(4,4) * constitutive_RKCK45dotState(4,g,i,e)%p + crystallite_dotTemperature(g,i,e) = a(1,4) * RKCK45dotTemperature(1,g,i,e) & + + a(2,4) * RKCK45dotTemperature(2,g,i,e) & + + a(3,4) * RKCK45dotTemperature(3,g,i,e) & + + a(4,4) * RKCK45dotTemperature(4,g,i,e) + elseif (n == 5) then + constitutive_dotState(g,i,e)%p = a(1,5) * constitutive_RKCK45dotState(1,g,i,e)%p & + + a(2,5) * constitutive_RKCK45dotState(2,g,i,e)%p & + + a(3,5) * constitutive_RKCK45dotState(3,g,i,e)%p & + + a(4,5) * constitutive_RKCK45dotState(4,g,i,e)%p & + + a(5,5) * constitutive_RKCK45dotState(5,g,i,e)%p + crystallite_dotTemperature(g,i,e) = a(1,5) * RKCK45dotTemperature(1,g,i,e) & + + a(2,5) * RKCK45dotTemperature(2,g,i,e) & + + a(3,5) * RKCK45dotTemperature(3,g,i,e) & + + a(4,5) * RKCK45dotTemperature(4,g,i,e) & + + a(5,5) * RKCK45dotTemperature(5,g,i,e) + endif + endif + enddo; enddo; enddo + !$OMP ENDDO + !$OMP DO + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + mySizeDotState = constitutive_sizeDotState(g,i,e) constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_subState0(g,i,e)%p(1:mySizeDotState) & + constitutive_dotState(g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e) crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) & @@ -1243,9 +1325,9 @@ do n = 1,5 ! --- update dependent states --- !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), & + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states endif constitutive_dotState(g,i,e)%p = 0.0_pReal ! reset dotState to zero @@ -1256,7 +1338,7 @@ do n = 1,5 ! --- stress integration --- !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then if (.not. crystallite_integrateStress(g,i,e,c(n))) then ! fraction of original time step if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... @@ -1281,13 +1363,19 @@ do n = 1,5 !$OMP END SINGLE endif !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_Fe, crystallite_Fp, & + call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & crystallite_Temperature(g,i,e), c(n)*crystallite_subdt(g,i,e), & ! fraction of original timestep crystallite_orientation, g,i,e) - crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), & + crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Temperature(g,i,e),g,i,e) + endif + enddo; enddo; enddo + !$OMP ENDDO + !$OMP DO + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState .or. crystallite_dotTemperature(g,i,e)/=crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... @@ -1307,55 +1395,92 @@ enddo ! --- STATE UPDATE WITH ERROR ESTIMATE FOR STATE AND TEMPERATURE --- -stateResiduum = 0.0_pReal -temperatureResiduum = 0.0_pReal relStateResiduum = 0.0_pReal relTemperatureResiduum = 0.0_pReal !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then mySizeDotState = constitutive_sizeDotState(g,i,e) constitutive_RKCK45dotState(6,g,i,e)%p = constitutive_dotState(g,i,e)%p ! store Runge-Kutta dotState RKCK45dotTemperature(6,g,i,e) = crystallite_dotTemperature(g,i,e) ! store Runge-Kutta dotTemperature + endif + enddo; enddo; enddo +!$OMP ENDDO - +!$OMP DO + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + mySizeDotState = constitutive_sizeDotState(g,i,e) + ! --- absolute residuum in state and temperature --- + ! NEED TO DO THE ADDITION IN THIS LENGTHY WAY BECAUSE OF PARALLELIZATION + ! CAN'T USE A REDUCTION CLAUSE ON A POINTER OR USER DEFINED TYPE - do j = 1,6 - stateResiduum(1:mySizeDotState,g,i,e) = stateResiduum(1:mySizeDotState,g,i,e) & - + db(j) * constitutive_RKCK45dotState(j,g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e) - temperatureResiduum(g,i,e) = temperatureResiduum(g,i,e) + db(j) * RKCK45dotTemperature(j,g,i,e) * crystallite_subdt(g,i,e) - enddo - - + stateResiduum(1:mySizeDotState,g,i,e) = ( db(1) * constitutive_RKCK45dotState(1,g,i,e)%p(1:mySizeDotState) & + + db(2) * constitutive_RKCK45dotState(2,g,i,e)%p(1:mySizeDotState) & + + db(3) * constitutive_RKCK45dotState(3,g,i,e)%p(1:mySizeDotState) & + + db(4) * constitutive_RKCK45dotState(4,g,i,e)%p(1:mySizeDotState) & + + db(5) * constitutive_RKCK45dotState(5,g,i,e)%p(1:mySizeDotState) & + + db(6) * constitutive_RKCK45dotState(6,g,i,e)%p(1:mySizeDotState)) & + * crystallite_subdt(g,i,e) + temperatureResiduum(g,i,e) = ( db(1) * RKCK45dotTemperature(1,g,i,e) & + + db(2) * RKCK45dotTemperature(2,g,i,e) & + + db(3) * RKCK45dotTemperature(3,g,i,e) & + + db(4) * RKCK45dotTemperature(4,g,i,e) & + + db(5) * RKCK45dotTemperature(5,g,i,e) & + + db(6) * RKCK45dotTemperature(6,g,i,e)) & + * crystallite_subdt(g,i,e) + ! --- dot state and dot temperature --- - - constitutive_dotState(g,i,e)%p = 0.0_pReal - crystallite_dotTemperature(g,i,e) = 0.0_pReal - do j = 1,6 - constitutive_dotState(g,i,e)%p = constitutive_dotState(g,i,e)%p + b(j) * constitutive_RKCK45dotState(j,g,i,e)%p - crystallite_dotTemperature(g,i,e) = crystallite_dotTemperature(g,i,e) + b(j) * RKCK45dotTemperature(j,g,i,e) - enddo + constitutive_dotState(g,i,e)%p = b(1) * constitutive_RKCK45dotState(1,g,i,e)%p & + + b(2) * constitutive_RKCK45dotState(2,g,i,e)%p & + + b(3) * constitutive_RKCK45dotState(3,g,i,e)%p & + + b(4) * constitutive_RKCK45dotState(4,g,i,e)%p & + + b(5) * constitutive_RKCK45dotState(5,g,i,e)%p & + + b(6) * constitutive_RKCK45dotState(6,g,i,e)%p + crystallite_dotTemperature(g,i,e) = b(1) * RKCK45dotTemperature(1,g,i,e) & + + b(2) * RKCK45dotTemperature(2,g,i,e) & + + b(3) * RKCK45dotTemperature(3,g,i,e) & + + b(4) * RKCK45dotTemperature(4,g,i,e) & + + b(5) * RKCK45dotTemperature(5,g,i,e) & + + b(6) * RKCK45dotTemperature(6,g,i,e) + endif + enddo; enddo; enddo +!$OMP ENDDO - ! --- state and temperature update and relative residui --- +! --- state and temperature update --- +!$OMP DO + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + mySizeDotState = constitutive_sizeDotState(g,i,e) constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_subState0(g,i,e)%p(1:mySizeDotState) & + constitutive_dotState(g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e) crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) & + crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e) + endif + enddo; enddo; enddo +!$OMP ENDDO + +! --- relative residui and state convergence --- + +!$OMP DO + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + mySizeDotState = constitutive_sizeDotState(g,i,e) forall (s = 1:mySizeDotState, abs(constitutive_state(g,i,e)%p(s)) > 0.0_pReal) & relStateResiduum(s,g,i,e) = stateResiduum(s,g,i,e) / constitutive_state(g,i,e)%p(s) if (crystallite_Temperature(g,i,e) > 0) & relTemperatureResiduum(g,i,e) = temperatureResiduum(g,i,e) / crystallite_Temperature(g,i,e) + !$OMP FLUSH(relStateResiduum,relTemperatureResiduum) - ! --- state convergence --- crystallite_todo(g,i,e) = & ( all( abs(relStateResiduum(:,g,i,e)) < rTol_crystalliteState & .or. abs(stateResiduum(1:mySizeDotState,g,i,e)) < constitutive_aTolState(g,i,e)%p(1:mySizeDotState) ) & .and. abs(relTemperatureResiduum(g,i,e)) < rTol_crystalliteTemperature ) - + if (verboseDebugger .and. (e == debug_e .and. i == debug_i .and. g == debug_g)) then !$OMP CRITICAL (write2out) write(6,*) '::: updateState',g,i,e @@ -1386,9 +1511,9 @@ relTemperatureResiduum = 0.0_pReal ! --- UPDATE DEPENDENT STATES IF RESIDUUM BELOW TOLERANCE --- !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), & + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states endif enddo; enddo; enddo @@ -1398,7 +1523,7 @@ relTemperatureResiduum = 0.0_pReal ! --- FINAL STRESS INTEGRATION STEP IF RESIDUUM BELOW TOLERANCE --- !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then if (crystallite_integrateStress(g,i,e)) then crystallite_converged(g,i,e) = .true. ! ... converged per definitionem @@ -1421,7 +1546,11 @@ relTemperatureResiduum = 0.0_pReal ! --- nonlocal convergence check --- -if (verboseDebugger .and. numerics_integrationMode == 1) write(6,*) 'crystallite_converged',crystallite_converged +if (verboseDebugger .and. numerics_integrationMode == 1) then + !$OMP CRITICAL (write2out) + write(6,*) 'crystallite_converged',crystallite_converged + !$OMP END CRITICAL (write2out) +endif if (.not. singleRun) then ! if not requesting Integration of just a single IP if ( any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) then ! any non-local not yet converged (or broken)... crystallite_converged = crystallite_converged .and. crystallite_localConstitution ! ...restart all non-local as not converged @@ -1443,7 +1572,6 @@ subroutine crystallite_integrateStateAdaptiveEuler(gg,ii,ee) use prec, only: pInt, & pReal use debug, only: debugger, & - selectiveDebugger, & verboseDebugger, & debug_e, & debug_i, & @@ -1505,14 +1633,14 @@ logical singleRun ! flag i if (present(ee) .and. present(ii) .and. present(gg)) then eIter = ee - iIter(:,ee) = ii - gIter(:,ee) = gg + iIter(1:2,ee) = ii + gIter(1:2,ee) = gg singleRun = .true. else eIter = FEsolving_execElem(1:2) do e = eIter(1),eIter(2) - iIter(:,e) = FEsolving_execIP(1:2,e) - gIter(:,e) = (/1,homogenization_Ngrains(mesh_element(3,e))/) + iIter(1:2,e) = FEsolving_execIP(1:2,e) + gIter(1:2,e) = (/1,homogenization_Ngrains(mesh_element(3,e))/) enddo singleRun = .false. endif @@ -1523,9 +1651,9 @@ endif !$OMP PARALLEL PRIVATE(mySizeDotState) !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), & + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states endif constitutive_dotState(g,i,e)%p = 0.0_pReal ! reset dotState to zero @@ -1537,14 +1665,20 @@ endif stateResiduum = 0.0_pReal !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_Fe, crystallite_Fp, & + call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_orientation, g,i,e) - crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), & + crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Temperature(g,i,e),g,i,e) - if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState - .or. crystallite_dotTemperature(g,i,e)/=crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature + endif + enddo; enddo; enddo +!$OMP ENDDO +!$OMP DO + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + if ( any(constitutive_dotState(g,i,e)%p /= constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState + .or. crystallite_dotTemperature(g,i,e) /= crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped @@ -1561,7 +1695,7 @@ stateResiduum = 0.0_pReal ! --- STATE UPDATE (EULER INTEGRATION) --- !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then mySizeDotState = constitutive_sizeDotState(g,i,e) stateResiduum(1:mySizeDotState,g,i,e) = - 0.5_pReal * constitutive_dotState(g,i,e)%p * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state and temperature @@ -1578,9 +1712,9 @@ stateResiduum = 0.0_pReal ! --- UPDATE DEPENDENT STATES (EULER INTEGRATION) --- !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), & + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states endif constitutive_dotState(g,i,e)%p = 0.0_pReal ! reset dotState to zero @@ -1591,7 +1725,7 @@ stateResiduum = 0.0_pReal ! --- STRESS INTEGRATION (EULER INTEGRATION) --- !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then if (.not. crystallite_integrateStress(g,i,e)) then if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... @@ -1610,14 +1744,20 @@ stateResiduum = 0.0_pReal ! --- DOT STATE AND TEMPERATURE (HEUN METHOD) --- !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_Fe, crystallite_Fp, & + call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_orientation, g,i,e) - crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), & + crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Temperature(g,i,e),g,i,e) - if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState - .or. crystallite_dotTemperature(g,i,e)/=crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature + endif + enddo; enddo; enddo +!$OMP ENDDO +!$OMP DO + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then + if ( any(constitutive_dotState(g,i,e)%p /= constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState + .or. crystallite_dotTemperature(g,i,e) /= crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped @@ -1636,7 +1776,7 @@ stateResiduum = 0.0_pReal relStateResiduum = 0.0_pReal relTemperatureResiduum = 0.0_pReal !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then mySizeDotState = constitutive_sizeDotState(g,i,e) @@ -1647,7 +1787,7 @@ relTemperatureResiduum = 0.0_pReal + 0.5_pReal * constitutive_dotState(g,i,e)%p * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state and temperature temperatureResiduum(g,i,e) = temperatureResiduum(g,i,e) & + 0.5_pReal * crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e) - + !$OMP FLUSH(stateResiduum,temperatureResiduum) ! --- relative residui --- @@ -1655,6 +1795,7 @@ relTemperatureResiduum = 0.0_pReal relStateResiduum(s,g,i,e) = stateResiduum(s,g,i,e) / constitutive_state(g,i,e)%p(s) if (crystallite_Temperature(g,i,e) > 0) & relTemperatureResiduum(g,i,e) = temperatureResiduum(g,i,e) / crystallite_Temperature(g,i,e) + !$OMP FLUSH(relStateResiduum,relTemperatureResiduum) if (verboseDebugger .and. (e == debug_e .and. i == debug_i .and. g == debug_g)) then !$OMP CRITICAL (write2out) @@ -1697,7 +1838,11 @@ relTemperatureResiduum = 0.0_pReal ! --- NONLOCAL CONVERGENCE CHECK --- -if (verboseDebugger .and. numerics_integrationMode==1_pInt) write(6,*) 'crystallite_converged',crystallite_converged +if (verboseDebugger .and. numerics_integrationMode==1_pInt) then + !$OMP CRITICAL (write2out) + write(6,*) 'crystallite_converged',crystallite_converged + !$OMP END CRITICAL (write2out) +endif if (.not. singleRun) then ! if not requesting Integration of just a single IP if ( any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) then ! any non-local not yet converged (or broken)... crystallite_converged = crystallite_converged .and. crystallite_localConstitution ! ...restart all non-local as not converged @@ -1761,14 +1906,14 @@ logical singleRun ! flag i if (present(ee) .and. present(ii) .and. present(gg)) then eIter = ee - iIter(:,ee) = ii - gIter(:,ee) = gg + iIter(1:2,ee) = ii + gIter(1:2,ee) = gg singleRun = .true. else eIter = FEsolving_execElem(1:2) do e = eIter(1),eIter(2) - iIter(:,e) = FEsolving_execIP(1:2,e) - gIter(:,e) = (/1,homogenization_Ngrains(mesh_element(3,e))/) + iIter(1:2,e) = FEsolving_execIP(1:2,e) + gIter(1:2,e) = (/1,homogenization_Ngrains(mesh_element(3,e))/) enddo singleRun = .false. endif @@ -1779,9 +1924,9 @@ endif !$OMP PARALLEL PRIVATE(mySizeDotState) !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), & + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states endif constitutive_dotState(g,i,e)%p = 0.0_pReal ! reset dotState to zero @@ -1792,12 +1937,18 @@ endif ! --- DOT STATE AND TEMPERATURE --- !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_Fe, crystallite_Fp, & + call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_orientation, g,i,e) - crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), & + crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Temperature(g,i,e),g,i,e) + endif + enddo; enddo; enddo +!$OMP ENDDO +!$OMP DO + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + if (crystallite_todo(g,i,e)) then if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState .or. crystallite_dotTemperature(g,i,e)/=crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... @@ -1816,35 +1967,34 @@ endif ! --- UPDATE STATE AND TEMPERATURE --- !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then mySizeDotState = constitutive_sizeDotState(g,i,e) constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_subState0(g,i,e)%p(1:mySizeDotState) & + constitutive_dotState(g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e) crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) & + crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e) - - if (verboseDebugger .and. (e == debug_e .and. i == debug_i .and. g == debug_g)) then - !$OMP CRITICAL (write2out) - write(6,*) '::: updateState',g,i,e - write(6,*) - write(6,'(a,/,12(e12.5,x))') 'updateState: dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState) - write(6,*) - write(6,'(a,/,12(e12.5,x))') 'updateState: new state', constitutive_state(g,i,e)%p(1:mySizeDotState) - write(6,*) - !$OMP END CRITICAL (write2out) - endif endif enddo; enddo; enddo !$OMP ENDDO +if (verboseDebugger .and. selectiveDebugger) then + !$OMP CRITICAL (write2out) + write(6,*) '::: updateState', debug_g, debug_i, debug_e + write(6,*) + write(6,'(a,/,12(e12.5,x))') 'updateState: dotState', constitutive_dotState(debug_g,debug_i,debug_e)%p(1:mySizeDotState) + write(6,*) + write(6,'(a,/,12(e12.5,x))') 'updateState: new state', constitutive_state(debug_g,debug_i,debug_e)%p(1:mySizeDotState) + write(6,*) + !$OMP END CRITICAL (write2out) +endif ! --- UPDATE DEPENDENT STATES --- !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), & + call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states endif enddo; enddo; enddo @@ -1854,7 +2004,7 @@ endif ! --- STRESS INTEGRATION --- !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then if (crystallite_integrateStress(g,i,e)) then crystallite_converged(g,i,e) = .true. @@ -1899,7 +2049,6 @@ subroutine crystallite_integrateStateFPI(gg,ii,ee) use prec, only: pInt, & pReal use debug, only: debugger, & - selectiveDebugger, & verboseDebugger, & debug_e, & debug_i, & @@ -1940,19 +2089,21 @@ integer(pInt), dimension(2,mesh_NcpElems) :: iIter, & ! bounds gIter ! bounds for grain iteration real(pReal) dot_prod12, & dot_prod22 -logical singleRun ! flag indicating computation for single (g,i,e) triple +logical singleRun, & ! flag indicating computation for single (g,i,e) triple + stateConverged, & ! flag indicating convergence of state integration + temperatureConverged ! flag indicating convergence of temperature integration if (present(ee) .and. present(ii) .and. present(gg)) then eIter = ee - iIter(:,ee) = ii - gIter(:,ee) = gg + iIter(1:2,ee) = ii + gIter(1:2,ee) = gg singleRun = .true. else eIter = FEsolving_execElem(1:2) do e = eIter(1),eIter(2) - iIter(:,e) = FEsolving_execIP(1:2,e) - gIter(:,e) = (/1,homogenization_Ngrains(mesh_element(3,e))/) + iIter(1:2,e) = FEsolving_execIP(1:2,e) + gIter(1:2,e) = (/1,homogenization_Ngrains(mesh_element(3,e))/) enddo singleRun = .false. endif @@ -1962,10 +2113,10 @@ endif ! --- RESET DEPENDENT STATES AND DOTSTATE --- -!$OMP PARALLEL +!$OMP PARALLEL PRIVATE(stateConverged,temperatureConverged) !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states @@ -1980,9 +2131,9 @@ endif ! --- DOT STATES --- !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_Fe, crystallite_Fp, & + call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_orientation, g,i,e) endif enddo; enddo; enddo @@ -1995,10 +2146,12 @@ endif crystallite_statedamper = 1.0_pReal !$OMP END SINGLE !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - crystallite_stateConverged(g,i,e) = crystallite_updateState(g,i,e) ! update state - crystallite_temperatureConverged(g,i,e) = crystallite_updateTemperature(g,i,e) ! update temperature + call crystallite_updateState(crystallite_todo(g,i,e), stateConverged, g,i,e) ! update state + !$OMP FLUSH(crystallite_todo) + call crystallite_updateTemperature(crystallite_todo(g,i,e), temperatureConverged, g,i,e) ! update temperature + !$OMP FLUSH(crystallite_todo) if ( .not. crystallite_localConstitution(g,i,e) .and. .not. crystallite_todo(g,i,e) ) then ! if updateState or updateTemperature signals broken non-local... !$OMP CRITICAL (checkTodo) crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped @@ -2012,7 +2165,7 @@ endif ! --- UPDATE DEPENDENT STATES --- !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states @@ -2039,13 +2192,16 @@ do while (any(crystallite_todo) .and. NiterationState < nState ) ! --- STRESS INTEGRATION --- !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e) - if ( .not. crystallite_localConstitution(g,i,e) .and. .not. crystallite_todo(g,i,e)) then ! if broken non-local... - !$OMP CRITICAL (checkTodo) - crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped - !$OMP END CRITICAL (checkTodo) + if (.not. crystallite_integrateStress(g,i,e)) then ! if broken ... + if (.not. crystallite_localConstitution(g,i,e)) then ! ... and non-local... + !$OMP CRITICAL (checkTodo) + crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ... then all non-locals skipped + !$OMP END CRITICAL (checkTodo) + else ! ... and local... + crystallite_todo(g,i,e) = .false. ! ... then skip only me + endif endif endif enddo; enddo; enddo @@ -2063,9 +2219,9 @@ do while (any(crystallite_todo) .and. NiterationState < nState ) ! --- DOT STATES --- !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_Fe, crystallite_Fp, & + call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_orientation, g,i,e) endif enddo; enddo; enddo @@ -2078,7 +2234,7 @@ do while (any(crystallite_todo) .and. NiterationState < nState ) crystallite_statedamper = 1.0_pReal !$OMP END SINGLE !$OMP DO PRIVATE(dot_prod12,dot_prod22) - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then ! --- state damper --- @@ -2089,19 +2245,23 @@ do while (any(crystallite_todo) .and. NiterationState < nState ) constitutive_previousDotState(g,i,e)%p - constitutive_previousDotState2(g,i,e)%p ) if ( dot_prod22 > 0.0_pReal & .and. ( dot_prod12 < 0.0_pReal & - .or. dot_product(constitutive_dotState(g,i,e)%p, constitutive_previousDotState(g,i,e)%p) < 0.0_pReal) ) & + .or. dot_product(constitutive_dotState(g,i,e)%p, constitutive_previousDotState(g,i,e)%p) < 0.0_pReal) ) then crystallite_statedamper(g,i,e) = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) - + !$OMP FLUSH(crystallite_statedamper) + endif + ! --- updates --- - crystallite_stateConverged(g,i,e) = crystallite_updateState(g,i,e) ! update state - crystallite_temperatureConverged(g,i,e) = crystallite_updateTemperature(g,i,e) ! update temperature - crystallite_converged(g,i,e) = crystallite_stateConverged(g,i,e) .and. crystallite_temperatureConverged(g,i,e) + call crystallite_updateState(crystallite_todo(g,i,e), stateConverged, g,i,e) ! update state + !$OMP FLUSH(crystallite_todo) + call crystallite_updateTemperature(crystallite_todo(g,i,e), temperatureConverged, g,i,e) ! update temperature + !$OMP FLUSH(crystallite_todo) + crystallite_converged(g,i,e) = stateConverged .and. temperatureConverged if ( .not. crystallite_localConstitution(g,i,e) .and. .not. crystallite_todo(g,i,e)) then ! if updateState or updateTemperature signals broken non-local... !$OMP CRITICAL (checkTodo) crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped !$OMP END CRITICAL (checkTodo) - elseif (crystallite_converged(g,i,e)) then + elseif (stateConverged .and. temperatureConverged) then ! check (private) logicals "stateConverged" and "temperatureConverged" instead of (shared) "crystallite_converged", so no need to flush the "crystallite_converged" array !$OMP CRITICAL (distributionState) debug_StateLoopDistribution(NiterationState,numerics_integrationMode) = & debug_StateLoopDistribution(NiterationState,numerics_integrationMode) + 1 @@ -2115,7 +2275,7 @@ do while (any(crystallite_todo) .and. NiterationState < nState ) ! --- UPDATE DEPENDENT STATES --- !$OMP DO - do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states @@ -2140,12 +2300,16 @@ do while (any(crystallite_todo) .and. NiterationState < nState ) ! --- CONVERGENCE CHECK --- if (.not. singleRun) then ! if not requesting Integration of just a single IP - if (any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) then ! any non-local not yet converged (or broken)... - crystallite_converged = crystallite_converged .and. crystallite_localConstitution ! ...restart all non-local as not converged - endif + !$OMP SINGLE + if (any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) then ! any non-local not yet converged (or broken)... + crystallite_converged = crystallite_converged .and. crystallite_localConstitution ! ...restart all non-local as not converged + endif + !$OMP END SINGLE endif - crystallite_todo = crystallite_todo .and. .not. crystallite_converged ! skip all converged + !$OMP SINGLE + crystallite_todo = crystallite_todo .and. .not. crystallite_converged ! skip all converged + !$OMP END SINGLE if (verboseDebugger .and. numerics_integrationMode == 1) then !$OMP SINGLE @@ -2169,7 +2333,7 @@ endsubroutine ! update the internal state of the constitutive law ! and tell whether state has converged !******************************************************************** -function crystallite_updateState(g,i,e) +subroutine crystallite_updateState(done, converged, g, i, e) !*** variables and functions from other modules ***! use prec, only: pReal, & @@ -2191,30 +2355,39 @@ use debug, only: debugger, & verboseDebugger !*** input variables ***! -integer(pInt), intent(in):: e, & ! element index - i, & ! integration point index - g ! grain index +integer(pInt), intent(in):: e, & ! element index + i, & ! integration point index + g ! grain index !*** output variables ***! -logical crystallite_updateState ! flag indicating if integration suceeded +logical, intent(out) :: converged, & ! flag indicating if state converged + done ! flag indicating if state was updated !*** local variables ***! -real(pReal), dimension(constitutive_sizeDotState(g,i,e)) :: residuum ! residuum from evolution of microstructure +real(pReal), dimension(constitutive_sizeDotState(g,i,e)) :: & + residuum, & ! residuum from evolution of microstructure + dotState, & ! local copy of dotState + state ! local copy of relevant part of state integer(pInt) mySize +!* start out as not done and not converged +!* and make local copies of state and dotState (needed for parallelization, since no chance of flushing a pointer) + +done = .false. +converged = .false. mySize = constitutive_sizeDotState(g,i,e) +dotState(1:mySize) = constitutive_dotState(g,i,e)%p(1:mySize) +state(1:mySize) = constitutive_state(g,i,e)%p(1:mySize) -! correct my dotState -constitutive_dotState(g,i,e)%p(1:mySize) = constitutive_dotState(g,i,e)%p(1:mySize) * crystallite_statedamper(g,i,e) & - + constitutive_previousDotState(g,i,e)%p(1:mySize) * (1.0_pReal-crystallite_statedamper(g,i,e)) +!* correct dotState, calculate residuum and check if it is valid + +dotState(1:mySize) = constitutive_dotState(g,i,e)%p(1:mySize) * crystallite_statedamper(g,i,e) & + + constitutive_previousDotState(g,i,e)%p(1:mySize) * (1.0_pReal - crystallite_statedamper(g,i,e)) residuum = constitutive_state(g,i,e)%p(1:mySize) - constitutive_subState0(g,i,e)%p(1:mySize) & - - constitutive_dotState(g,i,e)%p(1:mySize) * crystallite_subdt(g,i,e) - -if (any(residuum/=residuum)) then ! if NaN occured then return without changing the state... - crystallite_updateState = .false. ! ...indicate state update failed - crystallite_todo(g,i,e) = .false. ! ...no need to calculate any further + - dotState(1:mySize) * crystallite_subdt(g,i,e) +if (any(residuum /= residuum)) then ! if NaN occured then return without changing the state if (verboseDebugger) then !$OMP CRITICAL (write2out) write(6,*) '::: updateState encountered NaN',g,i,e @@ -2223,15 +2396,19 @@ if (any(residuum/=residuum)) then ! if NaN oc return endif -constitutive_state(g,i,e)%p(1:mySize) = constitutive_state(g,i,e)%p(1:mySize) - residuum -! setting flag to true if residuum is below relative/absolute tolerance, otherwise set it to false -crystallite_updateState = all( abs(residuum) < constitutive_aTolState(g,i,e)%p(1:mySize) & - .or. abs(residuum) < rTol_crystalliteState*abs(constitutive_state(g,i,e)%p(1:mySize)) ) +!* update state and set convergence flag to true if residuum is below relative/absolute tolerance, otherwise set it to false -if (verboseDebugger .and. (e == debug_e .and. i == debug_i .and. g == debug_g) .and. numerics_integrationMode == 1_pInt) then +state(1:mySize) = state(1:mySize) - residuum +done = .true. +converged = all( abs(residuum) < constitutive_aTolState(g,i,e)%p(1:mySize) & + .or. abs(residuum) < rTol_crystalliteState * abs(state(1:mySize)) ) + +if ( verboseDebugger & + .and. (e == debug_e .and. i == debug_i .and. g == debug_g) & + .and. numerics_integrationMode == 1_pInt) then !$OMP CRITICAL (write2out) - if (crystallite_updateState) then + if (converged) then write(6,*) '::: updateState converged',g,i,e else write(6,*) '::: updateState did not converge',g,i,e @@ -2239,17 +2416,20 @@ if (verboseDebugger .and. (e == debug_e .and. i == debug_i .and. g == debug_g) . write(6,*) write(6,'(a,f6.1)') 'updateState: crystallite_statedamper',crystallite_statedamper(g,i,e) write(6,*) - write(6,'(a,/,12(e12.5,x))') 'updateState: dotState',constitutive_dotState(g,i,e)%p(1:mySize) + write(6,'(a,/,12(e12.5,x))') 'updateState: dotState',dotState(1:mySize) write(6,*) - write(6,'(a,/,12(e12.5,x))') 'updateState: new state',constitutive_state(g,i,e)%p(1:mySize) + write(6,'(a,/,12(e12.5,x))') 'updateState: new state',state(1:mySize) write(6,*) -! write(6,'(a,/,12(f12.1,x))') 'updateState: relative residuum tolerance', abs(residuum / rTol_crystalliteState & -! / constitutive_state(g,i,e)%p(1:mySize)) -! write(6,*) !$OMP END CRITICAL (write2out) endif -endfunction + +!* sync global state and dotState with local copies + +constitutive_dotState(g,i,e)%p(1:mySize) = dotState(1:mySize) +constitutive_state(g,i,e)%p(1:mySize) = state(1:mySize) + +endsubroutine @@ -2257,56 +2437,58 @@ endfunction ! update the temperature of the grain ! and tell whether it has converged !******************************************************************** - function crystallite_updateTemperature(& - g,& ! grain number - i,& ! integration point number - e & ! element number - ) + subroutine crystallite_updateTemperature(done, converged, g, i, e) - !*** variables and functions from other modules ***! - use prec, only: pReal, & - pInt, & - pLongInt - use numerics, only: rTol_crystalliteTemperature - use constitutive, only: constitutive_dotTemperature - use debug, only: debugger +!*** variables and functions from other modules ***! +use prec, only: pReal, & + pInt, & + pLongInt +use numerics, only: rTol_crystalliteTemperature +use constitutive, only: constitutive_dotTemperature +use debug, only: debugger - !*** input variables ***! - integer(pInt), intent(in):: e, & ! element index - i, & ! integration point index - g ! grain index +!*** input variables ***! +integer(pInt), intent(in):: e, & ! element index + i, & ! integration point index + g ! grain index - !*** output variables ***! - logical crystallite_updateTemperature ! flag indicating if integration suceeded +!*** output variables ***! +logical, intent(out) :: converged, & ! flag indicating if temperature converged + done ! flag indicating if temperature was updated - !*** local variables ***! - real(pReal) residuum ! residuum from evolution of temperature - - ! calculate the residuum - residuum = crystallite_Temperature(g,i,e) - crystallite_subTemperature0(g,i,e) - & - crystallite_subdt(g,i,e) * & - constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e),crystallite_Temperature(g,i,e),g,i,e) - - ! if NaN occured then return without changing the state - if (residuum/=residuum) then - crystallite_updateTemperature = .false. ! indicate update failed - crystallite_todo(g,i,e) = .false. ! ...no need to calculate any further - !$OMP CRITICAL (write2out) - write(6,*) '::: updateTemperature encountered NaN',g,i,e - !$OMP END CRITICAL (write2out) - return - endif - - ! update the microstructure - crystallite_Temperature(g,i,e) = crystallite_Temperature(g,i,e) - residuum - - ! setting flag to true if residuum is below relative tolerance (or zero Kelvin), otherwise set it to false - crystallite_updateTemperature = crystallite_Temperature(g,i,e) == 0.0_pReal .or. & - abs(residuum) < rTol_crystalliteTemperature*crystallite_Temperature(g,i,e) +!*** local variables ***! +real(pReal) residuum ! residuum from evolution of temperature + +!* start out as not done and not converged + +done = .false. +converged = .false. + + +!* correct dotState, calculate residuum and check if it is valid +!* (if NaN occured then return without changing the temperature) + +residuum = crystallite_Temperature(g,i,e) - crystallite_subTemperature0(g,i,e) & + - constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e),crystallite_Temperature(g,i,e),g,i,e) & + * crystallite_subdt(g,i,e) +if (residuum /= residuum) then + !$OMP CRITICAL (write2out) + write(6,*) '::: updateTemperature encountered NaN',g,i,e + !$OMP END CRITICAL (write2out) return +endif + - endfunction +!* update temperature and set convergence flag to true if residuum is below relative tolerance (or zero Kelvin), otherwise set it to false + +crystallite_Temperature(g,i,e) = crystallite_Temperature(g,i,e) - residuum +done = .true. +!$OMP FLUSH(crystallite_Temperature) +converged = ( crystallite_Temperature(g,i,e) == 0.0_pReal & + .or. abs(residuum) < rTol_crystalliteTemperature * crystallite_Temperature(g,i,e)) + +endsubroutine @@ -2323,342 +2505,365 @@ function crystallite_integrateStress(& ) - !*** variables and functions from other modules ***! - use prec, only: pReal, & - pInt, & - pLongInt - use numerics, only: nStress, & - aTol_crystalliteStress, & - rTol_crystalliteStress, & - iJacoLpresiduum, & - relevantStrain, & - numerics_integrationMode - use debug, only: debugger, & - debug_g, & - debug_i, & - debug_e, & - verboseDebugger, & - debug_cumLpCalls, & - debug_cumLpTicks, & - debug_StressLoopDistribution, & - debug_LeapfrogBreakDistribution - use constitutive, only: constitutive_homogenizedC, & - constitutive_LpAndItsTangent - use math, only: math_mul33x33, & - math_mul66x6, & - math_mul99x99, & - math_transpose3x3, & - math_inv3x3, & - math_invert3x3, & - math_invert, & - math_det3x3, & - math_I3, & - math_identity2nd, & - math_Mandel66to3333, & - math_Mandel6to33, & - math_mandel33to6 +!*** variables and functions from other modules ***! +use prec, only: pReal, & + pInt, & + pLongInt +use numerics, only: nStress, & + aTol_crystalliteStress, & + rTol_crystalliteStress, & + iJacoLpresiduum, & + relevantStrain, & + numerics_integrationMode +use debug, only: debugger, & + debug_g, & + debug_i, & + debug_e, & + verboseDebugger, & + debug_cumLpCalls, & + debug_cumLpTicks, & + debug_StressLoopDistribution, & + debug_LeapfrogBreakDistribution +use constitutive, only: constitutive_homogenizedC, & + constitutive_LpAndItsTangent +use math, only: math_mul33x33, & + math_mul66x6, & + math_mul99x99, & + math_transpose3x3, & + math_inv3x3, & + math_invert3x3, & + math_invert, & + math_det3x3, & + math_I3, & + math_identity2nd, & + math_Mandel66to3333, & + math_Mandel6to33, & + math_mandel33to6 - implicit none +implicit none - !*** input variables ***! - integer(pInt), intent(in):: e, & ! element index - i, & ! integration point index - g ! grain index - real(pReal), optional, intent(in) :: fraction ! fraction of timestep +!*** input variables ***! +integer(pInt), intent(in):: e, & ! element index + i, & ! integration point index + g ! grain index +real(pReal), optional, intent(in) :: fraction ! fraction of timestep - !*** output variables ***! - logical crystallite_integrateStress ! flag indicating if integration suceeded +!*** output variables ***! +logical crystallite_integrateStress ! flag indicating if integration suceeded + +!*** local variables ***! +real(pReal), dimension(3,3):: Fg_new, & ! deformation gradient at end of timestep + Fp_current, & ! plastic deformation gradient at start of timestep + Fp_new, & ! plastic deformation gradient at end of timestep + Fe_new, & ! elastic deformation gradient at end of timestep + invFp_new, & ! inverse of Fp_new + invFp_current, & ! inverse of Fp_current + Lpguess, & ! current guess for plastic velocity gradient + Lpguess_old, & ! known last good guess for plastic velocity gradient + Lp_constitutive, & ! plastic velocity gradient resulting from constitutive law + residuum, & ! current residuum of plastic velocity gradient + residuum_old, & ! last residuum of plastic velocity gradient + A, & + B, & + BT, & + AB, & + BTA +real(pReal), dimension(6):: Tstar_v ! 2nd Piola-Kirchhoff Stress in Mandel-Notation +real(pReal), dimension(9,9):: dLpdT_constitutive, & ! partial derivative of plastic velocity gradient calculated by constitutive law + dTdLp, & ! partial derivative of 2nd Piola-Kirchhoff stress + dRdLp, & ! partial derivative of residuum (Jacobian for NEwton-Raphson scheme) + invdRdLp ! inverse of dRdLp +real(pReal), dimension(3,3,3,3):: C ! 4th rank elasticity tensor +real(pReal), dimension(6,6):: C_66 ! simplified 2nd rank elasticity tensor +real(pReal) p_hydro, & ! volumetric part of 2nd Piola-Kirchhoff Stress + det, & ! determinant + leapfrog, & ! acceleration factor for Newton-Raphson scheme + maxleap, & ! maximum acceleration factor + dt ! time increment +logical error ! flag indicating an error +integer(pInt) NiterationStress, & ! number of stress integrations + dummy, & + h, & + j, & + k, & + l, & + m, & + n, & + jacoCounter ! counter to check for Jacobian update +integer(pLongInt) tick, & + tock, & + tickrate, & + maxticks - !*** local variables ***! - real(pReal), dimension(3,3):: Fg_new, & ! deformation gradient at end of timestep - Fp_current, & ! plastic deformation gradient at start of timestep - Fp_new, & ! plastic deformation gradient at end of timestep - Fe_new, & ! elastic deformation gradient at end of timestep - invFp_new, & ! inverse of Fp_new - invFp_current, & ! inverse of Fp_current - Lpguess, & ! current guess for plastic velocity gradient - Lpguess_old, & ! known last good guess for plastic velocity gradient - Lp_constitutive, & ! plastic velocity gradient resulting from constitutive law - residuum, & ! current residuum of plastic velocity gradient - residuum_old, & ! last residuum of plastic velocity gradient - A, & - B, & - BT, & - AB, & - BTA - real(pReal), dimension(6):: Tstar_v ! 2nd Piola-Kirchhoff Stress in Mandel-Notation - real(pReal), dimension(9,9):: dLpdT_constitutive, & ! partial derivative of plastic velocity gradient calculated by constitutive law - dTdLp, & ! partial derivative of 2nd Piola-Kirchhoff stress - dRdLp, & ! partial derivative of residuum (Jacobian for NEwton-Raphson scheme) - invdRdLp ! inverse of dRdLp - real(pReal), dimension(3,3,3,3):: C ! 4th rank elasticity tensor - real(pReal), dimension(6,6):: C_66 ! simplified 2nd rank elasticity tensor - real(pReal) p_hydro, & ! volumetric part of 2nd Piola-Kirchhoff Stress - det, & ! determinant - leapfrog, & ! acceleration factor for Newton-Raphson scheme - maxleap, & ! maximum acceleration factor - dt ! time increment - logical error ! flag indicating an error - integer(pInt) NiterationStress, & ! number of stress integrations - dummy, & - h, & - j, & - k, & - l, & - m, & - n, & - jacoCounter ! counter to check for Jacobian update - integer(pLongInt) tick, & - tock, & - tickrate, & - maxticks - - ! be pessimistic - crystallite_integrateStress = .false. +!* be pessimistic - ! only integrate over fraction of timestep? - if (present(fraction)) then - dt = crystallite_subdt(g,i,e) * fraction - Fg_new = crystallite_subF0(:,:,g,i,e) + (crystallite_subF(:,:,g,i,e) - crystallite_subF0(:,:,g,i,e)) * fraction - else - dt = crystallite_subdt(g,i,e) - Fg_new = crystallite_subF(:,:,g,i,e) - endif +crystallite_integrateStress = .false. - ! feed local variables - Fp_current = crystallite_subFp0(:,:,g,i,e) - Tstar_v = crystallite_Tstar_v(:,g,i,e) - Lpguess_old = crystallite_Lp(:,:,g,i,e) ! consider present Lp good (i.e. worth remembering) ... - Lpguess = crystallite_Lp(:,:,g,i,e) ! ... and take it as first guess + +!* only integrate over fraction of timestep? + +if (present(fraction)) then + dt = crystallite_subdt(g,i,e) * fraction + Fg_new = crystallite_subF0(1:3,1:3,g,i,e) + (crystallite_subF(1:3,1:3,g,i,e) - crystallite_subF0(1:3,1:3,g,i,e)) * fraction +else + dt = crystallite_subdt(g,i,e) + Fg_new = crystallite_subF(1:3,1:3,g,i,e) +endif - ! inversion of Fp_current... - invFp_current = math_inv3x3(Fp_current) - if (all(invFp_current == 0.0_pReal)) then ! ... failed? - if (verboseDebugger .and. (e == debug_e .and. i == debug_i .and. g == debug_g)) then - !$OMP CRITICAL (write2out) - write(6,*) '::: integrateStress failed on invFp_current inversion',g,i,e - write(6,*) - write(6,'(a11,i3,x,i2,x,i5,/,3(3(f12.7,x)/))') 'invFp_new at ',g,i,e,math_transpose3x3(invFp_new(:,:)) - !$OMP END CRITICAL (write2out) - endif - return - endif +!* feed local variables + +Fp_current = crystallite_subFp0(1:3,1:3,g,i,e) +Tstar_v = crystallite_Tstar_v(1:6,g,i,e) +Lpguess_old = crystallite_Lp(1:3,1:3,g,i,e) ! consider present Lp good (i.e. worth remembering) ... +Lpguess = crystallite_Lp(1:3,1:3,g,i,e) ! ... and take it as first guess + + +!* inversion of Fp_current... + +invFp_current = math_inv3x3(Fp_current) +if (all(invFp_current == 0.0_pReal)) then ! ... failed? + if (verboseDebugger .and. (e == debug_e .and. i == debug_i .and. g == debug_g)) then + !$OMP CRITICAL (write2out) + write(6,*) '::: integrateStress failed on invFp_current inversion',g,i,e + write(6,*) + write(6,'(a11,i3,x,i2,x,i5,/,3(3(f12.7,x)/))') 'invFp_new at ',g,i,e,math_transpose3x3(invFp_new(1:3,1:3)) + !$OMP END CRITICAL (write2out) + endif + return +endif +A = math_mul33x33(transpose(invFp_current), math_mul33x33(transpose(Fg_new),math_mul33x33(Fg_new,invFp_current))) - A = math_mul33x33(transpose(invFp_current), math_mul33x33(transpose(Fg_new),math_mul33x33(Fg_new,invFp_current))) - ! get elasticity tensor - C_66 = constitutive_homogenizedC(g,i,e) -! if (debugger) write(6,'(a,/,6(6(f10.4,x)/))') 'elasticity',transpose(C_66(1:6,:))/1e9 - C = math_Mandel66to3333(C_66) - - ! start LpLoop with no acceleration - NiterationStress = 0_pInt - leapfrog = 1.0_pReal - maxleap = 16.0_pReal - jacoCounter = 0_pInt +!* get elasticity tensor + +C_66 = constitutive_homogenizedC(g,i,e) +! if (debugger) write(6,'(a,/,6(6(f10.4,x)/))') 'elasticity',transpose(C_66(1:6,1:6))/1e9 +C = math_Mandel66to3333(C_66) + + +!* start LpLoop with no acceleration + +NiterationStress = 0_pInt +leapfrog = 1.0_pReal +maxleap = 16.0_pReal +jacoCounter = 0_pInt LpLoop: do + NiterationStress = NiterationStress + 1 - ! increase loop counter - NiterationStress = NiterationStress + 1 + + !* too many loops required ? + + if (NiterationStress > nStress) then + if (verboseDebugger) then + !$OMP CRITICAL (write2out) + write(6,*) '::: integrateStress reached loop limit at ',g,i,e + write(6,*) + !$OMP END CRITICAL (write2out) + endif + return + endif - ! too many loops required ? - if (NiterationStress > nStress) then - if (verboseDebugger) then - !$OMP CRITICAL (write2out) - write(6,*) '::: integrateStress reached loop limit at ',g,i,e - write(6,*) - !$OMP END CRITICAL (write2out) - endif - return - endif + B = math_I3 - dt*Lpguess + BT = math_transpose3x3(B) + AB = math_mul33x33(A,B) + BTA = math_mul33x33(BT,A) - B = math_I3 - dt*Lpguess - BT = math_transpose3x3(B) - AB = math_mul33x33(A,B) - BTA = math_mul33x33(BT,A) + + !* calculate 2nd Piola-Kirchhoff stress tensor + + Tstar_v = 0.5_pReal * math_mul66x6(C_66,math_mandel33to6(math_mul33x33(BT,AB) - math_I3)) + p_hydro = sum(Tstar_v(1:3)) / 3.0_pReal + forall(n=1:3) Tstar_v(n) = Tstar_v(n) - p_hydro ! get deviatoric stress tensor - ! calculate 2nd Piola-Kirchhoff stress tensor - Tstar_v = 0.5_pReal*math_mul66x6(C_66,math_mandel33to6(math_mul33x33(BT,AB)-math_I3)) - p_hydro = sum(Tstar_v(1:3))/3.0_pReal - forall(n=1:3) Tstar_v(n) = Tstar_v(n) - p_hydro ! get deviatoric stress tensor + + !* calculate plastic velocity gradient and its tangent according to constitutive law + + call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) + call constitutive_LpAndItsTangent(Lp_constitutive, dLpdT_constitutive, Tstar_v, crystallite_Temperature(g,i,e), g, i, e) + call system_clock(count=tock,count_rate=tickrate,count_max=maxticks) + !$OMP CRITICAL (debugTimingLpTangent) + debug_cumLpCalls = debug_cumLpCalls + 1_pInt + debug_cumLpTicks = debug_cumLpTicks + tock-tick + if (tock < tick) debug_cumLpTicks = debug_cumLpTicks + maxticks + !$OMP END CRITICAL (debugTimingLpTangent) - ! calculate plastic velocity gradient and its tangent according to constitutive law - call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) - call constitutive_LpAndItsTangent(Lp_constitutive, dLpdT_constitutive, Tstar_v, crystallite_Temperature(g,i,e), g, i, e) - call system_clock(count=tock,count_rate=tickrate,count_max=maxticks) - !$OMP CRITICAL (debugTimingLpTangent) - debug_cumLpCalls = debug_cumLpCalls + 1_pInt - debug_cumLpTicks = debug_cumLpTicks + tock-tick - if (tock < tick) debug_cumLpTicks = debug_cumLpTicks + maxticks - !$OMP END CRITICAL (debugTimingLpTangent) - - if (verboseDebugger .and. (e == debug_e .and. i == debug_i .and. g == debug_g) .and. numerics_integrationMode == 1_pInt) then - !$OMP CRITICAL (write2out) - write(6,'(a,i3,x,i2,x,i5,x,a,x,i3)') '::: integrateStress at ' ,g,i,e, ' ; iteration ', NiterationStress - write(6,*) - write(6,'(a,/,3(3(e20.7,x)/))') 'Lp_constitutive', math_transpose3x3(Lp_constitutive(:,:)) - write(6,'(a,/,3(3(e20.7,x)/))') 'Lpguess', math_transpose3x3(Lpguess(:,:)) - !$OMP END CRITICAL (write2out) - endif + if (verboseDebugger .and. (e == debug_e .and. i == debug_i .and. g == debug_g) .and. numerics_integrationMode == 1_pInt) then + !$OMP CRITICAL (write2out) + write(6,'(a,i3,x,i2,x,i5,x,a,x,i3)') '::: integrateStress at ' ,g,i,e, ' ; iteration ', NiterationStress + write(6,*) + write(6,'(a,/,3(3(e20.7,x)/))') 'Lp_constitutive', math_transpose3x3(Lp_constitutive) + write(6,'(a,/,3(3(e20.7,x)/))') 'Lpguess', math_transpose3x3(Lpguess) + !$OMP END CRITICAL (write2out) + endif - ! update current residuum - residuum = Lpguess - Lp_constitutive - ! Check for convergence of loop - if (.not.(any(residuum/=residuum)) .and. & ! exclude any NaN in residuum + !* update current residuum and check for convergence of loop + + residuum = Lpguess - Lp_constitutive + if (.not.(any(residuum/=residuum)) .and. & ! exclude any NaN in residuum ( maxval(abs(residuum)) < aTol_crystalliteStress .or. & ! below absolute tolerance .or. ( any(abs(dt*Lpguess) > relevantStrain) .and. & ! worth checking? .and. maxval(abs(residuum/Lpguess), abs(dt*Lpguess) > relevantStrain) < rTol_crystalliteStress & ! below relative tolerance ) & ) & - ) & - exit LpLoop + ) then + exit LpLoop + endif + + + !* NaN occured at regular speed? -> return + + if (any(residuum/=residuum) .and. leapfrog == 1.0) then + if (debugger) then + !$OMP CRITICAL (write2out) + write(6,'(a,i3,x,i2,x,i5,x,a,i3,x,a)') '::: integrateStress encountered NaN at ',g,i,e,& + '; iteration ', NiterationStress, & + '>> returning..!' + !$OMP END CRITICAL (write2out) + endif + return + - ! NaN occured at regular speed? - if (any(residuum/=residuum) .and. leapfrog == 1.0) then - if (debugger) then - !$OMP CRITICAL (write2out) - write(6,'(a,i3,x,i2,x,i5,x,a,i3,x,a)') '::: integrateStress encountered NaN at ',g,i,e,& - '; iteration ', NiterationStress, & - '>> returning..!' - !$OMP END CRITICAL (write2out) - endif - return - - ! something went wrong at accelerated speed? - elseif (leapfrog > 1.0_pReal .and. & ! at fast pace .and. - ( sum(residuum*residuum) > sum(residuum_old*residuum_old) .or. & ! worse residuum .or. - sum(residuum*residuum_old) < 0.0_pReal .or. & ! residuum changed sign (overshoot) .or. - any(residuum/=residuum) & ! NaN occured - ) & - ) then - if (verboseDebugger) then - !$OMP CRITICAL (write2out) - write(6,'(a,i3,x,i2,x,i5,x,a,i3)') '::: integrateStress encountered high-speed crash at ',g,i,e,& - '; iteration ', NiterationStress - !$OMP END CRITICAL (write2out) - endif - maxleap = 0.5_pReal * leapfrog ! limit next acceleration - leapfrog = 1.0_pReal ! grinding halt - jacoCounter = 0_pInt ! reset counter for Jacobian update (we want to do an update next time!) - - ! restore old residuum and Lp - Lpguess = Lpguess_old - residuum = residuum_old - - !$OMP CRITICAL (distributionLeapfrogBreak) - debug_LeapfrogBreakDistribution(NiterationStress,numerics_integrationMode) = & - debug_LeapfrogBreakDistribution(NiterationStress,numerics_integrationMode) + 1 - !$OMP END CRITICAL (distributionLeapfrogBreak) - - ! residuum got better - else - ! calculate Jacobian for correction term - if (mod(jacoCounter, iJacoLpresiduum) == 0_pInt) then - dTdLp = 0.0_pReal - do h=1,3; do j=1,3; do k=1,3; do l=1,3; do m=1,3 -! forall (h=1:3,j=1:3,k=1:3,l=1:3,m=1:3) & - dTdLp(3*(h-1)+j,3*(k-1)+l) = dTdLp(3*(h-1)+j,3*(k-1)+l) + C(h,j,l,m)*AB(k,m)+C(h,j,m,l)*BTA(m,k) - enddo; enddo; enddo; enddo; enddo - dTdLp = -0.5_pReal*dt*dTdLp - dRdLp = math_identity2nd(9) - math_mul99x99(dLpdT_constitutive,dTdLp) - invdRdLp = 0.0_pReal - call math_invert(9,dRdLp,invdRdLp,dummy,error) ! invert dR/dLp --> dLp/dR - if (error) then - if (verboseDebugger .and. (e == debug_e .and. i == debug_i .and. g == debug_g)) then - !$OMP CRITICAL (write2out) - write(6,'(a,i3,x,i2,x,i5,x,a,i3)') '::: integrateStress failed on dR/dLp inversion at ',g,i,e, & - '; iteration ', NiterationStress - write(6,*) - write(6,'(a,/,9(9(e15.3,x)/))') 'dRdLp',transpose(dRdLp(:,:)) - write(6,'(a,/,9(9(e15.3,x)/))') 'dLpdT_constitutive',transpose(dLpdT_constitutive(:,:)) - write(6,'(a,/,3(3(e20.7,x)/))') 'Lp_constitutive',math_transpose3x3(Lp_constitutive(:,:)) - write(6,'(a,/,3(3(e20.7,x)/))') 'Lpguess',math_transpose3x3(Lpguess(:,:)) - !$OMP END CRITICAL (write2out) - endif - return - else - if (verboseDebugger .and. (e == debug_e .and. i == debug_i .and. g == debug_g) .and. numerics_integrationMode==1_pInt) then - !$OMP CRITICAL (write2out) - write(6,'(a,i3,x,i2,x,i5,x,a,i3)') '::: integrateStress did dR/dLp inversion at ',g,i,e, & - '; iteration ', NiterationStress - write(6,*) - write(6,'(a,/,9(9(e15.3,x)/))') 'dRdLp',transpose(dRdLp(:,:)) - write(6,'(a,/,9(9(e15.3,x)/))') 'dLpdT_constitutive',transpose(dLpdT_constitutive(:,:)) - !$OMP END CRITICAL (write2out) - endif + !* something went wrong at accelerated speed? -> restore old residuum and Lp + + elseif (leapfrog > 1.0_pReal .and. & ! at fast pace .and. + ( sum(residuum*residuum) > sum(residuum_old*residuum_old) .or. & ! worse residuum .or. + sum(residuum*residuum_old) < 0.0_pReal .or. & ! residuum changed sign (overshoot) .or. + any(residuum/=residuum) & ! NaN occured + ) & + ) then + if (verboseDebugger) then + !$OMP CRITICAL (write2out) + write(6,'(a,i3,x,i2,x,i5,x,a,i3)') '::: integrateStress encountered high-speed crash at ',g,i,e,& + '; iteration ', NiterationStress + !$OMP END CRITICAL (write2out) + endif + maxleap = 0.5_pReal * leapfrog ! limit next acceleration + leapfrog = 1.0_pReal ! grinding halt + jacoCounter = 0_pInt ! reset counter for Jacobian update (we want to do an update next time!) + Lpguess = Lpguess_old + residuum = residuum_old + !$OMP CRITICAL (distributionLeapfrogBreak) + debug_LeapfrogBreakDistribution(NiterationStress,numerics_integrationMode) = & + debug_LeapfrogBreakDistribution(NiterationStress,numerics_integrationMode) + 1 + !$OMP END CRITICAL (distributionLeapfrogBreak) + + + !* residuum got better? -> calculate Jacobian for correction term and remember current residuum and Lpguess + + else + if (mod(jacoCounter, iJacoLpresiduum) == 0_pInt) then + dTdLp = 0.0_pReal + do h=1,3; do j=1,3; do k=1,3; do l=1,3; do m=1,3 + dTdLp(3*(h-1)+j,3*(k-1)+l) = dTdLp(3*(h-1)+j,3*(k-1)+l) + C(h,j,l,m)*AB(k,m)+C(h,j,m,l)*BTA(m,k) + enddo; enddo; enddo; enddo; enddo + dTdLp = -0.5_pReal*dt*dTdLp + dRdLp = math_identity2nd(9) - math_mul99x99(dLpdT_constitutive,dTdLp) + invdRdLp = 0.0_pReal + call math_invert(9,dRdLp,invdRdLp,dummy,error) ! invert dR/dLp --> dLp/dR + if (error) then + if (verboseDebugger .and. (e == debug_e .and. i == debug_i .and. g == debug_g)) then + !$OMP CRITICAL (write2out) + write(6,'(a,i3,x,i2,x,i5,x,a,i3)') '::: integrateStress failed on dR/dLp inversion at ',g,i,e, & + '; iteration ', NiterationStress + write(6,*) + write(6,'(a,/,9(9(e15.3,x)/))') 'dRdLp',transpose(dRdLp) + write(6,'(a,/,9(9(e15.3,x)/))') 'dLpdT_constitutive',transpose(dLpdT_constitutive) + write(6,'(a,/,3(3(e20.7,x)/))') 'Lp_constitutive',math_transpose3x3(Lp_constitutive) + write(6,'(a,/,3(3(e20.7,x)/))') 'Lpguess',math_transpose3x3(Lpguess) + !$OMP END CRITICAL (write2out) endif - endif - jacoCounter = jacoCounter + 1_pInt ! increase counter for jaco update + return + else + if (verboseDebugger .and. (e == debug_e .and. i == debug_i .and. g == debug_g) .and. numerics_integrationMode==1_pInt) then + !$OMP CRITICAL (write2out) + write(6,'(a,i3,x,i2,x,i5,x,a,i3)') '::: integrateStress did dR/dLp inversion at ',g,i,e, & + '; iteration ', NiterationStress + write(6,*) + write(6,'(a,/,9(9(e15.3,x)/))') 'dRdLp',transpose(dRdLp(:,:)) + write(6,'(a,/,9(9(e15.3,x)/))') 'dLpdT_constitutive',transpose(dLpdT_constitutive) + !$OMP END CRITICAL (write2out) + endif + endif + endif + jacoCounter = jacoCounter + 1_pInt ! increase counter for jaco update + residuum_old = residuum + Lpguess_old = Lpguess - ! remember current residuum and Lpguess - residuum_old = residuum - Lpguess_old = Lpguess - - ! accelerate? - if (NiterationStress > 1 .and. leapfrog+1.0_pReal <= maxleap) leapfrog = leapfrog + 1.0_pReal - endif + + !* accelerate? + + if (NiterationStress > 1 .and. leapfrog+1.0_pReal <= maxleap) leapfrog = leapfrog + 1.0_pReal + + endif - ! leapfrog to updated Lp - do k=1,3; do l=1,3; do m=1,3; do n=1,3 - Lpguess(k,l) = Lpguess(k,l) - leapfrog*invdRdLp(3*(k-1)+l,3*(m-1)+n)*residuum(m,n) - enddo; enddo; enddo; enddo - enddo LpLoop + + !* leapfrog to updated Lp + + do k=1,3; do l=1,3; do m=1,3; do n=1,3 + Lpguess(k,l) = Lpguess(k,l) - leapfrog * invdRdLp(3*(k-1)+l,3*(m-1)+n) * residuum(m,n) + enddo; enddo; enddo; enddo - ! calculate new plastic and elastic deformation gradient - invFp_new = math_mul33x33(invFp_current,B) - invFp_new = invFp_new/math_det3x3(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize by det - call math_invert3x3(invFp_new,Fp_new,det,error) - if (error) then - if (verboseDebugger .and. (e == debug_e .and. i == debug_i .and. g == debug_g)) then - !$OMP CRITICAL (write2out) - write(6,'(a,i3,x,i2,x,i5,x,a,x,i3)') '::: integrateStress failed on invFp_new inversion at ',g,i,e, & - ' ; iteration ', NiterationStress - write(6,*) - write(6,'(a11,3(i3,x),/,3(3(f12.7,x)/))') 'invFp_new at ',g,i,e,math_transpose3x3(invFp_new(:,:)) - !$OMP END CRITICAL (write2out) - endif - return - endif - Fe_new = math_mul33x33(Fg_new,invFp_new) ! calc resulting Fe +enddo LpLoop - ! add volumetric component to 2nd Piola-Kirchhoff stress - forall (n=1:3) Tstar_v(n) = Tstar_v(n) + p_hydro + +!* calculate new plastic and elastic deformation gradient + +invFp_new = math_mul33x33(invFp_current,B) +invFp_new = invFp_new/math_det3x3(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize by det +call math_invert3x3(invFp_new,Fp_new,det,error) +if (error) then + if (verboseDebugger .and. (e == debug_e .and. i == debug_i .and. g == debug_g)) then + !$OMP CRITICAL (write2out) + write(6,'(a,i3,x,i2,x,i5,x,a,x,i3)') '::: integrateStress failed on invFp_new inversion at ',g,i,e, & + ' ; iteration ', NiterationStress + write(6,*) + write(6,'(a11,3(i3,x),/,3(3(f12.7,x)/))') 'invFp_new at ',g,i,e,math_transpose3x3(invFp_new) + !$OMP END CRITICAL (write2out) + endif + return +endif +Fe_new = math_mul33x33(Fg_new,invFp_new) ! calc resulting Fe + + +!* add volumetric component to 2nd Piola-Kirchhoff stress and calculate 1st Piola-Kirchhoff stress + +forall (n=1:3) Tstar_v(n) = Tstar_v(n) + p_hydro +crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(Fe_new, math_mul33x33(math_Mandel6to33(Tstar_v), math_transpose3x3(invFp_new))) - ! calculate 1st Piola-Kirchhoff stress - crystallite_P(:,:,g,i,e) = math_mul33x33(Fe_new,math_mul33x33(math_Mandel6to33(Tstar_v),transpose(invFp_new))) - - ! store local values in global variables - crystallite_Lp(:,:,g,i,e) = Lpguess - crystallite_Tstar_v(:,g,i,e) = Tstar_v - crystallite_Fp(:,:,g,i,e) = Fp_new - crystallite_Fe(:,:,g,i,e) = Fe_new - crystallite_invFp(:,:,g,i,e) = invFp_new - ! set return flag to true - crystallite_integrateStress = .true. - if (verboseDebugger .and. (e == debug_e .and. i == debug_i .and. g == debug_g) .and. numerics_integrationMode == 1_pInt) then - !$OMP CRITICAL (write2out) - write(6,'(a,i3,x,i2,x,i5,x,a,x,i3)') '::: integrateStress converged at ',g,i,e,' ; iteration ', NiterationStress - write(6,*) - write(6,'(a,/,3(3(f12.7,x)/))') 'P / MPa',math_transpose3x3(crystallite_P(1:3,:,g,i,e))/1e6 - write(6,'(a,/,3(3(f12.7,x)/))') 'Cauchy / MPa', math_mul33x33(crystallite_P(1:3,1:3,g,i,e),transpose(Fg_new)) & - / 1e6 / math_det3x3(Fg_new) - write(6,'(a,/,3(3(f12.7,x)/))') 'Fe Lp Fe^-1',math_transpose3x3( & - math_mul33x33(Fe_new, math_mul33x33(crystallite_Lp(1:3,1:3,g,i,e), & - math_inv3x3(Fe_new)))) ! transpose to get correct print out order - write(6,'(a,/,3(3(f12.7,x)/))') 'Fp',math_transpose3x3(crystallite_Fp(:,:,g,i,e)) - !$OMP END CRITICAL (write2out) - endif +!* store local values in global variables - !$OMP CRITICAL (distributionStress) - debug_StressLoopDistribution(NiterationStress,numerics_integrationMode) = & - debug_StressLoopDistribution(NiterationStress,numerics_integrationMode) + 1 - !$OMP END CRITICAL (distributionStress) +crystallite_Lp(1:3,1:3,g,i,e) = Lpguess +crystallite_Tstar_v(1:6,g,i,e) = Tstar_v +crystallite_Fp(1:3,1:3,g,i,e) = Fp_new +crystallite_Fe(1:3,1:3,g,i,e) = Fe_new +crystallite_invFp(1:3,1:3,g,i,e) = invFp_new - return + +!* set return flag to true + +crystallite_integrateStress = .true. +if (verboseDebugger .and. (e == debug_e .and. i == debug_i .and. g == debug_g) .and. numerics_integrationMode == 1_pInt) then + !$OMP CRITICAL (write2out) + write(6,'(a,i3,x,i2,x,i5,x,a,x,i3)') '::: integrateStress converged at ',g,i,e,' ; iteration ', NiterationStress + write(6,*) + write(6,'(a,/,3(3(f12.7,x)/))') 'P / MPa',math_transpose3x3(crystallite_P(1:3,1:3,g,i,e))/1e6 + write(6,'(a,/,3(3(f12.7,x)/))') 'Cauchy / MPa', math_mul33x33(crystallite_P(1:3,1:3,g,i,e), math_transpose3x3(Fg_new)) & + / 1e6 / math_det3x3(Fg_new) + write(6,'(a,/,3(3(f12.7,x)/))') 'Fe Lp Fe^-1',math_transpose3x3( & + math_mul33x33(Fe_new, math_mul33x33(crystallite_Lp(1:3,1:3,g,i,e), & + math_inv3x3(Fe_new)))) ! transpose to get correct print out order + write(6,'(a,/,3(3(f12.7,x)/))') 'Fp',math_transpose3x3(crystallite_Fp(1:3,1:3,g,i,e)) + !$OMP END CRITICAL (write2out) +endif + +!$OMP CRITICAL (distributionStress) + debug_StressLoopDistribution(NiterationStress,numerics_integrationMode) = & + debug_StressLoopDistribution(NiterationStress,numerics_integrationMode) + 1 +!$OMP END CRITICAL (distributionStress) endfunction @@ -2690,8 +2895,7 @@ use mesh, only: mesh_element, & FE_NipNeighbors use debug, only: debugger, & debug_e, debug_i, debug_g, & - verboseDebugger, & - selectiveDebugger + verboseDebugger use constitutive_nonlocal, only: constitutive_nonlocal_structure, & constitutive_nonlocal_updateCompatibility @@ -2715,11 +2919,12 @@ integer(pInt) e, & ! element index myStructure, & ! lattice structure neighboringStructure real(pReal), dimension(3,3) :: U, R +real(pReal), dimension(4) :: orientation logical error ! --- CALCULATE ORIENTATION AND LATTICE ROTATION --- -!$OMP PARALLEL DO PRIVATE(error,U,R) +!$OMP PARALLEL DO PRIVATE(error,U,R,orientation) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do g = 1,homogenization_Ngrains(mesh_element(3,e)) @@ -2727,16 +2932,14 @@ logical error call math_pDecomposition(crystallite_Fe(1:3,1:3,g,i,e), U, R, error) ! polar decomposition of Fe if (error) then call IO_warning(650, e, i, g) - crystallite_orientation(:,g,i,e) = (/1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal/) ! fake orientation + orientation = (/1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal/) ! fake orientation else - crystallite_orientation(:,g,i,e) = math_RtoQuaternion(transpose(R)) + orientation = math_RtoQuaternion(transpose(R)) endif - - crystallite_rotation(:,g,i,e) = & - math_QuaternionDisorientation( math_qConj(crystallite_orientation(1:4,g,i,e)), & ! calculate grainrotation - math_qConj(crystallite_orientation0(1:4,g,i,e)), & - 0_pInt ) ! we don't want symmetry here - + crystallite_rotation(1:4,g,i,e) = math_QuaternionDisorientation(math_qConj(orientation), & + math_qConj(orientation), & + 0_pInt ) ! we don't want symmetry here + crystallite_orientation(1:4,g,i,e) = orientation enddo enddo enddo @@ -2845,64 +3048,65 @@ function crystallite_postResults(& crystallite_postResults = 0.0_pReal c = 0_pInt - crystallite_postResults(c+1) = crystallite_sizePostResults(crystID); c = c+1_pInt ! size of results from cryst + crystallite_postResults(c+1) = crystallite_sizePostResults(crystID) ! size of results from cryst + c = c + 1_pInt do o = 1,crystallite_Noutput(crystID) select case(crystallite_output(o,crystID)) case ('phase') - crystallite_postResults(c+1) = material_phase(g,i,e) ! phaseID of grain + crystallite_postResults(c+1) = material_phase(g,i,e) ! phaseID of grain c = c + 1_pInt case ('volume') - crystallite_postResults(c+1) = material_volume(g,i,e) ! grain volume (not fraction but absolute, right?) + crystallite_postResults(c+1) = material_volume(g,i,e) ! grain volume (not fraction but absolute, right?) c = c + 1_pInt case ('orientation') - crystallite_postResults(c+1:c+4) = & - crystallite_orientation(:,g,i,e) ! grain orientation as quaternion + crystallite_postResults(c+1:c+4) = crystallite_orientation(1:4,g,i,e) ! grain orientation as quaternion c = c + 4_pInt case ('eulerangles') - crystallite_postResults(c+1:c+3) = inDeg * & - math_QuaternionToEuler(crystallite_orientation(:,g,i,e)) ! grain orientation as Euler angles in degree + crystallite_postResults(c+1:c+3) = inDeg * math_QuaternionToEuler(crystallite_orientation(1:4,g,i,e)) ! grain orientation as Euler angles in degree c = c + 3_pInt case ('grainrotation') - crystallite_postResults(c+1:c+4) = & - math_QuaternionToAxisAngle(crystallite_rotation(1:4,g,i,e)) ! grain rotation away from initial orientation as axis-angle - crystallite_postResults(c+4) = inDeg * crystallite_postResults(c+4) ! angle in degree + crystallite_postResults(c+1:c+4) = math_QuaternionToAxisAngle(crystallite_rotation(1:4,g,i,e)) ! grain rotation away from initial orientation as axis-angle + crystallite_postResults(c+4) = inDeg * crystallite_postResults(c+4) ! angle in degree c = c + 4_pInt case ('defgrad','f') mySize = 9_pInt - crystallite_postResults(c+1:c+1+mySize) = reshape(math_transpose3x3(crystallite_partionedF(:,:,g,i,e)),(/mySize/)) + crystallite_postResults(c+1:c+1+mySize) = reshape(math_transpose3x3(crystallite_partionedF(1:3,1:3,g,i,e)),(/mySize/)) c = c + mySize case ('fe') mySize = 9_pInt - crystallite_postResults(c+1:c+1+mySize) = reshape(math_transpose3x3(crystallite_Fe(:,:,g,i,e)),(/mySize/)) + crystallite_postResults(c+1:c+1+mySize) = reshape(math_transpose3x3(crystallite_Fe(1:3,1:3,g,i,e)),(/mySize/)) c = c + mySize case ('ee') - Ee = 0.5_pReal * (math_mul33x33(math_transpose3x3(crystallite_Fe(:,:,g,i,e)), crystallite_Fe(:,:,g,i,e)) - math_I3) + Ee = 0.5_pReal * (math_mul33x33(math_transpose3x3(crystallite_Fe(1:3,1:3,g,i,e)), crystallite_Fe(1:3,1:3,g,i,e)) - math_I3) mySize = 9_pInt - crystallite_postResults(c+1:c+1+mySize) = reshape(Ee(:,:),(/mySize/)) + crystallite_postResults(c+1:c+1+mySize) = reshape(Ee,(/mySize/)) c = c + mySize case ('fp') mySize = 9_pInt - crystallite_postResults(c+1:c+1+mySize) = reshape(math_transpose3x3(crystallite_Fp(:,:,g,i,e)),(/mySize/)) + crystallite_postResults(c+1:c+1+mySize) = reshape(math_transpose3x3(crystallite_Fp(1:3,1:3,g,i,e)),(/mySize/)) c = c + mySize case ('lp') mySize = 9_pInt - crystallite_postResults(c+1:c+1+mySize) = reshape(math_transpose3x3(crystallite_Lp(:,:,g,i,e)),(/mySize/)) + crystallite_postResults(c+1:c+1+mySize) = reshape(math_transpose3x3(crystallite_Lp(1:3,1:3,g,i,e)),(/mySize/)) c = c + mySize case ('p','firstpiola','1stpiola') mySize = 9_pInt - crystallite_postResults(c+1:c+1+mySize) = reshape(math_transpose3x3(crystallite_P(:,:,g,i,e)),(/mySize/)) + crystallite_postResults(c+1:c+1+mySize) = reshape(math_transpose3x3(crystallite_P(1:3,1:3,g,i,e)),(/mySize/)) c = c + mySize case ('s','tstar','secondpiola','2ndpiola') mySize = 9_pInt - crystallite_postResults(c+1:c+1+mySize) = reshape(math_Mandel6to33(crystallite_Tstar_v(:,g,i,e)),(/mySize/)) + crystallite_postResults(c+1:c+1+mySize) = reshape(math_Mandel6to33(crystallite_Tstar_v(1:6,g,i,e)),(/mySize/)) c = c + mySize end select enddo - crystallite_postResults(c+1) = constitutive_sizePostResults(g,i,e); c = c+1_pInt ! size of constitutive results - crystallite_postResults(c+1:c+constitutive_sizePostResults(g,i,e)) = & - constitutive_postResults(crystallite_Tstar_v(:,g,i,e), crystallite_Fe(:,:,g,i,e), crystallite_Temperature(g,i,e), dt, g, i, e) + crystallite_postResults(c+1) = constitutive_sizePostResults(g,i,e) ! size of constitutive results + c = c + 1_pInt + crystallite_postResults(c+1:c+constitutive_sizePostResults(g,i,e)) = constitutive_postResults(crystallite_Tstar_v(1:6,g,i,e), & + crystallite_Fe(1:3,1:3,g,i,e), & + crystallite_Temperature(g,i,e), & + dt, g, i, e) c = c + constitutive_sizePostResults(g,i,e) return diff --git a/code/debug.f90 b/code/debug.f90 index 882dbdc29..e7efbdf98 100644 --- a/code/debug.f90 +++ b/code/debug.f90 @@ -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 diff --git a/code/homogenization.f90 b/code/homogenization.f90 index a339f08cd..6715a0183 100644 --- a/code/homogenization.f90 +++ b/code/homogenization.f90 @@ -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 <<>> + 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 <<>> + 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 <<>> +!* 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 ! <> @@ -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)) ! <> + !$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 - ! <> + !$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 <<>> +!* 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 <<>> + 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 <<>> +!* 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 <<>> use homogenization_isostrain - use homogenization_RGC ! RGC homogenization added <<>> + 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 <<>> + 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 <<>> + 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 <<>> + 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 <<>> + 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 <<>> + 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 <<>> + 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 <<>> +!* RGC homogenization case (homogenization_RGC_label) homogenization_postResults = homogenization_RGC_postResults(homogenization_state(ip,el),ip,el) end select diff --git a/code/homogenization_RGC.f90 b/code/homogenization_RGC.f90 index 9fae5cf28..52ae8aae6 100644 --- a/code/homogenization_RGC.f90 +++ b/code/homogenization_RGC.f90 @@ -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 diff --git a/code/homogenization_isostrain.f90 b/code/homogenization_isostrain.f90 index 03e6458b1..50c9e578d 100644 --- a/code/homogenization_isostrain.f90 +++ b/code/homogenization_isostrain.f90 @@ -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 <<>> + 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 diff --git a/code/lattice.f90 b/code/lattice.f90 index 9db7c60d7..61338a213 100644 --- a/code/lattice.f90 +++ b/code/lattice.f90 @@ -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 diff --git a/code/material.f90 b/code/material.f90 index a8f17b8d4..7da65e97a 100644 --- a/code/material.f90 +++ b/code/material.f90 @@ -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) diff --git a/code/math.f90 b/code/math.f90 index 5a98ec0c3..bb1c5e4e0 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -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 diff --git a/code/mesh.f90 b/code/mesh.f90 index eee8c7db6..5f54e10f3 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -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... diff --git a/code/mpie_cpfem_marc.f90 b/code/mpie_cpfem_marc.f90 index af56dfd9d..27257e9cf 100644 --- a/code/mpie_cpfem_marc.f90 +++ b/code/mpie_cpfem_marc.f90 @@ -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 diff --git a/code/numerics.f90 b/code/numerics.f90 index bee72aad0..0d9df4b32 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -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 diff --git a/code/prec.f90 b/code/prec.f90 index a3b21a273..90579e7a3 100644 --- a/code/prec.f90 +++ b/code/prec.f90 @@ -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 diff --git a/code/prec_single.f90 b/code/prec_single.f90 index 428fe59a9..cedceb0a1 100644 --- a/code/prec_single.f90 +++ b/code/prec_single.f90 @@ -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