From b25396374a983781d9b58bc1598047b4b1abe7fb Mon Sep 17 00:00:00 2001 From: Denny Tjahjanto Date: Thu, 27 Aug 2009 12:10:06 +0000 Subject: [PATCH] homogenization_RGC.f90 >>> adding some lines, mostly for debugging purpose. no critical change. constitutive_phenopowerlaw.f90 >>> adding new parameter: constitutive_phenopowerlaw_w0_slip, i.e., the hardening rate exponent. homogenization.f90 >>> most important change is to add an if-else statement (line 379-380) to switch crystallite_requeted = .false. for already converged material point iteration (el/ip). the rest of the changes are cosmetics and debugging stuffs. crystallite.f90 >>> similar to homogenization.f90, the most important change is to add additional if-else statement (line 574) in the jacobian (perturbation) loop. now the jacobian calculation will only be performed when crystallite_requested = .true.. the rest is only cosmetic. --- code/constitutive_phenopowerlaw.f90 | 28 ++- code/crystallite.f90 | 184 +++++++++--------- code/homogenization.f90 | 39 ++-- code/homogenization_RGC.f90 | 289 ++++++++++++++++++---------- 4 files changed, 325 insertions(+), 215 deletions(-) diff --git a/code/constitutive_phenopowerlaw.f90 b/code/constitutive_phenopowerlaw.f90 index f09661727..3436739ef 100644 --- a/code/constitutive_phenopowerlaw.f90 +++ b/code/constitutive_phenopowerlaw.f90 @@ -103,8 +103,9 @@ MODULE constitutive_phenopowerlaw real(pReal), dimension(:,:,:), allocatable :: constitutive_phenopowerlaw_hardeningMatrix_slipslip real(pReal), dimension(:,:,:), allocatable :: constitutive_phenopowerlaw_hardeningMatrix_sliptwin real(pReal), dimension(:,:,:), allocatable :: constitutive_phenopowerlaw_hardeningMatrix_twinslip - real(pReal), dimension(:,:,:), allocatable :: constitutive_phenopowerlaw_hardeningMatrix_twintwin + real(pReal), dimension(:,:,:), allocatable :: constitutive_phenopowerlaw_hardeningMatrix_twintwin + real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_w0_slip CONTAINS !**************************************** @@ -203,6 +204,8 @@ subroutine constitutive_phenopowerlaw_init(file) constitutive_phenopowerlaw_interaction_twinslip = 0.0_pReal constitutive_phenopowerlaw_interaction_twintwin = 0.0_pReal + allocate(constitutive_phenopowerlaw_w0_slip(maxNinstance)) ; constitutive_phenopowerlaw_w0_slip = 0.0_pReal + rewind(file) line = '' section = 0 @@ -288,25 +291,39 @@ subroutine constitutive_phenopowerlaw_init(file) case ('interaction_twintwin') forall (j = 1:lattice_maxNinteraction) & constitutive_phenopowerlaw_interaction_twintwin(j,i) = IO_floatValue(line,positions,1+j) + case ('w0_slip') + constitutive_phenopowerlaw_w0_slip(i) = IO_floatValue(line,positions,2) end select endif enddo + write(6,*)'constitutive_phenopowerlaw_structureName = ',constitutive_phenopowerlaw_structureName + write(6,*)'constitutive_phenopowerlaw_structure = ',constitutive_phenopowerlaw_structure + 100 do i = 1,maxNinstance + write(6,*)'1constitutive_phenopowerlaw_gdot0_slip = ',constitutive_phenopowerlaw_gdot0_slip(i) + write(6,*)'1constitutive_phenopowerlaw_w0_slip = ',constitutive_phenopowerlaw_w0_slip(i) constitutive_phenopowerlaw_structure(i) = lattice_initializeStructure(constitutive_phenopowerlaw_structureName(i), & ! get structure constitutive_phenopowerlaw_CoverA(i)) + write(6,*)'>constitutive_phenopowerlaw_structure = ',constitutive_phenopowerlaw_structure(i) + write(6,*)'2constitutive_phenopowerlaw_gdot0_slip = ',constitutive_phenopowerlaw_gdot0_slip(i) constitutive_phenopowerlaw_Nslip(:,i) = min(lattice_NslipSystem(:,constitutive_phenopowerlaw_structure(i)),& ! limit active slip systems per family to max constitutive_phenopowerlaw_Nslip(:,i)) + write(6,*)'3constitutive_phenopowerlaw_gdot0_slip = ',constitutive_phenopowerlaw_gdot0_slip(i) constitutive_phenopowerlaw_Ntwin(:,i) = min(lattice_NtwinSystem(:,constitutive_phenopowerlaw_structure(i)),& ! limit active twin systems per family to max constitutive_phenopowerlaw_Ntwin(:,i)) + write(6,*)'4constitutive_phenopowerlaw_gdot0_slip = ',constitutive_phenopowerlaw_gdot0_slip(i) constitutive_phenopowerlaw_totalNslip(i) = sum(constitutive_phenopowerlaw_Nslip(:,i)) ! how many slip systems altogether + write(6,*)'5constitutive_phenopowerlaw_gdot0_slip = ',constitutive_phenopowerlaw_gdot0_slip(i) constitutive_phenopowerlaw_totalNtwin(i) = sum(constitutive_phenopowerlaw_Ntwin(:,i)) ! how many twin systems altogether + write(6,*)'6constitutive_phenopowerlaw_gdot0_slip = ',constitutive_phenopowerlaw_gdot0_slip(i) - if (constitutive_phenopowerlaw_structure(i) < 1 .or. & ! sanity checks - constitutive_phenopowerlaw_structure(i) > 3) call IO_error(205) - + if (constitutive_phenopowerlaw_structure(i) < 1 ) call IO_error(205) + write(6,*)'7constitutive_phenopowerlaw_gdot0_slip = ',constitutive_phenopowerlaw_gdot0_slip(i) if (any(constitutive_phenopowerlaw_tau0_slip(:,i) < 0.0_pReal .and. & constitutive_phenopowerlaw_Nslip(:,i) > 0)) call IO_error(210) + write(6,*)'8constitutive_phenopowerlaw_gdot0_slip = ',constitutive_phenopowerlaw_gdot0_slip(i) + write(6,*)'8constitutive_phenopowerlaw_w0_slip = ',constitutive_phenopowerlaw_w0_slip(i) if (constitutive_phenopowerlaw_gdot0_slip(i) <= 0.0_pReal) call IO_error(211) if (constitutive_phenopowerlaw_n_slip(i) <= 0.0_pReal) call IO_error(212) if (any(constitutive_phenopowerlaw_tausat_slip(:,i) <= 0.0_pReal .and. & @@ -704,7 +721,8 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el constitutive_phenopowerlaw_spr(matID)*dsqrt(state(ipc,ip,el)%p(index_F)) 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)/ssat) ! system-dependent prefactor for slip--slip interaction + h_slipslip(j) = c_slipslip*(1.0_pReal-state(ipc,ip,el)%p(j)/ssat)**constitutive_phenopowerlaw_w0_slip(matID) + ! system-dependent prefactor for slip--slip interaction h_sliptwin(j) = c_sliptwin ! no system-dependent part diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 1a05340e8..0216332e1 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -196,7 +196,7 @@ subroutine crystallite_init(Temperature) write(6,'(a35,x,7(i5,x))') 'crystallite_subdt: ', shape(crystallite_subdt) write(6,'(a35,x,7(i5,x))') 'crystallite_subFrac: ', shape(crystallite_subFrac) write(6,'(a35,x,7(i5,x))') 'crystallite_subStep: ', shape(crystallite_subStep) - write(6,'(a35,x,7(i5,x))') 'crystallite_localConstitution: ', shape(crystallite_localConstitution) + write(6,'(a35,x,7(i5,x))') 'crystallite_localConstitution: ', shape(crystallite_localConstitution) write(6,'(a35,x,7(i5,x))') 'crystallite_requested: ', shape(crystallite_requested) write(6,'(a35,x,7(i5,x))') 'crystallite_onTrack: ', shape(crystallite_onTrack) write(6,'(a35,x,7(i5,x))') 'crystallite_converged: ', shape(crystallite_converged) @@ -239,7 +239,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) math_mul66x6, & math_Mandel6to33, & math_Mandel33to6, & - math_I3 + math_I3, & + math_Plain3333to99 use FEsolving, only: FEsolving_execElem, & FEsolving_execIP, & theInc @@ -292,6 +293,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) stateConverged, & ! flag indicating if state converged converged ! flag indicating if iteration converged + real(pReal), dimension(9,9) :: dPdF99 ! ------ initialize to starting condition ------ @@ -347,7 +349,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains - debugger = (e == 1 .and. i == 1 .and. g == 1) + debugger = .false. !(e == 1 .and. i == 1 .and. g == 1) if (crystallite_converged(g,i,e)) then if (debugger) then !$OMP CRITICAL (write2out) @@ -395,12 +397,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) crystallite_subStep(g,i,e) * & (crystallite_partionedF(:,:,g,i,e) - crystallite_partionedF0(:,:,g,i,e)) crystallite_subdt(g,i,e) = crystallite_subStep(g,i,e) * crystallite_dt(g,i,e) - if (debugger) then - !$OMP CRITICAL (write2out) - write(6,'(a36,e8.3)') 'current timestep crystallite_subdt: ',crystallite_subdt(g,i,e) - write(6,*) - !$OMPEND CRITICAL (write2out) - endif crystallite_converged(g,i,e) = .false. ! start out non-converged endif enddo @@ -505,7 +501,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains - debugger = (e == 1 .and. i == 1 .and. g == 1) if (crystallite_nonfinished(g,i,e)) & ! all undone crystallites call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_Fp(:,:,g,i,e), & crystallite_invFp(:,:,g,i,e), crystallite_Temperature(g,i,e), g, i, e) @@ -516,7 +511,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains - debugger = (e == 1 .and. i == 1 .and. g == 1) if (crystallite_nonfinished(g,i,e)) then ! all undone crystallites crystallite_stateConverged(g,i,e) = crystallite_updateState(g,i,e) ! update state crystallite_temperatureConverged(g,i,e) = crystallite_updateTemperature(g,i,e) ! update temperature @@ -533,7 +527,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo !$OMPEND PARALLEL DO - write(6,*) 'NiterationState: ',NiterationState +! write(6,*) 'NiterationState: ',NiterationState crystallite_nonfinished = crystallite_nonfinished .and. crystallite_onTrack .and. .not. crystallite_converged enddo ! crystallite convergence loop @@ -552,8 +546,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains - if (.not. crystallite_converged(g,i,e)) then ! respond fully elastically - call IO_warning(600,e,i,g) + if (.not. crystallite_converged(g,i,e)) then ! respond fully elastically (might be not required due to becoming terminally ill anyway) +! call IO_warning(600,e,i,g) invFp = math_inv3x3(crystallite_partionedFp0(:,:,g,i,e)) Fe_guess = math_mul33x33(crystallite_partionedF(:,:,g,i,e),invFp) Tstar = math_Mandel6to33( & @@ -576,95 +570,95 @@ subroutine crystallite_stressAndItsTangent(updateJaco) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains -! debugger = (g == 1 .and. i == 1 .and. e == 1) - if (crystallite_converged(g,i,e)) then ! grain converged in above iteration - mySizeState = constitutive_sizeState(g,i,e) ! number of state variables for this grain - mySizeDotState = constitutive_sizeDotState(g,i,e) ! number of dotStates for this grain - myState(1:mySizeState) = constitutive_state(g,i,e)%p ! remember unperturbed, converged state, ... - myDotState(1:mySizeDotState) = constitutive_dotState(g,i,e)%p ! ... dotStates, ... - myTemperature = crystallite_Temperature(g,i,e) ! ... Temperature, ... - myF = crystallite_subF(:,:,g,i,e) ! ... and kinematics - myFp = crystallite_Fp(:,:,g,i,e) - myInvFp = crystallite_invFp(:,:,g,i,e) - myFe = crystallite_Fe(:,:,g,i,e) - myLp = crystallite_Lp(:,:,g,i,e) - myTstar_v = crystallite_Tstar_v(:,g,i,e) - myP = crystallite_P(:,:,g,i,e) - if (debugger) then - !$OMP CRITICAL (write2out) - write (6,*) '#############' - write (6,*) 'central solution' - write (6,*) '#############' - write (6,'(a,/,3(3(f12.4,x)/))') ' P of 1 1 1',myP(1:3,:)/1e6 - write (6,'(a,/,3(3(f12.8,x)/))') ' Fp of 1 1 1',myFp(1:3,:) - write (6,'(a,/,3(3(f12.8,x)/))') ' Lp of 1 1 1',myLp(1:3,:) - write (6,'(a,/,16(6(e12.4,x)/),2(f12.4,x))') 'state of 1 1 1',myState/1e6 - !$OMPEND CRITICAL (write2out) - endif - do k = 1,3 ! perturbation... - do l = 1,3 ! ...components - crystallite_subF(:,:,g,i,e) = myF ! initialize perturbed F to match converged - crystallite_subF(k,l,g,i,e) = crystallite_subF(k,l,g,i,e) + pert_Fg ! perturb single component - if (debugger) then - !$OMP CRITICAL (write2out) - write (6,*) '=============' - write (6,'(i1,x,i1)') k,l - write (6,*) '=============' - write (6,'(a,/,3(3(f12.6,x)/))') 'pertF of 1 1 1',crystallite_subF(1:3,:,g,i,e) - !$OMPEND CRITICAL (write2out) - endif - onTrack = .true. - converged = .false. - NiterationState = 0_pInt - do while(.not. converged .and. onTrack .and. NiterationState < nState) ! keep cycling until done (potentially non-converged) - NiterationState = NiterationState + 1_pInt - onTrack = crystallite_integrateStress(g,i,e) ! stress of perturbed situation (overwrites _P,_Tstar_v,_Fp,_Lp,_Fe) - if (onTrack) then - call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_Fp(:,:,g,i,e), & - crystallite_invFp(:,:,g,i,e), crystallite_Temperature(g,i,e), g, i, e) - - stateConverged = crystallite_updateState(g,i,e) ! update state - temperatureConverged = crystallite_updateTemperature(g,i,e) ! update temperature - converged = stateConverged .and. temperatureConverged - endif + if (crystallite_requested(g,i,e)) then ! first check whether is requested at all! + if (crystallite_converged(g,i,e)) then ! grain converged in above iteration + mySizeState = constitutive_sizeState(g,i,e) ! number of state variables for this grain + mySizeDotState = constitutive_sizeDotState(g,i,e) ! number of dotStates for this grain + myState(1:mySizeState) = constitutive_state(g,i,e)%p ! remember unperturbed, converged state, ... + myDotState(1:mySizeDotState) = constitutive_dotState(g,i,e)%p ! ... dotStates, ... + myTemperature = crystallite_Temperature(g,i,e) ! ... Temperature, ... + myF = crystallite_subF(:,:,g,i,e) ! ... and kinematics + myFp = crystallite_Fp(:,:,g,i,e) + myInvFp = crystallite_invFp(:,:,g,i,e) + myFe = crystallite_Fe(:,:,g,i,e) + myLp = crystallite_Lp(:,:,g,i,e) + myTstar_v = crystallite_Tstar_v(:,g,i,e) + myP = crystallite_P(:,:,g,i,e) + if (debugger) then + !$OMP CRITICAL (write2out) + write (6,*) '#############' + write (6,*) 'central solution of cryst_StressAndTangent' + write (6,*) '#############' + write (6,'(a,/,3(3(f12.4,x)/))') ' P of 1 1 1',myP(1:3,:)/1e6 + write (6,'(a,/,3(3(f12.8,x)/))') ' Fp of 1 1 1',myFp(1:3,:) + write (6,'(a,/,3(3(f12.8,x)/))') ' Lp of 1 1 1',myLp(1:3,:) + write (6,'(a,/,16(6(e12.4,x)/),2(f12.4,x))') 'state of 1 1 1',myState/1e6 + !$OMPEND CRITICAL (write2out) + endif + do k = 1,3 ! perturbation... + do l = 1,3 ! ...components + crystallite_subF(:,:,g,i,e) = myF ! initialize perturbed F to match converged + crystallite_subF(k,l,g,i,e) = crystallite_subF(k,l,g,i,e) + pert_Fg ! perturb single component if (debugger) then !$OMP CRITICAL (write2out) - write (6,*) '-------------' - write (6,'(l,x,l)') onTrack,converged - write (6,'(a,/,3(3(f12.4,x)/))') 'pertP of 1 1 1',crystallite_P(1:3,:,g,i,e)/1e6 - write (6,'(a,/,3(3(f12.4,x)/))') 'DP of 1 1 1',(crystallite_P(1:3,:,g,i,e)-myP(1:3,:))/1e6 - write (6,'(a,/,16(6(e12.4,x)/),/,2(f12.4,x))') 'state of 1 1 1',constitutive_state(g,i,e)%p/1e6 - write (6,'(a,/,16(6(e12.4,x)/),/,2(f12.4,x))') 'Dstate of 1 1 1',(constitutive_state(g,i,e)%p-myState)/1e6 + write (6,*) '=============' + write (6,'(i1,x,i1)') k,l + write (6,*) '=============' + write (6,'(a,/,3(3(f12.6,x)/))') 'pertF of 1 1 1',crystallite_subF(1:3,:,g,i,e) !$OMPEND CRITICAL (write2out) endif + onTrack = .true. + converged = .false. + NiterationState = 0_pInt + do while(.not. converged .and. onTrack .and. NiterationState < nState) ! keep cycling until done (potentially non-converged) + NiterationState = NiterationState + 1_pInt + onTrack = crystallite_integrateStress(g,i,e) ! stress of perturbed situation (overwrites _P,_Tstar_v,_Fp,_Lp,_Fe) + if (onTrack) then + call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_Fp(:,:,g,i,e), & + crystallite_invFp(:,:,g,i,e), crystallite_Temperature(g,i,e), g, i, e) + + stateConverged = crystallite_updateState(g,i,e) ! update state + temperatureConverged = crystallite_updateTemperature(g,i,e) ! update temperature + converged = stateConverged .and. temperatureConverged + endif + if (debugger) then + !$OMP CRITICAL (write2out) + write (6,*) '-------------' + write (6,'(a,x,l,x,l)') 'ontrack + converged:',onTrack,converged + write (6,'(a,/,3(3(f12.4,x)/))') 'pertP of 1 1 1',crystallite_P(1:3,:,g,i,e)/1e6 + write (6,'(a,/,3(3(f12.4,x)/))') 'DP of 1 1 1',(crystallite_P(1:3,:,g,i,e)-myP(1:3,:))/1e6 + write (6,'(a,/,16(6(e12.4,x)/),/,2(f12.4,x))') 'state of 1 1 1',constitutive_state(g,i,e)%p/1e6 + write (6,'(a,/,16(6(e12.4,x)/),/,2(f12.4,x))') 'Dstate of 1 1 1',(constitutive_state(g,i,e)%p-myState)/1e6 + !$OMPEND CRITICAL (write2out) + endif + enddo + if (converged) & ! converged state warrants stiffness update + crystallite_dPdF(:,:,k,l,g,i,e) = (crystallite_P(:,:,g,i,e) - myP)/pert_Fg ! tangent dP_ij/dFg_kl + constitutive_state(g,i,e)%p = myState ! restore unperturbed, converged state, ... + constitutive_dotState(g,i,e)%p = myDotState ! ... dotState, ... + crystallite_Temperature(g,i,e) = myTemperature ! ... temperature, ... + crystallite_Fp(:,:,g,i,e) = myFp ! ... and kinematics + crystallite_invFp(:,:,g,i,e) = myInvFp + crystallite_Fe(:,:,g,i,e) = myFe + crystallite_Lp(:,:,g,i,e) = myLp + crystallite_Tstar_v(:,g,i,e) = myTstar_v + crystallite_P(:,:,g,i,e) = myP + !$OMP CRITICAL (out) + debug_StiffnessStateLoopDistribution(NiterationState) = & + debug_StiffnessstateLoopDistribution(NiterationState) + 1 + !$OMPEND CRITICAL (out) enddo - if (converged) & ! converged state warrants stiffness update - crystallite_dPdF(:,:,k,l,g,i,e) = (crystallite_P(:,:,g,i,e) - myP)/pert_Fg ! tangent dP_ij/dFg_kl - constitutive_state(g,i,e)%p = myState ! restore unperturbed, converged state, ... - constitutive_dotState(g,i,e)%p = myDotState ! ... dotState, ... - crystallite_Temperature(g,i,e) = myTemperature ! ... temperature, ... - crystallite_Fp(:,:,g,i,e) = myFp ! ... and kinematics - crystallite_invFp(:,:,g,i,e) = myInvFp - crystallite_Fe(:,:,g,i,e) = myFe - crystallite_Lp(:,:,g,i,e) = myLp - crystallite_Tstar_v(:,g,i,e) = myTstar_v - crystallite_P(:,:,g,i,e) = myP - !$OMP CRITICAL (out) - debug_StiffnessStateLoopDistribution(NiterationState) = & - debug_StiffnessstateLoopDistribution(NiterationState) + 1 - !$OMPEND CRITICAL (out) enddo - enddo - - else ! grain did not converged - crystallite_dPdF(:,:,:,:,g,i,e) = crystallite_fallbackdPdF(:,:,:,:,g,i,e) ! use fallback - endif - enddo - enddo - enddo + else ! grain did not converged + crystallite_dPdF(:,:,:,:,g,i,e) = crystallite_fallbackdPdF(:,:,:,:,g,i,e) ! use (elastic) fallback + endif ! grain convergence + endif ! grain request + enddo ! grain loop + enddo ! ip loop + enddo ! element loop !$OMPEND PARALLEL DO - endif + endif ! jacobian calculation endsubroutine diff --git a/code/homogenization.f90 b/code/homogenization.f90 index 1ccfb98e3..dca439242 100644 --- a/code/homogenization.f90 +++ b/code/homogenization.f90 @@ -266,28 +266,28 @@ subroutine materialpoint_stressAndItsTangent(& do while (any(materialpoint_subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMin)) ! cutback loop for material points +! write(6,'(a,/,125(8(f8.5,x),/))') 'mp_subSteps',materialpoint_subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed - ! debugger = (e == 1 .and. i == 1) + debugger = (e == 1 .and. i == 1) ! if our materialpoint converged then we are either finished or have to wind forward if (materialpoint_converged(i,e)) then + if (debugger) then + !$OMP CRITICAL (write2out) + write(6,'(a21,f10.8,a34,f10.8,a37,/)') 'winding forward from ', & + materialpoint_subFrac(i,e), ' to current materialpoint_subFrac ', & + materialpoint_subFrac(i,e)+materialpoint_subStep(i,e),' in materialpoint_stressAndItsTangent' + !$OMPEND CRITICAL (write2out) + endif ! calculate new subStep and new subFrac materialpoint_subFrac(i,e) = materialpoint_subFrac(i,e) + materialpoint_subStep(i,e) materialpoint_subStep(i,e) = min(1.0_pReal-materialpoint_subFrac(i,e), 2.0_pReal * materialpoint_subStep(i,e)) - - if (debugger) then - !$OMP CRITICAL (write2out) - write(6,'(a21,f10.8,a34,f10.8,a37,/)') 'winding forward from ', & - materialpoint_subFrac(i,e) - materialpoint_subStep(i,e),' to current materialpoint_subFrac ', & - materialpoint_subFrac(i,e),' in materialpoint_stressAndItsTangent' - !$OMPEND CRITICAL (write2out) - endif - + ! still stepping needed if (materialpoint_subStep(i,e) > subStepMin) then @@ -357,6 +357,9 @@ subroutine materialpoint_stressAndItsTangent(& ) .and. NiterationMPstate < nMPstate) ! convergence loop for materialpoint NiterationMPstate = NiterationMPstate + 1 +! write(6,'(a,/,125(8(l,x),/))') 'material point request and not done', & +! materialpoint_requested .and. .not. materialpoint_doneAndHappy(1,:,:) + ! --+>> deformation partitioning <<+-- ! ! based on materialpoint_subF0,.._subF, @@ -373,19 +376,22 @@ subroutine materialpoint_stressAndItsTangent(& call homogenization_partitionDeformation(i,e) ! partition deformation onto constituents crystallite_dt(1:myNgrains,i,e) = materialpoint_subdt(i,e) ! propagate materialpoint dt to grains crystallite_requested(1:myNgrains,i,e) = .true. ! request calculation for constituents + else + crystallite_requested(1:myNgrains,i,e) = .false. ! calculation for constituents not required anymore endif enddo enddo !$OMP END PARALLEL DO - +! write(6,'(a,/,125(8(8(l,x),2x),/))') 'crystallite request with updated partitioning', crystallite_requested ! --+>> crystallite integration <<+-- ! ! based on crystallite_partionedF0,.._partionedF ! incrementing by crystallite_dt - call crystallite_stressAndItsTangent(updateJaco) ! request stress and tangent calculation for constituent grains + +! write(6,'(a,/,125(8(8(l,x),2x),/))') 'crystallite converged', crystallite_converged ! --+>> state update <<+-- @@ -407,6 +413,8 @@ subroutine materialpoint_stressAndItsTangent(& enddo enddo !$OMP END PARALLEL DO +! write(6,'(a,/,125(8(l,x),/))') 'material point done', materialpoint_doneAndHappy(1,:,:) +! write(6,'(a,/,125(8(l,x),/))') 'material point converged', materialpoint_converged enddo ! homogenization convergence loop @@ -417,7 +425,6 @@ subroutine materialpoint_stressAndItsTangent(& ! check for non-performer: any(.not. converged) ! replace everybody with odd response ? - !$OMP PARALLEL DO elementLoop: do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed @@ -432,9 +439,9 @@ elementLoop: do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate enddo elementLoop !$OMP END PARALLEL DO - write (6,*) 'Material Point finished' - - ! how to deal with stiffness? + write (6,*) + write (6,*) 'Material Point end' + write (6,*) return endsubroutine diff --git a/code/homogenization_RGC.f90 b/code/homogenization_RGC.f90 index 475842456..081a40795 100644 --- a/code/homogenization_RGC.f90 +++ b/code/homogenization_RGC.f90 @@ -25,7 +25,6 @@ MODULE homogenization_RGC homogenization_RGC_ciAlpha character(len=64), dimension(:,:), allocatable :: homogenization_RGC_output - CONTAINS !**************************************** !* - homogenization_RGC_init @@ -147,8 +146,8 @@ function homogenization_RGC_stateInit(myInstance) integer(pInt), intent(in) :: myInstance real(pReal), dimension(homogenization_RGC_sizeState(myInstance)) :: homogenization_RGC_stateInit -!* Open a debugging file << not used at the moment >> -! open(1978,file='homogenization_RGC_debugging.out') +!* Open a debugging file + open(1978,file='homogenization_RGC_debugging.out') homogenization_RGC_stateInit = 0.0_pReal return @@ -184,42 +183,42 @@ subroutine homogenization_RGC_partitionDeformation(& integer(pInt), dimension (4) :: intFace integer(pInt), dimension (3) :: iGrain3 integer(pInt) homID, iGrain,iFace,i,j + logical RGCdebug ! integer(pInt), parameter :: nFace = 6 - homID = homogenization_typeInstance(mesh_element(3,el)) - F = 0.0_pReal + RGCdebug = (el == 1 .and. ip == 1) + !* Debugging the overall deformation gradient -! if (ip == 1 .and. el == 1) then -! write(1978,'(x,a32)')'Overall deformation gradient: ' -! do i = 1,3 -! write(1978,'(x,3(e10.4,x))')(avgF(i,j), j = 1,3) -! enddo -! endif + if (RGCdebug) then + write(1978,'(x,a32,i3,i3)')'Overall deformation gradient: ' + do i = 1,3 + write(1978,'(x,3(e10.4,x))')(avgF(i,j), j = 1,3) + enddo + write(1978,*)' ' + endif !* Compute the deformation gradient of individual grains due to relaxations + homID = homogenization_typeInstance(mesh_element(3,el)) + F = 0.0_pReal do iGrain = 1,homogenization_Ngrains(mesh_element(3,el)) call homogenization_RGC_grain1to3(iGrain3,iGrain,homID) do iFace = 1,nFace call homogenization_RGC_getInterface(intFace,iFace,iGrain3) call homogenization_RGC_relaxationVector(aVect,intFace,state,homID) -!* Debugging the grain relaxation vectors -! if (ip == 1 .and. el == 1) then -! write(1978,'(x,a32,x,i3)')'Relaxation vector of interface: ',iFace -! write(1978,'(x,3(e10.4,x))')(aVect(j), j = 1,3) -! endif call homogenization_RGC_interfaceNormal(nVect,intFace) forall (i=1:3,j=1:3) & F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations enddo F(:,:,iGrain) = F(:,:,iGrain) + avgF(:,:) ! relaxed deformation gradient !* Debugging the grain deformation gradients -! if (ip == 1 .and. el == 1) then -! write(1978,'(x,a32,x,i3)')'Deformation gradient of grain: ',iGrain -! do i = 1,3 -! write(1978,'(x,3(e10.4,x))')(F(i,j,iGrain), j = 1,3) -! enddo -! endif + if (RGCdebug) then + write(1978,'(x,a32,x,i3)')'Deformation gradient of grain: ',iGrain + do i = 1,3 + write(1978,'(x,3(e10.4,x))')(F(i,j,iGrain), j = 1,3) + enddo + write(1978,*)' ' + endif enddo return @@ -265,26 +264,49 @@ function homogenization_RGC_updateState(& real(pReal), dimension (homogenization_maxNgrains) :: NN,pNN real(pReal), dimension (3) :: normP,normN,mornP,mornN real(pReal) residMax,stresMax,constitutiveWork,penaltyEnergy - logical error + logical error,RGCdebug,RGCdebugJacobi ! integer(pInt), parameter :: nFace = 6 ! real(pReal), dimension(:,:), allocatable :: tract,jmatrix,jnverse,smatrix,pmatrix real(pReal), dimension(:), allocatable :: resid,relax,p_relax,p_resid + RGCdebug = (el == 1 .and. ip == 1) + RGCdebugJacobi = .false. + !* Get the dimension of the cluster (grains and interfaces) homID = homogenization_typeInstance(mesh_element(3,el)) nGDim = homogenization_RGC_Ngrains(:,homID) nIntFaceTot = (nGDim(1)-1)*nGDim(2)*nGDim(3) + nGDim(1)*(nGDim(2)-1)*nGDim(3) & + nGDim(1)*nGDim(2)*(nGDim(3)-1) -!* Allocate the size of the arrays/matrices depending on the size of the cluster +!* Allocate the size of the arrays/matrices depending on the size of the cluster allocate(resid(3*nIntFaceTot)); resid = 0.0_pReal allocate(tract(nIntFaceTot,3)); tract = 0.0_pReal allocate(relax(3*nIntFaceTot)); relax = state%p(1:3*nIntFaceTot) +!* Debugging the obtained state + if (RGCdebug) then + write(1978,'(x,a30)')'Obtained state: ' + do i = 1,3*nIntFaceTot + write(1978,'(x,2(e10.4,x))')state%p(i) + enddo + write(1978,*)' ' + endif !* Stress-like penalty related to mismatch or incompatibility at interfaces call homogenization_RGC_stressPenalty(R,NN,F,ip,el,homID) +!* Debugging the mismatch, stress and penalty of grains + if (RGCdebug) then + do iGrain = 1,homogenization_Ngrains(mesh_element(3,el)) + write(1978,'(x,a30,x,i3,x,a4,x,e10.4)')'Mismatch magnitude of grain(',iGrain,') :',NN(iGrain) + write(1978,*)' ' + write(1978,'(x,a30,x,i3)')'Stress and penalty of grain: ',iGrain + do i = 1,3 + write(1978,'(x,3(e10.4,x),x,3(e10.4,x))')(P(i,j,iGrain), j = 1,3),(R(i,j,iGrain), j = 1,3) + enddo + write(1978,*)' ' + enddo + endif !* Compute the residual stress from the balance of traction at all (interior) interfaces do iNum = 1,nIntFaceTot @@ -300,13 +322,6 @@ function homogenization_RGC_updateState(& call homogenization_RGC_grain3to1(iGrP,iGr3P,homID) call homogenization_RGC_getInterface(intFaceP,2*faceID(1)-1,iGr3P) call homogenization_RGC_interfaceNormal(normP,intFaceP) ! get the interface normal -!* Debugging the grains and their stresses -! if (ip == 1 .and. el == 1) then -! write(1978,'(x,a30,2(x,i3))')'Stresses of grains: ',iGrN(iNum),iGrP(iNum) -! do i = 1,3 -! write(1978,'(x,3(e10.4,x),x,3(e10.4,x))')(P(i,j,iGrN(iNum)), j = 1,3),(P(i,j,iGrP(iNum)), j = 1,3) -! enddo -! endif do i = 1,3 ! compute the traction balance at the interface do j = 1,3 tract(iNum,i) = tract(iNum,i) + (P(i,j,iGrP) + R(i,j,iGrP))*normP(j) & @@ -315,10 +330,11 @@ function homogenization_RGC_updateState(& enddo enddo !* Debugging the residual stress -! if (ip == 1 .and. el == 1) then -! write(1978,'(x,a30,x,i3)')'Traction difference: ',iNum -! write(1978,'(x,3(e10.4,x))')(tract(iNum,j), j = 1,3) -! endif + if (RGCdebug) then + write(1978,'(x,a30,x,i3)')'Traction difference: ',iNum + write(1978,'(x,3(e10.4,x))')(tract(iNum,j), j = 1,3) + write(1978,*)' ' + endif enddo !* Convergence check for stress residual @@ -326,28 +342,27 @@ function homogenization_RGC_updateState(& stresLoc = maxloc(P) residMax = maxval(tract) residLoc = maxloc(tract) -!* Temporary debugging statement << not used at the moment >> -! if (ip == 1 .and. el == 1) then -! write(1978,'(x,a)')' ' -! write(1978,'(x,a)')'Residual check ...' -! write(1978,'(x,a15,x,e10.4,x,a7,i3,x,a12,i2,i2)')'Max stress: ',stresMax, & -! '@ grain',stresLoc(3),'in component',stresLoc(1),stresLoc(2) -! write(1978,'(x,a15,x,e10.4,x,a7,i3,x,a12,i2)')'Max residual: ',residMax, & -! '@ iface',residLoc(1),'in direction',residLoc(2) -! endif +!* Debugging the convergent criteria + if (RGCdebug) then + write(1978,'(x,a)')' ' + write(1978,'(x,a)')'Residual check ...' + write(1978,'(x,a15,x,e10.4,x,a7,i3,x,a12,i2,i2)')'Max stress: ',stresMax, & + '@ grain',stresLoc(3),'in component',stresLoc(1),stresLoc(2) + write(1978,'(x,a15,x,e10.4,x,a7,i3,x,a12,i2)')'Max residual: ',residMax, & + '@ iface',residLoc(1),'in direction',residLoc(2) + endif homogenization_RGC_updateState = .false. -!* If convergence reached => done and happy +!*** If convergence reached => done and happy if (residMax < relTol_RGC*stresMax .or. residMax < absTol_RGC) then homogenization_RGC_updateState = .true. -!* Temporary debugging statement << not used at the moment >> -! if (ip == 1 .and. el == 1) then -! write(1978,'(x,a55)')'... done and happy' -! endif -!* Compute/update the state for postResult, i.e., ... -!* ... the (bulk) constitutive work, + if (RGCdebug) then + write(1978,'(x,a55)')'... done and happy' + write(1978,*)' ' + endif + write(6,'(x,a,x,i3,x,a6,x,i3,x,a12)')'RGC_updateState: ip',ip,'| el',el,'converged :)' +!* Then compute/update the state for postResult, i.e., ... +!* ... the (bulk) constitutive work and penalty energy constitutiveWork = state%p(3*nIntFaceTot+1) - state%p(3*nIntFaceTot+1) = constitutiveWork -!* ... the penalty energy, and penaltyEnergy = state%p(3*nIntFaceTot+2) do iGrain = 1,homogenization_Ngrains(mesh_element(3,el)) do i = 1,3 @@ -357,26 +372,34 @@ function homogenization_RGC_updateState(& enddo enddo enddo + state%p(3*nIntFaceTot+1) = constitutiveWork state%p(3*nIntFaceTot+2) = penaltyEnergy !* ... the overall mismatch state%p(3*nIntFaceTot+3) = sum(NN) + if (RGCdebug) then + write(1978,'(x,a30,x,e10.4)')'Constitutive work: ',constitutiveWork + write(1978,'(x,a30,x,e10.4)')'Penalty energy: ',penaltyEnergy + write(1978,'(x,a30,x,e10.4)')'Magnitude mismatch: ',sum(NN) + write(1978,*)' ' + endif deallocate(tract,resid,relax) return -!* If residual blows-up => done but unhappy +!*** If residual blows-up => done but unhappy elseif (residMax > relMax_RGC*stresMax .or. residMax > absMax_RGC) then - homogenization_RGC_updateState(1) = .true. -!* Temporary debugging statement << not used at the moment >> -! if (ip == 1 .and. el == 1) then -! write(1978,'(x,a55)')'... done but not happy' -! endif + homogenization_RGC_updateState = (/.true.,.false./) + if (RGCdebug) then + write(1978,'(x,a55)')'... broken' + write(1978,*)' ' + endif + write(6,'(x,a,x,i3,x,a6,x,i3,x,a9)')'RGC_updateState: ip',ip,'| el',el,'broken :(' deallocate(tract,resid,relax) return -!* Otherwise, proceed with computing the state update +!*** Otherwise, proceed with computing the Jacobian and state update else -!* Temporary debugging statement << not used at the moment >> -! if (ip == 1 .and. el == 1) then -! write(1978,'(x,a55)')'... not done' -! endif + if (RGCdebug) then + write(1978,'(x,a55)')'... not yet done' + write(1978,*)' ' + endif endif !* Construct the Jacobian matrix of the constitutive stress tangent from dPdF @@ -414,12 +437,13 @@ function homogenization_RGC_updateState(& enddo enddo !* Debugging the global Jacobian matrix of stress tangent -! if (ip == 1 .and. el == 1) then -! write(1978,'(x,a24)')'Jacobian matrix of stress' -! do i = 1,3*nIntFaceTot -! write(1978,'(x,400(e10.4,x))')(smatrix(i,j), j = 1,3*nIntFaceTot) -! enddo -! endif + if (RGCdebugJacobi) then + write(1978,'(x,a30)')'Jacobian matrix of stress' + do i = 1,3*nIntFaceTot + write(1978,'(x,100(e10.4,x))')(smatrix(i,j), j = 1,3*nIntFaceTot) + enddo + write(1978,*)' ' + endif !* Compute the Jacobian of the stress-like penalty (penalty tangent) using perturbation technique allocate(pmatrix(3*nIntFaceTot,3*nIntFaceTot)); pmatrix = 0.0_pReal @@ -429,14 +453,7 @@ function homogenization_RGC_updateState(& p_relax = relax p_relax(ipert) = relax(ipert) + pPert_RGC ! perturb the relaxation vector state%p(1:3*nIntFaceTot) = p_relax -!* Debugging the perturbed state -! if (ip == 1 .and. el == 1) then -! write(1978,'(x,a32)')'State and perturbed state: ' -! do i = 1,3*nIntFaceTot -! write(1978,'(x,2(e10.4,x))')relax(i),pelax(i) -! enddo -! endif - call homogenization_RGC_partitionDeformation(pF,F0,avgF,state,ip,el) + call homogenization_RGC_grainDeformation(pF,F0,avgF,state,el) call homogenization_RGC_stressPenalty(pR,pNN,pF,ip,el,homID) p_resid = 0.0_pReal do iNum = 1,nIntFaceTot @@ -463,35 +480,45 @@ function homogenization_RGC_updateState(& pmatrix(:,ipert) = p_resid/pPert_RGC enddo !* Debugging the global Jacobian matrix of penalty tangent -! if (ip == 1 .and. el == 1) then -! write(1978,'(x,a24)')'Jacobian matrix of penalty' -! do i = 1,3*nIntFaceTot -! write(1978,'(x,400(e10.4,x))')(pmatrix(i,j), j = 1,3*nIntFaceTot) -! enddo -! endif + if (RGCdebugJacobi) then + write(1978,'(x,a30)')'Jacobian matrix of penalty' + do i = 1,3*nIntFaceTot + write(1978,'(x,100(e10.4,x))')(pmatrix(i,j), j = 1,3*nIntFaceTot) + enddo + write(1978,*)' ' + endif !* The overall Jacobian matrix (due to constitutive and penalty tangents) allocate(jmatrix(3*nIntFaceTot,3*nIntFaceTot)); jmatrix = smatrix + pmatrix + if (RGCdebugJacobi) then + write(1978,'(x,a30)')'Jacobian matrix (total)' + do i = 1,3*nIntFaceTot + write(1978,'(x,100(e10.4,x))')(jmatrix(i,j), j = 1,3*nIntFaceTot) + enddo + write(1978,*)' ' + endif allocate(jnverse(3*nIntFaceTot,3*nIntFaceTot)); jnverse = 0.0_pReal call math_invert(3*nIntFaceTot,jmatrix,jnverse,ival,error) !* Debugging the inverse Jacobian matrix -! if (ip == 1 .and. el == 1) then -! write(1978,'(x,a20)')'Jacobian inverse' -! do i = 1,3*nIntFaceTot -! write(1978,'(x,400(e10.4,x))')(jnverse(i,j), j = 1,3*nIntFaceTot) -! enddo -! endif + if (RGCdebugJacobi) then + write(1978,'(x,a30)')'Jacobian inverse' + do i = 1,3*nIntFaceTot + write(1978,'(x,100(e10.4,x))')(jnverse(i,j), j = 1,3*nIntFaceTot) + enddo + write(1978,*)' ' + endif !* Calculate the state update (i.e., new relaxation vectors) forall(i=1:3*nIntFaceTot,j=1:3*nIntFaceTot) relax(i) = relax(i) - jnverse(i,j)*resid(j) state%p(1:3*nIntFaceTot) = relax !* Debugging the return state -! if (ip == 1 .and. el == 1) then -! write(1978,'(x,a32)')'Returned state: ' -! do i = 1,3*nIntFaceTot -! write(1978,'(x,2(e10.4,x))')state%p(i) -! enddo -! endif + if (RGCdebugJacobi) then + write(1978,'(x,a30)')'Returned state: ' + do i = 1,3*nIntFaceTot + write(1978,'(x,2(e10.4,x))')state%p(i) + enddo + write(1978,*)' ' + endif deallocate(tract,resid,jmatrix,jnverse,relax,pmatrix,smatrix,p_relax,p_resid) return @@ -514,6 +541,7 @@ subroutine homogenization_RGC_averageStressAndItsTangent(& use prec, only: pReal,pInt,p_vec use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips use material, only: homogenization_maxNgrains, homogenization_Ngrains,homogenization_typeInstance + use math, only: math_Plain3333to99 implicit none !* Definition of variables @@ -521,12 +549,26 @@ subroutine homogenization_RGC_averageStressAndItsTangent(& real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: P real(pReal), dimension (3,3,3,3,homogenization_maxNgrains), intent(in) :: dPdF + real(pReal), dimension (9,9) :: dPdF99 integer(pInt), intent(in) :: ip,el ! - logical homogenization_RGC_stateUpdate - integer(pInt) homID, i, Ngrains + logical homogenization_RGC_stateUpdate,RGCdebug + integer(pInt) homID, i, j, Ngrains, iGrain -! homID = homogenization_typeInstance(mesh_element(3,el)) ! <> + RGCdebug = .false. !(ip == 1 .and. el == 1) + + homID = homogenization_typeInstance(mesh_element(3,el)) +!* Debugging the grain tangent + if (RGCdebug) then + do iGrain = 1,homogenization_Ngrains(mesh_element(3,el)) + dPdF99 = math_Plain3333to99(dPdF(:,:,:,:,iGrain)) + write(1978,'(x,a30,x,i3)')'Stress tangent of grain: ',iGrain + do i = 1,9 + write(1978,'(x,(e10.4,x))') (dPdF99(i,j), j = 1,9) + enddo + write(1978,*)' ' + enddo + endif Ngrains = homogenization_Ngrains(mesh_element(3,el)) avgP = sum(P,3)/dble(Ngrains) dAvgPdAvgF = sum(dPdF,5)/dble(Ngrains) @@ -993,4 +1035,53 @@ subroutine homogenization_RGC_interface1to4(& endsubroutine +!******************************************************************** +! calculating the grain deformation gradient +!******************************************************************** +subroutine homogenization_RGC_grainDeformation(& + F, & ! partioned def grad per grain +! + F0, & ! initial partioned def grad per grain + avgF, & ! my average def grad + state, & ! my state + el & ! my element + ) + use prec, only: pReal,pInt,p_vec + use mesh, only: mesh_element + use material, only: homogenization_maxNgrains,homogenization_Ngrains,homogenization_typeInstance + implicit none + +!* Definition of variables + real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F + real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: F0 + real(pReal), dimension (3,3), intent(in) :: avgF + type(p_vec), intent(in) :: state + integer(pInt), intent(in) :: el +! + real(pReal), dimension (3) :: aVect,nVect + integer(pInt), dimension (4) :: intFace + integer(pInt), dimension (3) :: iGrain3 + integer(pInt) homID, iGrain,iFace,i,j +! + integer(pInt), parameter :: nFace = 6 + +!* Compute the deformation gradient of individual grains due to relaxations + homID = homogenization_typeInstance(mesh_element(3,el)) + F = 0.0_pReal + do iGrain = 1,homogenization_Ngrains(mesh_element(3,el)) + call homogenization_RGC_grain1to3(iGrain3,iGrain,homID) + do iFace = 1,nFace + call homogenization_RGC_getInterface(intFace,iFace,iGrain3) + call homogenization_RGC_relaxationVector(aVect,intFace,state,homID) + call homogenization_RGC_interfaceNormal(nVect,intFace) + forall (i=1:3,j=1:3) & + F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations + enddo + F(:,:,iGrain) = F(:,:,iGrain) + avgF(:,:) ! relaxed deformation gradient + enddo + + return + +endsubroutine + END MODULE