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

* new global integer variable "numerics_integrationMode" (1 indicating integration of central solution, 2 indicating integration of perturbed state) 

* combined "integrator" and "integratorStiffness" in new global variable "numerics_integrator"
This commit is contained in:
Christoph Kords 2011-02-23 08:29:51 +00:00
parent d835380bc0
commit 96d3682d5e
3 changed files with 80 additions and 77 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

38
code/constitutive.f90
View File

@ -50,7 +50,7 @@ subroutine constitutive_init()
!************************************** !**************************************
use prec, only: pReal,pInt use prec, only: pReal,pInt
use debug, only: debugger, selectiveDebugger, debug_e, debug_i, debug_g use debug, only: debugger, selectiveDebugger, debug_e, debug_i, debug_g
use numerics, only: integrator, integratorStiffness use numerics, only: numerics_integrator
use IO, only: IO_error, IO_open_file, IO_open_jobFile use IO, only: IO_error, IO_open_file, IO_open_jobFile
use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips
use material use material
@ -129,13 +129,13 @@ subroutine constitutive_init()
allocate(constitutive_sizeDotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizeDotState = 0_pInt allocate(constitutive_sizeDotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizeDotState = 0_pInt
allocate(constitutive_sizeState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizeState = 0_pInt allocate(constitutive_sizeState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizeState = 0_pInt
allocate(constitutive_sizePostResults(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)); constitutive_sizePostResults = 0_pInt allocate(constitutive_sizePostResults(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)); constitutive_sizePostResults = 0_pInt
if (integrator == 1 .or. integratorStiffness == 1) then if (any(numerics_integrator == 1)) then
allocate(constitutive_previousDotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(constitutive_previousDotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(constitutive_previousDotState2(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(constitutive_previousDotState2(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
endif endif
if (integrator == 4 .or. integratorStiffness == 4) & if (any(numerics_integrator == 4)) &
allocate(constitutive_RK4dotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(constitutive_RK4dotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
if (integrator == 5 .or. integratorStiffness == 5) & if (any(numerics_integrator == 5)) &
allocate(constitutive_RKCK45dotState(6,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(constitutive_RKCK45dotState(6,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
do e = 1,mesh_NcpElems ! loop over elements do e = 1,mesh_NcpElems ! loop over elements
@ -154,13 +154,13 @@ subroutine constitutive_init()
allocate(constitutive_aTolState(g,i,e)%p(constitutive_j2_sizeState(myInstance))) allocate(constitutive_aTolState(g,i,e)%p(constitutive_j2_sizeState(myInstance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
if (integrator == 1 .or. integratorStiffness == 1) then if (any(numerics_integrator == 1)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) allocate(constitutive_previousDotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
endif endif
if (integrator == 4 .or. integratorStiffness == 4) & if (any(numerics_integrator == 4)) &
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
if (integrator == 5 .or. integratorStiffness == 5) then if (any(numerics_integrator == 5)) then
do s = 1,6 do s = 1,6
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
enddo enddo
@ -180,13 +180,13 @@ subroutine constitutive_init()
allocate(constitutive_aTolState(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) allocate(constitutive_aTolState(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
if (integrator == 1 .or. integratorStiffness == 1) then if (any(numerics_integrator == 1)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) allocate(constitutive_previousDotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
endif endif
if (integrator == 4 .or. integratorStiffness == 4) & if (any(numerics_integrator == 4)) &
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
if (integrator == 5 .or. integratorStiffness == 5) then if (any(numerics_integrator == 5)) then
do s = 1,6 do s = 1,6
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
enddo enddo
@ -206,13 +206,13 @@ subroutine constitutive_init()
allocate(constitutive_aTolState(g,i,e)%p(constitutive_titanmod_sizeState(myInstance))) allocate(constitutive_aTolState(g,i,e)%p(constitutive_titanmod_sizeState(myInstance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
if (integrator == 1 .or. integratorStiffness == 1) then if (any(numerics_integrator == 1)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) allocate(constitutive_previousDotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
endif endif
if (integrator == 4 .or. integratorStiffness == 4) & if (any(numerics_integrator == 4)) &
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
if (integrator == 5 .or. integratorStiffness == 5) then if (any(numerics_integrator == 5)) then
do s = 1,6 do s = 1,6
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
enddo enddo
@ -232,13 +232,13 @@ subroutine constitutive_init()
allocate(constitutive_aTolState(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance))) allocate(constitutive_aTolState(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
if (integrator == 1 .or. integratorStiffness == 1) then if (any(numerics_integrator == 1)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) allocate(constitutive_previousDotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
endif endif
if (integrator == 4 .or. integratorStiffness == 4) & if (any(numerics_integrator == 4)) &
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
if (integrator == 5 .or. integratorStiffness == 5) then if (any(numerics_integrator == 5)) then
do s = 1,6 do s = 1,6
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
enddo enddo
@ -258,13 +258,13 @@ subroutine constitutive_init()
allocate(constitutive_aTolState(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) allocate(constitutive_aTolState(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
if (integrator == 1 .or. integratorStiffness == 1) then if (any(numerics_integrator == 1)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) allocate(constitutive_previousDotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
endif endif
if (integrator == 4 .or. integratorStiffness == 4) & if (any(numerics_integrator == 4)) &
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
if (integrator == 5 .or. integratorStiffness == 5) then if (any(numerics_integrator == 5)) then
do s = 1,6 do s = 1,6
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
enddo enddo

2
code/constitutive_nonlocal.f90
View File

@ -1664,7 +1664,7 @@ if (.not. phase_localConstitution(material_phase(g,ip,el))) then
* constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el) ** 2.0_pReal * constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el) ** 2.0_pReal
where (constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el) < 0.0_pReal) & ! ..negative compatibility... where (constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el) < 0.0_pReal) & ! ..negative compatibility...
rhoDotFlux(1:ns,topp) = rhoDotFlux(1:ns,topp) + lineLength / mesh_ipVolume(ip,el) & ! ... transferring to opposite signed dislocation type rhoDotFlux(1:ns,topp) = rhoDotFlux(1:ns,topp) + lineLength / mesh_ipVolume(ip,el) & ! ... transferring to opposite signed dislocation type
* constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el) ** 2.0_pReal * constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el) ** 2.0_pReal
endif endif
enddo enddo
enddo enddo

117
code/crystallite.f90
View File

@ -88,9 +88,7 @@ use prec, only: pInt, &
pReal pReal
use debug, only: debug_info, & use debug, only: debug_info, &
debug_reset debug_reset
use numerics, only: integrator, & use numerics, only: subStepSizeCryst, &
integratorStiffness, &
subStepSizeCryst, &
stepIncreaseCryst stepIncreaseCryst
use math, only: math_I3, & use math, only: math_I3, &
math_EulerToR, & math_EulerToR, &
@ -425,8 +423,8 @@ use numerics, only: subStepMinCryst, &
pert_method, & pert_method, &
nCryst, & nCryst, &
iJacoStiffness, & iJacoStiffness, &
integratorStiffness, & numerics_integrator, &
integrator numerics_integrationMode
use debug, only: debugger, & use debug, only: debugger, &
selectiveDebugger, & selectiveDebugger, &
verboseDebugger, & verboseDebugger, &
@ -540,6 +538,7 @@ crystallite_subStep = 0.0_pReal
! --+>> CRYSTALLITE CUTBACK LOOP <<+-- ! --+>> CRYSTALLITE CUTBACK LOOP <<+--
NiterationCrystallite = 0_pInt NiterationCrystallite = 0_pInt
numerics_integrationMode = 1_pInt
do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMinCryst)) ! cutback loop for crystallites do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMinCryst)) ! cutback loop for crystallites
!$OMP PARALLEL DO PRIVATE(myNgrains,formerSubStep) !$OMP PARALLEL DO PRIVATE(myNgrains,formerSubStep)
@ -617,17 +616,17 @@ do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2
! --- integrate --- ! --- integrate ---
if (any(crystallite_todo)) then if (any(crystallite_todo)) then
select case(integrator) select case(numerics_integrator(numerics_integrationMode))
case(1) case(1)
call crystallite_integrateStateFPI(1) call crystallite_integrateStateFPI()
case(2) case(2)
call crystallite_integrateStateEuler(1) call crystallite_integrateStateEuler()
case(3) case(3)
call crystallite_integrateStateAdaptiveEuler(1) call crystallite_integrateStateAdaptiveEuler()
case(4) case(4)
call crystallite_integrateStateRK4(1) call crystallite_integrateStateRK4()
case(5) case(5)
call crystallite_integrateStateRKCK45(1) call crystallite_integrateStateRKCK45()
endselect endselect
endif endif
@ -669,6 +668,7 @@ enddo
! --+>> STIFFNESS CALCULATION <<+-- ! --+>> STIFFNESS CALCULATION <<+--
if(updateJaco) then ! Jacobian required if(updateJaco) then ! Jacobian required
numerics_integrationMode = 2_pInt
! --- BACKUP --- ! --- BACKUP ---
@ -714,17 +714,17 @@ if(updateJaco) then
crystallite_todo = crystallite_requested .and. crystallite_converged crystallite_todo = crystallite_requested .and. crystallite_converged
where (crystallite_todo) crystallite_converged = .false. ! start out non-converged where (crystallite_todo) crystallite_converged = .false. ! start out non-converged
select case(integratorStiffness) select case(numerics_integrator(numerics_integrationMode))
case(1) case(1)
call crystallite_integrateStateFPI(2) call crystallite_integrateStateFPI()
case(2) case(2)
call crystallite_integrateStateEuler(2) call crystallite_integrateStateEuler()
case(3) case(3)
call crystallite_integrateStateAdaptiveEuler(2) call crystallite_integrateStateAdaptiveEuler()
case(4) case(4)
call crystallite_integrateStateRK4(2) call crystallite_integrateStateRK4()
case(5) case(5)
call crystallite_integrateStateRKCK45(2) call crystallite_integrateStateRKCK45()
end select end select
!OMP PARALLEL DO PRIVATE(myNgrains) !OMP PARALLEL DO PRIVATE(myNgrains)
@ -804,11 +804,12 @@ endsubroutine
! integrate stress, state and Temperature with ! integrate stress, state and Temperature with
! 4h order explicit Runge Kutta method ! 4h order explicit Runge Kutta method
!******************************************************************** !********************************************************************
subroutine crystallite_integrateStateRK4(mode,gg,ii,ee) subroutine crystallite_integrateStateRK4(gg,ii,ee)
!*** variables and functions from other modules ***! !*** variables and functions from other modules ***!
use prec, only: pInt, & use prec, only: pInt, &
pReal pReal
use numerics, only: numerics_integrationMode
use debug, only: debugger, & use debug, only: debugger, &
selectiveDebugger, & selectiveDebugger, &
verboseDebugger, & verboseDebugger, &
@ -838,7 +839,6 @@ real(pReal), dimension(4), parameter :: timeStepFraction = (/0.5_pReal, 0.
real(pReal), dimension(4), parameter :: weight = (/1.0_pReal, 2.0_pReal, 2.0_pReal, 1.0_pReal/) ! factor giving the fraction of the original timestep used for Runge Kutta Integration real(pReal), dimension(4), parameter :: weight = (/1.0_pReal, 2.0_pReal, 2.0_pReal, 1.0_pReal/) ! factor giving the fraction of the original timestep used for Runge Kutta Integration
!*** input variables ***! !*** input variables ***!
integer(pInt), intent(in) :: mode ! mode of calculation; 1: central solution, 2: stiffness (by perturbation)
integer(pInt), optional, intent(in):: ee, & ! element index integer(pInt), optional, intent(in):: ee, & ! element index
ii, & ! integration point index ii, & ! integration point index
gg ! grain index gg ! grain index
@ -957,7 +957,7 @@ do n = 1,4
!$OMP DO !$OMP DO
do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
if (crystallite_integrateStress(mode,g,i,e,timeStepFraction(n))) then ! fraction of original times step if (crystallite_integrateStress(g,i,e,timeStepFraction(n))) then ! fraction of original times step
if (n == 4) then ! final integration step if (n == 4) then ! final integration step
if (verboseDebugger .and. selectiveDebugger .and. e == debug_e .and. i == debug_i .and. g == debug_g) then if (verboseDebugger .and. selectiveDebugger .and. e == debug_e .and. i == debug_i .and. g == debug_g) then
mySizeDotState = constitutive_sizeDotState(g,i,e) mySizeDotState = constitutive_sizeDotState(g,i,e)
@ -973,7 +973,7 @@ do n = 1,4
crystallite_converged(g,i,e) = .true. ! ... converged per definition crystallite_converged(g,i,e) = .true. ! ... converged per definition
crystallite_todo(g,i,e) = .false. ! ... integration done crystallite_todo(g,i,e) = .false. ! ... integration done
!$OMP CRITICAL (distributionState) !$OMP CRITICAL (distributionState)
debug_StateLoopDistribution(n,mode) = debug_StateLoopDistribution(n,mode) + 1 debug_StateLoopDistribution(n,numerics_integrationMode) = debug_StateLoopDistribution(n,numerics_integrationMode) + 1
!$OMP END CRITICAL (distributionState) !$OMP END CRITICAL (distributionState)
endif endif
else ! broken stress integration else ! broken stress integration
@ -1041,7 +1041,7 @@ endsubroutine
! 5th order Runge-Kutta Cash-Karp method with adaptive step size ! 5th order Runge-Kutta Cash-Karp method with adaptive step size
! (use 5th order solution to advance = "local extrapolation") ! (use 5th order solution to advance = "local extrapolation")
!******************************************************************** !********************************************************************
subroutine crystallite_integrateStateRKCK45(mode,gg,ii,ee) subroutine crystallite_integrateStateRKCK45(gg,ii,ee)
!*** variables and functions from other modules ***! !*** variables and functions from other modules ***!
use prec, only: pInt, & use prec, only: pInt, &
@ -1056,7 +1056,8 @@ use debug, only: debugger, &
use numerics, only: rTol_crystalliteState, & use numerics, only: rTol_crystalliteState, &
rTol_crystalliteTemperature, & rTol_crystalliteTemperature, &
subStepSizeCryst, & subStepSizeCryst, &
stepIncreaseCryst stepIncreaseCryst, &
numerics_integrationMode
use FEsolving, only: FEsolving_execElem, & use FEsolving, only: FEsolving_execElem, &
FEsolving_execIP, & FEsolving_execIP, &
theInc theInc
@ -1080,7 +1081,6 @@ implicit none
!*** input variables ***! !*** input variables ***!
integer(pInt), intent(in) :: mode ! mode of calculation; 1: central solution, 2: stiffness (by perturbation)
integer(pInt), optional, intent(in):: ee, & ! element index integer(pInt), optional, intent(in):: ee, & ! element index
ii, & ! integration point index ii, & ! integration point index
gg ! grain index gg ! grain index
@ -1261,7 +1261,7 @@ do n = 1,5
!$OMP DO !$OMP DO
do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
if (.not. crystallite_integrateStress(mode,g,i,e,c(n))) then ! fraction of original time step if (.not. crystallite_integrateStress(g,i,e,c(n))) then ! fraction of original time step
if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local...
!$OMP CRITICAL (checkTodo) !$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped
@ -1404,11 +1404,11 @@ relTemperatureResiduum = 0.0_pReal
!$OMP DO !$OMP DO
do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
if (crystallite_integrateStress(mode,g,i,e)) then if (crystallite_integrateStress(g,i,e)) then
crystallite_converged(g,i,e) = .true. ! ... converged per definitionem crystallite_converged(g,i,e) = .true. ! ... converged per definitionem
crystallite_todo(g,i,e) = .false. ! ... integration done crystallite_todo(g,i,e) = .false. ! ... integration done
!$OMP CRITICAL (distributionState) !$OMP CRITICAL (distributionState)
debug_StateLoopDistribution(6,mode) = debug_StateLoopDistribution(6,mode) + 1 debug_StateLoopDistribution(6,numerics_integrationMode) = debug_StateLoopDistribution(6,numerics_integrationMode) + 1
!$OMP END CRITICAL (distributionState) !$OMP END CRITICAL (distributionState)
else else
if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local...
@ -1425,7 +1425,7 @@ relTemperatureResiduum = 0.0_pReal
! --- nonlocal convergence check --- ! --- nonlocal convergence check ---
if (verboseDebugger .and. mode==1) write(6,*) 'crystallite_converged',crystallite_converged if (verboseDebugger .and. numerics_integrationMode == 1) write(6,*) 'crystallite_converged',crystallite_converged
if (.not. singleRun) then ! if not requesting Integration of just a single IP if (.not. singleRun) then ! if not requesting Integration of just a single IP
if ( any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) then ! any non-local not yet converged (or broken)... if ( any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) then ! any non-local not yet converged (or broken)...
crystallite_converged = crystallite_converged .and. crystallite_localConstitution ! ...restart all non-local as not converged crystallite_converged = crystallite_converged .and. crystallite_localConstitution ! ...restart all non-local as not converged
@ -1441,7 +1441,7 @@ endsubroutine
! integrate stress, state and Temperature with ! integrate stress, state and Temperature with
! 1nd order Euler method with adaptive step size ! 1nd order Euler method with adaptive step size
!******************************************************************** !********************************************************************
subroutine crystallite_integrateStateAdaptiveEuler(mode,gg,ii,ee) subroutine crystallite_integrateStateAdaptiveEuler(gg,ii,ee)
!*** variables and functions from other modules ***! !*** variables and functions from other modules ***!
use prec, only: pInt, & use prec, only: pInt, &
@ -1456,7 +1456,8 @@ use debug, only: debugger, &
use numerics, only: rTol_crystalliteState, & use numerics, only: rTol_crystalliteState, &
rTol_crystalliteTemperature, & rTol_crystalliteTemperature, &
subStepSizeCryst, & subStepSizeCryst, &
stepIncreaseCryst stepIncreaseCryst, &
numerics_integrationMode
use FEsolving, only: FEsolving_execElem, & use FEsolving, only: FEsolving_execElem, &
FEsolving_execIP FEsolving_execIP
use mesh, only: mesh_element, & use mesh, only: mesh_element, &
@ -1478,7 +1479,6 @@ implicit none
!*** input variables ***! !*** input variables ***!
integer(pInt), intent(in) :: mode ! mode of calculation; 1: central solution, 2: stiffness (by perturbation)
integer(pInt), optional, intent(in):: ee, & ! element index integer(pInt), optional, intent(in):: ee, & ! element index
ii, & ! integration point index ii, & ! integration point index
gg ! grain index gg ! grain index
@ -1598,7 +1598,7 @@ stateResiduum = 0.0_pReal
!$OMP DO !$OMP DO
do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
if (.not. crystallite_integrateStress(mode,g,i,e)) then if (.not. crystallite_integrateStress(g,i,e)) then
if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local...
!$OMP CRITICAL (checkTodo) !$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped
@ -1690,7 +1690,7 @@ relTemperatureResiduum = 0.0_pReal
crystallite_converged(g,i,e) = .true. ! ... converged per definitionem crystallite_converged(g,i,e) = .true. ! ... converged per definitionem
crystallite_todo(g,i,e) = .false. ! ... integration done crystallite_todo(g,i,e) = .false. ! ... integration done
!$OMP CRITICAL (distributionState) !$OMP CRITICAL (distributionState)
debug_StateLoopDistribution(2,mode) = debug_StateLoopDistribution(2,mode) + 1 debug_StateLoopDistribution(2,numerics_integrationMode) = debug_StateLoopDistribution(2,numerics_integrationMode) + 1
!$OMP END CRITICAL (distributionState) !$OMP END CRITICAL (distributionState)
endif endif
@ -1703,7 +1703,7 @@ relTemperatureResiduum = 0.0_pReal
! --- NONLOCAL CONVERGENCE CHECK --- ! --- NONLOCAL CONVERGENCE CHECK ---
if (verboseDebugger .and. mode==1) write(6,*) 'crystallite_converged',crystallite_converged if (verboseDebugger .and. numerics_integrationMode==1_pInt) write(6,*) 'crystallite_converged',crystallite_converged
if (.not. singleRun) then ! if not requesting Integration of just a single IP if (.not. singleRun) then ! if not requesting Integration of just a single IP
if ( any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) then ! any non-local not yet converged (or broken)... if ( any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) then ! any non-local not yet converged (or broken)...
crystallite_converged = crystallite_converged .and. crystallite_localConstitution ! ...restart all non-local as not converged crystallite_converged = crystallite_converged .and. crystallite_localConstitution ! ...restart all non-local as not converged
@ -1718,11 +1718,12 @@ endsubroutine
! integrate stress, state and Temperature with ! integrate stress, state and Temperature with
! 1st order explicit Euler method ! 1st order explicit Euler method
!******************************************************************** !********************************************************************
subroutine crystallite_integrateStateEuler(mode,gg,ii,ee) subroutine crystallite_integrateStateEuler(gg,ii,ee)
!*** variables and functions from other modules ***! !*** variables and functions from other modules ***!
use prec, only: pInt, & use prec, only: pInt, &
pReal pReal
use numerics, only: numerics_integrationMode
use debug, only: debugger, & use debug, only: debugger, &
selectiveDebugger, & selectiveDebugger, &
verboseDebugger, & verboseDebugger, &
@ -1746,7 +1747,6 @@ use constitutive, only: constitutive_sizeDotState, &
implicit none implicit none
!*** input variables ***! !*** input variables ***!
integer(pInt), intent(in) :: mode ! mode of calculation; 1: central solution, 2: stiffness (by perturbation)
integer(pInt), optional, intent(in):: ee, & ! element index integer(pInt), optional, intent(in):: ee, & ! element index
ii, & ! integration point index ii, & ! integration point index
gg ! grain index gg ! grain index
@ -1863,10 +1863,10 @@ endif
!$OMP DO !$OMP DO
do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
if (crystallite_integrateStress(mode,g,i,e)) then if (crystallite_integrateStress(g,i,e)) then
crystallite_converged(g,i,e) = .true. crystallite_converged(g,i,e) = .true.
!$OMP CRITICAL (distributionState) !$OMP CRITICAL (distributionState)
debug_StateLoopDistribution(1,mode) = debug_StateLoopDistribution(1,mode) + 1 debug_StateLoopDistribution(1,numerics_integrationMode) = debug_StateLoopDistribution(1,numerics_integrationMode) + 1
!$OMP END CRITICAL (distributionState) !$OMP END CRITICAL (distributionState)
else ! broken stress integration else ! broken stress integration
if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local...
@ -1900,7 +1900,7 @@ endsubroutine
! adaptive 1st order explicit Euler method ! adaptive 1st order explicit Euler method
! using Fixed Point Iteration to adapt the stepsize ! using Fixed Point Iteration to adapt the stepsize
!******************************************************************** !********************************************************************
subroutine crystallite_integrateStateFPI(mode,gg,ii,ee) subroutine crystallite_integrateStateFPI(gg,ii,ee)
!*** variables and functions from other modules ***! !*** variables and functions from other modules ***!
use prec, only: pInt, & use prec, only: pInt, &
@ -1912,7 +1912,8 @@ use debug, only: debugger, &
debug_i, & debug_i, &
debug_g, & debug_g, &
debug_StateLoopDistribution debug_StateLoopDistribution
use numerics, only: nState use numerics, only: nState, &
numerics_integrationMode
use FEsolving, only: FEsolving_execElem, & use FEsolving, only: FEsolving_execElem, &
FEsolving_execIP FEsolving_execIP
use mesh, only: mesh_element, & use mesh, only: mesh_element, &
@ -1930,7 +1931,6 @@ use constitutive, only: constitutive_sizeDotState, &
implicit none implicit none
!*** input variables ***! !*** input variables ***!
integer(pInt), intent(in) :: mode ! mode of calculation; 1: central solution, 2: stiffness (by perturbation)
integer(pInt), optional, intent(in):: ee, & ! element index integer(pInt), optional, intent(in):: ee, & ! element index
ii, & ! integration point index ii, & ! integration point index
gg ! grain index gg ! grain index
@ -2049,7 +2049,7 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
!$OMP DO !$OMP DO
do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
crystallite_todo(g,i,e) = crystallite_integrateStress(mode,g,i,e) crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e)
if ( .not. crystallite_localConstitution(g,i,e) .and. .not. crystallite_todo(g,i,e)) then ! if broken non-local... if ( .not. crystallite_localConstitution(g,i,e) .and. .not. crystallite_todo(g,i,e)) then ! if broken non-local...
!$OMP CRITICAL (checkTodo) !$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped
@ -2059,7 +2059,7 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
enddo; enddo; enddo enddo; enddo; enddo
!$OMP ENDDO !$OMP ENDDO
if (verboseDebugger .and. mode == 1) then if (verboseDebugger .and. numerics_integrationMode == 1) then
!$OMP SINGLE !$OMP SINGLE
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) count(crystallite_todo(:,:,:)),'grains todo after stress integration' write(6,*) count(crystallite_todo(:,:,:)),'grains todo after stress integration'
@ -2112,7 +2112,8 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
!$OMP END CRITICAL (checkTodo) !$OMP END CRITICAL (checkTodo)
elseif (crystallite_converged(g,i,e)) then elseif (crystallite_converged(g,i,e)) then
!$OMP CRITICAL (distributionState) !$OMP CRITICAL (distributionState)
debug_StateLoopDistribution(NiterationState,mode) = debug_StateLoopDistribution(NiterationState,mode) + 1 debug_StateLoopDistribution(NiterationState,numerics_integrationMode) = &
debug_StateLoopDistribution(NiterationState,numerics_integrationMode) + 1
!$OMP END CRITICAL (distributionState) !$OMP END CRITICAL (distributionState)
endif endif
endif endif
@ -2135,7 +2136,7 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
!$OMP ENDDO !$OMP ENDDO
if (verboseDebugger .and. mode == 1) then if (verboseDebugger .and. numerics_integrationMode == 1) then
!$OMP SINGLE !$OMP SINGLE
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) count(crystallite_converged(:,:,:)),'grains converged after state integration no.', NiterationState write(6,*) count(crystallite_converged(:,:,:)),'grains converged after state integration no.', NiterationState
@ -2155,7 +2156,7 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
crystallite_todo = crystallite_todo .and. .not. crystallite_converged ! skip all converged crystallite_todo = crystallite_todo .and. .not. crystallite_converged ! skip all converged
if (verboseDebugger .and. mode == 1) then if (verboseDebugger .and. numerics_integrationMode == 1) then
!$OMP SINGLE !$OMP SINGLE
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) count(crystallite_converged(:,:,:)),'grains converged after non-local check' write(6,*) count(crystallite_converged(:,:,:)),'grains converged after non-local check'
@ -2183,7 +2184,8 @@ function crystallite_updateState(g,i,e)
use prec, only: pReal, & use prec, only: pReal, &
pInt, & pInt, &
pLongInt pLongInt
use numerics, only: rTol_crystalliteState use numerics, only: rTol_crystalliteState, &
numerics_integrationMode
use constitutive, only: constitutive_dotState, & use constitutive, only: constitutive_dotState, &
constitutive_previousDotState, & constitutive_previousDotState, &
constitutive_sizeDotState, & constitutive_sizeDotState, &
@ -2236,7 +2238,7 @@ constitutive_state(g,i,e)%p(1:mySize) = constitutive_state(g,i,e)%p(1:mySize) -
crystallite_updateState = all( abs(residuum) < constitutive_aTolState(g,i,e)%p(1:mySize) & crystallite_updateState = all( abs(residuum) < constitutive_aTolState(g,i,e)%p(1:mySize) &
.or. abs(residuum) < rTol_crystalliteState*abs(constitutive_state(g,i,e)%p(1:mySize)) ) .or. abs(residuum) < rTol_crystalliteState*abs(constitutive_state(g,i,e)%p(1:mySize)) )
if (verboseDebugger .and. (e == debug_e .and. i == debug_i .and. g == debug_g)) then if (verboseDebugger .and. (e == debug_e .and. i == debug_i .and. g == debug_g) .and. numerics_integrationMode == 1_pInt) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
if (crystallite_updateState) then if (crystallite_updateState) then
write(6,*) '::: updateState converged',g,i,e write(6,*) '::: updateState converged',g,i,e
@ -2323,7 +2325,6 @@ endfunction
!*** acceleration of the Newton-Raphson correction *** !*** acceleration of the Newton-Raphson correction ***
!*********************************************************************** !***********************************************************************
function crystallite_integrateStress(& function crystallite_integrateStress(&
mode, & ! 1: central solution, 2: stiffness (by perturbation)
g,& ! grain number g,& ! grain number
i,& ! integration point number i,& ! integration point number
e,& ! element number e,& ! element number
@ -2339,7 +2340,8 @@ function crystallite_integrateStress(&
aTol_crystalliteStress, & aTol_crystalliteStress, &
rTol_crystalliteStress, & rTol_crystalliteStress, &
iJacoLpresiduum, & iJacoLpresiduum, &
relevantStrain relevantStrain, &
numerics_integrationMode
use debug, only: debugger, & use debug, only: debugger, &
debug_g, & debug_g, &
debug_i, & debug_i, &
@ -2368,8 +2370,7 @@ function crystallite_integrateStress(&
implicit none implicit none
!*** input variables ***! !*** input variables ***!
integer(pInt), intent(in):: mode, & ! 1 or 2 integer(pInt), intent(in):: e, & ! element index
e, & ! element index
i, & ! integration point index i, & ! integration point index
g ! grain index g ! grain index
real(pReal), optional, intent(in) :: fraction ! fraction of timestep real(pReal), optional, intent(in) :: fraction ! fraction of timestep
@ -2502,7 +2503,7 @@ LpLoop: do
if (tock < tick) debug_cumLpTicks = debug_cumLpTicks + maxticks if (tock < tick) debug_cumLpTicks = debug_cumLpTicks + maxticks
!$OMP END CRITICAL (debugTimingLpTangent) !$OMP END CRITICAL (debugTimingLpTangent)
if (verboseDebugger .and. (e == debug_e .and. i == debug_i .and. g == debug_g)) then if (verboseDebugger .and. (e == debug_e .and. i == debug_i .and. g == debug_g) .and. numerics_integrationMode == 1_pInt) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(a,i3,x,i2,x,i5,x,a,x,i3)') '::: integrateStress at ' ,g,i,e, ' ; iteration ', NiterationStress write(6,'(a,i3,x,i2,x,i5,x,a,x,i3)') '::: integrateStress at ' ,g,i,e, ' ; iteration ', NiterationStress
write(6,*) write(6,*)
@ -2557,7 +2558,8 @@ LpLoop: do
residuum = residuum_old residuum = residuum_old
!$OMP CRITICAL (distributionLeapfrogBreak) !$OMP CRITICAL (distributionLeapfrogBreak)
debug_LeapfrogBreakDistribution(NiterationStress,mode) = debug_LeapfrogBreakDistribution(NiterationStress,mode) + 1 debug_LeapfrogBreakDistribution(NiterationStress,numerics_integrationMode) = &
debug_LeapfrogBreakDistribution(NiterationStress,numerics_integrationMode) + 1
!$OMP END CRITICAL (distributionLeapfrogBreak) !$OMP END CRITICAL (distributionLeapfrogBreak)
! residuum got better ! residuum got better
@ -2587,7 +2589,7 @@ LpLoop: do
endif endif
return return
else else
if (verboseDebugger .and. (e == debug_e .and. i == debug_i .and. g == debug_g)) then if (verboseDebugger .and. (e == debug_e .and. i == debug_i .and. g == debug_g) .and. numerics_integrationMode==1_pInt) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(a,i3,x,i2,x,i5,x,a,i3)') '::: integrateStress did dR/dLp inversion at ',g,i,e, & write(6,'(a,i3,x,i2,x,i5,x,a,i3)') '::: integrateStress did dR/dLp inversion at ',g,i,e, &
'; iteration ', NiterationStress '; iteration ', NiterationStress
@ -2646,7 +2648,7 @@ LpLoop: do
! set return flag to true ! set return flag to true
crystallite_integrateStress = .true. crystallite_integrateStress = .true.
if (verboseDebugger .and. (e == debug_e .and. i == debug_i .and. g == debug_g)) then if (verboseDebugger .and. (e == debug_e .and. i == debug_i .and. g == debug_g) .and. numerics_integrationMode == 1_pInt) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(a,i3,x,i2,x,i5,x,a,x,i3)') '::: integrateStress converged at ',g,i,e,' ; iteration ', NiterationStress write(6,'(a,i3,x,i2,x,i5,x,a,x,i3)') '::: integrateStress converged at ',g,i,e,' ; iteration ', NiterationStress
write(6,*) write(6,*)
@ -2661,7 +2663,8 @@ LpLoop: do
endif endif
!$OMP CRITICAL (distributionStress) !$OMP CRITICAL (distributionStress)
debug_StressLoopDistribution(NiterationStress,mode) = debug_StressLoopDistribution(NiterationStress,mode) + 1 debug_StressLoopDistribution(NiterationStress,numerics_integrationMode) = &
debug_StressLoopDistribution(NiterationStress,numerics_integrationMode) + 1
!$OMP END CRITICAL (distributionStress) !$OMP END CRITICAL (distributionStress)
return return