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

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.
This commit is contained in:
Denny Tjahjanto 2009-08-27 12:10:06 +00:00
parent f70013ee9b
commit b25396374a
4 changed files with 325 additions and 215 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
code/constitutive_phenopowerlaw.f90
View File

@ -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

184
code/crystallite.f90
View File

@ -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

39
code/homogenization.f90
View File

@ -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

289
code/homogenization_RGC.f90
View File

@ -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)) ! <<not required at the moment>>
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