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