removed numerical perturbation calculation

This commit is contained in:
Martin Diehl 2016-07-25 20:07:12 +02:00
parent 8235bf3422
commit e4c590699f
1 changed files with 97 additions and 371 deletions

View File

@ -518,8 +518,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
nCryst, & nCryst, &
numerics_integrator, & numerics_integrator, &
numerics_integrationMode, & numerics_integrationMode, &
numerics_timeSyncing, & numerics_timeSyncing
analyticJaco
use debug, only: & use debug, only: &
debug_level, & debug_level, &
debug_crystallite, & debug_crystallite, &
@ -582,23 +581,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
invFp, & ! inverse of the plastic deformation gradient invFp, & ! inverse of the plastic deformation gradient
Fe_guess, & ! guess for elastic deformation gradient Fe_guess, & ! guess for elastic deformation gradient
Tstar ! 2nd Piola-Kirchhoff stress tensor Tstar ! 2nd Piola-Kirchhoff stress tensor
real(pReal), allocatable, dimension(:,:,:,:,:,:,:) :: &
dPdF_perturbation1, &
dPdF_perturbation2
real(pReal), allocatable, dimension(:,:,:,:,:) :: &
F_backup, &
Fp_backup, &
InvFp_backup, &
Fi_backup, &
InvFi_backup, &
Fe_backup, &
Lp_backup, &
Li_backup, &
P_backup
real(pReal), allocatable, dimension(:,:,:,:) :: &
Tstar_v_backup
logical, allocatable, dimension(:,:,:) :: &
convergenceFlag_backup
integer(pInt) :: & integer(pInt) :: &
NiterationCrystallite, & ! number of iterations in crystallite loop NiterationCrystallite, & ! number of iterations in crystallite loop
c, & !< counter in integration point component loop c, & !< counter in integration point component loop
@ -1137,10 +1119,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
! --+>> STIFFNESS CALCULATION <<+-- ! --+>> STIFFNESS CALCULATION <<+--
computeJacobian: if(updateJaco) then computeJacobian: if(updateJaco) then
jacobianMethod: if (analyticJaco) then
! --- ANALYTIC JACOBIAN ---
!$OMP PARALLEL DO PRIVATE(dSdF,dSdFe,dSdFi,dLpdS,dLpdFi,dFpinvdF,dLidS,dLidFi,dFidS,& !$OMP PARALLEL DO PRIVATE(dSdF,dSdFe,dSdFi,dLpdS,dLpdFi,dFpinvdF,dLidS,dLidFi,dFidS,&
!$OMP rhs_3333,lhs_3333,temp_99,temp_33,temp_3333,myNcomponents,error) !$OMP rhs_3333,lhs_3333,temp_99,temp_33,temp_3333,myNcomponents,error)
elementLooping6: do e = FEsolving_execElem(1),FEsolving_execElem(2) elementLooping6: do e = FEsolving_execElem(1),FEsolving_execElem(2)
@ -1250,258 +1228,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
enddo; enddo enddo; enddo
enddo elementLooping6 enddo elementLooping6
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
else jacobianMethod
! --- STANDARD (PERTURBATION METHOD) FOR JACOBIAN ---
numerics_integrationMode = 2_pInt
! --- BACKUP ---
allocate(dPdF_perturbation1(3,3,3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal)
allocate(dPdF_perturbation2(3,3,3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal)
allocate(F_backup (3,3, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal)
allocate(Fp_backup (3,3, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal)
allocate(InvFp_backup (3,3, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal)
allocate(Fi_backup (3,3, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal)
allocate(InvFi_backup (3,3, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal)
allocate(Fe_backup (3,3, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal)
allocate(Lp_backup (3,3, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal)
allocate(Li_backup (3,3, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal)
allocate(P_backup (3,3, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal)
allocate(Tstar_v_backup (6, homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = 0.0_pReal)
allocate(convergenceFlag_backup (homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source = .false.)
!$OMP PARALLEL DO PRIVATE(myNcomponents)
elementLooping7: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNcomponents = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e); do c = 1,myNcomponents
plasticState (phaseAt(c,i,e))%state_backup(:,phasememberAt(c,i,e)) = &
plasticState (phaseAt(c,i,e))%state( :,phasememberAt(c,i,e))
do mySource = 1_pInt, phase_Nsources(phaseAt(c,i,e))
sourceState(phaseAt(c,i,e))%p(mySource)%state_backup(:,phasememberAt(c,i,e)) = &
sourceState(phaseAt(c,i,e))%p(mySource)%state( :,phasememberAt(c,i,e))
enddo
plasticState (phaseAt(c,i,e))%dotState_backup(:,phasememberAt(c,i,e)) = &
plasticState (phaseAt(c,i,e))%dotState( :,phasememberAt(c,i,e))
do mySource = 1_pInt, phase_Nsources(phaseAt(c,i,e))
sourceState(phaseAt(c,i,e))%p(mySource)%dotState_backup(:,phasememberAt(c,i,e)) = &
sourceState(phaseAt(c,i,e))%p(mySource)%dotState( :,phasememberAt(c,i,e))
enddo
F_backup(1:3,1:3,c,i,e) = crystallite_subF(1:3,1:3,c,i,e) ! ... and kinematics
Fp_backup(1:3,1:3,c,i,e) = crystallite_Fp(1:3,1:3,c,i,e)
InvFp_backup(1:3,1:3,c,i,e) = crystallite_invFp(1:3,1:3,c,i,e)
Fi_backup(1:3,1:3,c,i,e) = crystallite_Fi(1:3,1:3,c,i,e)
InvFi_backup(1:3,1:3,c,i,e) = crystallite_invFi(1:3,1:3,c,i,e)
Fe_backup(1:3,1:3,c,i,e) = crystallite_Fe(1:3,1:3,c,i,e)
Lp_backup(1:3,1:3,c,i,e) = crystallite_Lp(1:3,1:3,c,i,e)
Li_backup(1:3,1:3,c,i,e) = crystallite_Li(1:3,1:3,c,i,e)
Tstar_v_backup(1:6,c,i,e) = crystallite_Tstar_v(1:6,c,i,e)
P_backup(1:3,1:3,c,i,e) = crystallite_P(1:3,1:3,c,i,e)
convergenceFlag_backup(c,i,e) = crystallite_converged(c,i,e)
enddo; enddo
enddo elementLooping7
!$END PARALLEL DO
! --- CALCULATE STATE AND STRESS FOR PERTURBATION ---
dPdF_perturbation1 = crystallite_dPdF0 ! initialize stiffness with known good values from last increment
dPdF_perturbation2 = crystallite_dPdF0 ! initialize stiffness with known good values from last increment
pertubationLoop: do perturbation = 1,2 ! forward and backward perturbation
if (iand(pert_method,perturbation) > 0_pInt) then ! mask for desired direction
myPert = -pert_Fg * (-1.0_pReal)**perturbation ! set perturbation step
do k = 1,3; do l = 1,3 ! ...alter individual components
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. c == debug_g) &
.or. .not. iand(debug_level(debug_crystallite),debug_levelSelective) /= 0_pInt)) &
write(6,'(a,2(1x,i1),1x,a,/)') '<< CRYST >> [[[[[[ Stiffness perturbation',k,l,']]]]]]'
! --- INITIALIZE UNPERTURBED STATE ---
select case(numerics_integrator(numerics_integrationMode))
case(1_pInt)
!why not OMP? ! Fix-point method: restore to last converged state at end of subinc, since this is probably closest to perturbed state
do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNcomponents = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e); do c = 1,myNcomponents
plasticState (phaseAt(c,i,e))%state( :,phasememberAt(c,i,e)) = &
plasticState (phaseAt(c,i,e))%state_backup(:,phasememberAt(c,i,e))
do mySource = 1_pInt, phase_Nsources(phaseAt(c,i,e))
sourceState(phaseAt(c,i,e))%p(mySource)%state( :,phasememberAt(c,i,e)) = &
sourceState(phaseAt(c,i,e))%p(mySource)%state_backup(:,phasememberAt(c,i,e))
enddo
plasticState (phaseAt(c,i,e))%dotState( :,phasememberAt(c,i,e)) = &
plasticState (phaseAt(c,i,e))%dotState_backup(:,phasememberAt(c,i,e))
do mySource = 1_pInt, phase_Nsources(phaseAt(c,i,e))
sourceState(phaseAt(c,i,e))%p(mySource)%dotState( :,phasememberAt(c,i,e)) = &
sourceState(phaseAt(c,i,e))%p(mySource)%dotState_backup(:,phasememberAt(c,i,e))
enddo
crystallite_Fp(1:3,1:3,c,i,e) = Fp_backup(1:3,1:3,c,i,e)
crystallite_invFp(1:3,1:3,c,i,e) = InvFp_backup(1:3,1:3,c,i,e)
crystallite_Fi(1:3,1:3,c,i,e) = Fi_backup(1:3,1:3,c,i,e)
crystallite_invFi(1:3,1:3,c,i,e) = InvFi_backup(1:3,1:3,c,i,e)
crystallite_Fe(1:3,1:3,c,i,e) = Fe_backup(1:3,1:3,c,i,e)
crystallite_Lp(1:3,1:3,c,i,e) = Lp_backup(1:3,1:3,c,i,e)
crystallite_Li(1:3,1:3,c,i,e) = Li_backup(1:3,1:3,c,i,e)
crystallite_Tstar_v(1:6,c,i,e) = Tstar_v_backup(1:6,c,i,e)
enddo; enddo
enddo
case(2_pInt,3_pInt) ! explicit Euler methods: nothing to restore (except for F), since we are only doing a stress integration step
case(4_pInt,5_pInt)
!why not OMP? ! explicit Runge-Kutta methods: restore to start of subinc, since we are doing a full integration of state and stress
do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNcomponents = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e); do c = 1,myNcomponents
plasticState (phaseAt(c,i,e))%state( :,phasememberAt(c,i,e)) = &
plasticState (phaseAt(c,i,e))%subState0(:,phasememberAt(c,i,e))
do mySource = 1_pInt, phase_Nsources(phaseAt(c,i,e))
sourceState(phaseAt(c,i,e))%p(mySource)%state( :,phasememberAt(c,i,e)) = &
sourceState(phaseAt(c,i,e))%p(mySource)%subState0(:,phasememberAt(c,i,e))
enddo
plasticState (phaseAt(c,i,e))%dotState( :,phasememberAt(c,i,e)) = &
plasticState (phaseAt(c,i,e))%dotState_backup(:,phasememberAt(c,i,e))
do mySource = 1_pInt, phase_Nsources(phaseAt(c,i,e))
sourceState(phaseAt(c,i,e))%p(mySource)%dotState( :,phasememberAt(c,i,e)) = &
sourceState(phaseAt(c,i,e))%p(mySource)%dotState_backup(:,phasememberAt(c,i,e))
enddo
crystallite_Fp(1:3,1:3,c,i,e) = crystallite_subFp0(1:3,1:3,c,i,e)
crystallite_Fi(1:3,1:3,c,i,e) = crystallite_subFi0(1:3,1:3,c,i,e)
crystallite_Fe(1:3,1:3,c,i,e) = crystallite_subFe0(1:3,1:3,c,i,e)
crystallite_Lp(1:3,1:3,c,i,e) = crystallite_subLp0(1:3,1:3,c,i,e)
crystallite_Li(1:3,1:3,c,i,e) = crystallite_subLi0(1:3,1:3,c,i,e)
crystallite_Tstar_v(1:6,c,i,e) = crystallite_subTstar0_v(1:6,c,i,e)
enddo; enddo
enddo
end select
! --- PERTURB EITHER FORWARD OR BACKWARD ---
!why not OMP?
do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNcomponents = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do c = 1,myNcomponents
crystallite_subF(1:3,1:3,c,i,e) = F_backup(1:3,1:3,c,i,e)
crystallite_subF(k,l,c,i,e) = crystallite_subF(k,l,c,i,e) + myPert
crystallite_todo(c,i,e) = crystallite_requested(c,i,e) &
.and. convergenceFlag_backup(c,i,e)
if (crystallite_todo(c,i,e)) crystallite_converged(c,i,e) = .false. ! start out non-converged
enddo; enddo; enddo
select case(numerics_integrator(numerics_integrationMode))
case(1_pInt)
call crystallite_integrateStateFPI()
case(2_pInt)
call crystallite_integrateStateEuler()
case(3_pInt)
call crystallite_integrateStateAdaptiveEuler()
case(4_pInt)
call crystallite_integrateStateRK4()
case(5_pInt)
call crystallite_integrateStateRKCK45()
end select
!why not OMP?
elementLooping8: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNcomponents = homogenization_Ngrains(mesh_element(3,e))
select case(perturbation)
case(1_pInt)
forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), c = 1:myNcomponents, &
crystallite_requested(c,i,e) .and. crystallite_converged(c,i,e)) & ! converged state warrants stiffness update
dPdF_perturbation1(1:3,1:3,k,l,c,i,e) = &
(crystallite_P(1:3,1:3,c,i,e) - P_backup(1:3,1:3,c,i,e)) / myPert ! tangent dP_ij/dFg_kl
case(2_pInt)
forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), c = 1:myNcomponents, &
crystallite_requested(c,i,e) .and. crystallite_converged(c,i,e)) & ! converged state warrants stiffness update
dPdF_perturbation2(1:3,1:3,k,l,c,i,e) = &
(crystallite_P(1:3,1:3,c,i,e) - P_backup(1:3,1:3,c,i,e)) / myPert ! tangent dP_ij/dFg_kl
end select
enddo elementLooping8
enddo; enddo ! k,l component perturbation loop
endif
enddo pertubationLoop
! --- STIFFNESS ACCORDING TO PERTURBATION METHOD AND CONVERGENCE ---
elementLooping9: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNcomponents = homogenization_Ngrains(mesh_element(3,e))
select case(pert_method)
case(1_pInt)
forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), c = 1:myNcomponents, &
crystallite_requested(c,i,e) .and. convergenceFlag_backup(c,i,e)) & ! perturbation mode 1: central solution converged
crystallite_dPdF(1:3,1:3,1:3,1:3,c,i,e) = dPdF_perturbation1(1:3,1:3,1:3,1:3,c,i,e)
case(2_pInt)
forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), c = 1:myNcomponents, &
crystallite_requested(c,i,e) .and. convergenceFlag_backup(c,i,e)) & ! perturbation mode 2: central solution converged
crystallite_dPdF(1:3,1:3,1:3,1:3,c,i,e) = dPdF_perturbation2(1:3,1:3,1:3,1:3,c,i,e)
case(3_pInt)
forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), c = 1:myNcomponents, &
crystallite_requested(c,i,e) .and. convergenceFlag_backup(c,i,e)) & ! perturbation mode 3: central solution converged
crystallite_dPdF(1:3,1:3,1:3,1:3,c,i,e) = 0.5_pReal* ( dPdF_perturbation1(1:3,1:3,1:3,1:3,c,i,e) &
+ dPdF_perturbation2(1:3,1:3,1:3,1:3,c,i,e))
end select
forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), c = 1:myNcomponents, &
crystallite_requested(c,i,e) .and. .not. convergenceFlag_backup(c,i,e)) & ! for any pertubation mode: if central solution did not converge...
crystallite_dPdF(1:3,1:3,1:3,1:3,c,i,e) = crystallite_fallbackdPdF(1:3,1:3,1:3,1:3,c,i,e) ! ...use (elastic) fallback
enddo elementLooping9
! --- RESTORE ---
!why not OMP?
elementLooping10: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNcomponents = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e); do c = 1,myNcomponents
plasticState (phaseAt(c,i,e))%state( :,phasememberAt(c,i,e)) = &
plasticState (phaseAt(c,i,e))%state_backup(:,phasememberAt(c,i,e))
do mySource = 1_pInt, phase_Nsources(phaseAt(c,i,e))
sourceState(phaseAt(c,i,e))%p(mySource)%state( :,phasememberAt(c,i,e)) = &
sourceState(phaseAt(c,i,e))%p(mySource)%state_backup(:,phasememberAt(c,i,e))
enddo
plasticState (phaseAt(c,i,e))%dotState( :,phasememberAt(c,i,e)) = &
plasticState (phaseAt(c,i,e))%dotState_backup(:,phasememberAt(c,i,e))
do mySource = 1_pInt, phase_Nsources(phaseAt(c,i,e))
sourceState(phaseAt(c,i,e))%p(mySource)%dotState( :,phasememberAt(c,i,e)) = &
sourceState(phaseAt(c,i,e))%p(mySource)%dotState_backup(:,phasememberAt(c,i,e))
enddo
crystallite_subF(1:3,1:3,c,i,e) = F_backup(1:3,1:3,c,i,e)
crystallite_Fp(1:3,1:3,c,i,e) = Fp_backup(1:3,1:3,c,i,e)
crystallite_invFp(1:3,1:3,c,i,e) = InvFp_backup(1:3,1:3,c,i,e)
crystallite_Fi(1:3,1:3,c,i,e) = Fi_backup(1:3,1:3,c,i,e)
crystallite_invFi(1:3,1:3,c,i,e) = InvFi_backup(1:3,1:3,c,i,e)
crystallite_Fe(1:3,1:3,c,i,e) = Fe_backup(1:3,1:3,c,i,e)
crystallite_Lp(1:3,1:3,c,i,e) = Lp_backup(1:3,1:3,c,i,e)
crystallite_Li(1:3,1:3,c,i,e) = Li_backup(1:3,1:3,c,i,e)
crystallite_Tstar_v(1:6,c,i,e) = Tstar_v_backup(1:6,c,i,e)
crystallite_P(1:3,1:3,c,i,e) = P_backup(1:3,1:3,c,i,e)
crystallite_converged(c,i,e) = convergenceFlag_backup(c,i,e)
enddo; enddo
enddo elementLooping10
deallocate(dPdF_perturbation1)
deallocate(dPdF_perturbation2)
deallocate(F_backup )
deallocate(Fp_backup )
deallocate(InvFp_backup )
deallocate(Fi_backup )
deallocate(InvFi_backup )
deallocate(Fe_backup )
deallocate(Lp_backup )
deallocate(Li_backup )
deallocate(P_backup )
deallocate(Tstar_v_backup )
deallocate(convergenceFlag_backup)
endif jacobianMethod
endif computeJacobian endif computeJacobian
!why not OMP? !why not OMP?