remved bug in phenopower law state initialization and a minor change in Makefile

This commit is contained in:
Luv Sharma 2014-07-07 14:21:58 +00:00
parent 7a478d646a
commit 7dfda9d8af
3 changed files with 31 additions and 50 deletions

View File

@ -365,7 +365,7 @@ endif
DAMASK_spectral.exe: DAMASK_spectral_driver.o DAMASK_spectral.exe: DAMASK_spectral_driver.o
$(PREFIX) $(COMPILERNAME) $(OPENMP_FLAG_$(F90)) $(LINK_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) \ $(PREFIX) $(COMPILERNAME) $(OPENMP_FLAG_$(F90)) $(LINK_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) \
-o DAMASK_spectral.exe DAMASK_spectral_driver.o \ -o DAMASK_spectral.exe DAMASK_spectral_driver.o \
$(COMPILED_FILES) $(LIB_DIRS) $(LIBRARIES) $(LIB_DIRS) $(RUN_PATH) $(SUFFIX) $(COMPILED_FILES) $(LIB_DIRS) $(RUN_PATH) $(LIBRARIES) $(SUFFIX)
DAMASK_spectral_driver.o: DAMASK_spectral_driver.f90 DAMASK_spectral_solverBasic.o $(PETSC_FILES) DAMASK_spectral_driver.o: DAMASK_spectral_driver.f90 DAMASK_spectral_solverBasic.o $(PETSC_FILES)
$(PREFIX) $(COMPILERNAME) $(COMPILE_MAXOPTI) -c DAMASK_spectral_driver.f90 $(SUFFIX) $(PREFIX) $(COMPILERNAME) $(COMPILE_MAXOPTI) -c DAMASK_spectral_driver.f90 $(SUFFIX)

View File

@ -514,12 +514,12 @@ subroutine constitutive_phenopowerlaw_init(fileUnit)
sizeDotState = sizeState sizeDotState = sizeState
plasticState(phase)%sizeDotState = sizeState plasticState(phase)%sizeDotState = sizeState
plasticState(phase)%sizePostResults = constitutive_phenopowerlaw_sizePostResults(instance) plasticState(phase)%sizePostResults = constitutive_phenopowerlaw_sizePostResults(instance)
allocate(plasticState(phase)%aTolState (sizeState), source=0.0_pReal) allocate(plasticState(phase)%aTolState ( sizeState), source=0.0_pReal)
allocate(plasticState(phase)%state0 (sizeState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%state0 ( sizeState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%partionedState0 (sizeState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%partionedState0 ( sizeState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%subState0 (sizeState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%subState0 ( sizeState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%state (sizeState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%state ( sizeState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%state_backup ( sizeState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal)
if (any(numerics_integrator == 1_pInt)) then if (any(numerics_integrator == 1_pInt)) then
@ -581,6 +581,7 @@ subroutine constitutive_phenopowerlaw_init(fileUnit)
enddo; enddo enddo; enddo
enddo; enddo enddo; enddo
call constitutive_phenopowerlaw_stateInit(phase,instance) call constitutive_phenopowerlaw_stateInit(phase,instance)
call constitutive_phenopowerlaw_aTolState(phase,instance) call constitutive_phenopowerlaw_aTolState(phase,instance)
endif myPhase2 endif myPhase2
@ -609,7 +610,7 @@ subroutine constitutive_phenopowerlaw_stateInit(ph,instance)
i i
real(pReal), dimension(plasticState(ph)%sizeState) :: & real(pReal), dimension(plasticState(ph)%sizeState) :: &
tempState tempState
tempState = 0.0_pReal
do i = 1_pInt,lattice_maxNslipFamily do i = 1_pInt,lattice_maxNslipFamily
tempState(1+sum(constitutive_phenopowerlaw_Nslip(1:i-1,instance)) : & tempState(1+sum(constitutive_phenopowerlaw_Nslip(1:i-1,instance)) : &
sum(constitutive_phenopowerlaw_Nslip(1:i ,instance))) = & sum(constitutive_phenopowerlaw_Nslip(1:i ,instance))) = &
@ -624,8 +625,7 @@ subroutine constitutive_phenopowerlaw_stateInit(ph,instance)
constitutive_phenopowerlaw_tau0_twin(i,instance) constitutive_phenopowerlaw_tau0_twin(i,instance)
enddo enddo
plasticState(ph)%state0 = spread(tempState,2,size(plasticState(ph)%state0(1,:))) plasticState(ph)%state0(:,:) = spread(tempState,2,size(plasticState(ph)%state0(1,:)))
end subroutine constitutive_phenopowerlaw_stateInit end subroutine constitutive_phenopowerlaw_stateInit
@ -875,7 +875,7 @@ subroutine constitutive_phenopowerlaw_dotState(Tstar_v,ipc,ip,el)
index_F = nSlip + nTwin + 2_pInt index_F = nSlip + nTwin + 2_pInt
offset_accshear_slip = nSlip + nTwin + 2_pInt offset_accshear_slip = nSlip + nTwin + 2_pInt
offset_accshear_twin = nSlip + nTwin + 2_pInt + nSlip offset_accshear_twin = nSlip + nTwin + 2_pInt + nSlip
plasticState(ph)%dotState = 0.0_pReal plasticState(ph)%dotState(:,of) = 0.0_pReal
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -982,6 +982,7 @@ subroutine constitutive_phenopowerlaw_dotState(Tstar_v,ipc,ip,el)
enddo enddo
enddo twinFamiliesLoop2 enddo twinFamiliesLoop2
end subroutine constitutive_phenopowerlaw_dotState end subroutine constitutive_phenopowerlaw_dotState

View File

@ -1197,7 +1197,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity)
endforall endforall
enddo elementLooping7 enddo elementLooping7
!$END PARALLEL DO !$END PARALLEL DO
! --- CALCULATE STATE AND STRESS FOR PERTURBATION --- ! --- CALCULATE STATE AND STRESS FOR PERTURBATION ---
dPdF_perturbation1 = crystallite_dPdF0 ! initialize stiffness with known good values from last increment dPdF_perturbation1 = crystallite_dPdF0 ! initialize stiffness with known good values from last increment
@ -1208,7 +1207,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity)
do k = 1,3; do l = 1,3 ! ...alter individual components do k = 1,3; do l = 1,3 ! ...alter individual components
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) & if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) &
write(6,'(a,2(1x,i1),1x,a,/)') '<< CRYST >> [[[[[[ Stiffness perturbation',k,l,']]]]]]' write(6,'(a,2(1x,i1),1x,a,/)') '<< CRYST >> [[[[[[ Stiffness perturbation',k,l,']]]]]]'
! --- INITIALIZE UNPERTURBED STATE --- ! --- INITIALIZE UNPERTURBED STATE ---
select case(numerics_integrator(numerics_integrationMode)) select case(numerics_integrator(numerics_integrationMode))
@ -3040,6 +3038,7 @@ subroutine crystallite_integrateStateFPI()
g,i,e) g,i,e)
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMP ENDDO !$OMP ENDDO
!$OMP DO PRIVATE(p,c) !$OMP DO PRIVATE(p,c)
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
@ -3049,7 +3048,8 @@ subroutine crystallite_integrateStateFPI()
c = mappingConstitutive(1,g,i,e) c = mappingConstitutive(1,g,i,e)
if ( any(plasticState(p)%dotState(:,c) /= plasticState(p)%dotState(:,c)) .or.& if ( any(plasticState(p)%dotState(:,c) /= plasticState(p)%dotState(:,c)) .or.&
any(damageState(p)%dotState(:,c) /= damageState(p)%dotState(:,c)) .or.& any(damageState(p)%dotState(:,c) /= damageState(p)%dotState(:,c)) .or.&
any(thermalState(p)%dotState(:,c) /= thermalState(p)%dotState(:,c))) then ! NaN occured in dotState any(thermalState(p)%dotState(:,c) /= thermalState(p)%dotState(:,c))) then !NaN occured in dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken is a non-local... if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken is a non-local...
!$OMP CRITICAL (checkTodo) !$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals done (and broken) crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals done (and broken)
@ -3062,7 +3062,6 @@ subroutine crystallite_integrateStateFPI()
enddo; enddo; enddo enddo; enddo; enddo
!$OMP ENDDO !$OMP ENDDO
! --- UPDATE STATE --- ! --- UPDATE STATE ---
!$OMP DO PRIVATE(mySizeDotState,mySizePlasticDotState,mySizeDamageDotState,mySizeThermalDotState,p,c) !$OMP DO PRIVATE(mySizeDotState,mySizePlasticDotState,mySizeDamageDotState,mySizeThermalDotState,p,c)
@ -3088,7 +3087,6 @@ subroutine crystallite_integrateStateFPI()
!$OMP END PARALLEL !$OMP END PARALLEL
! --+>> STATE LOOP <<+-- ! --+>> STATE LOOP <<+--
NiterationState = 0_pInt NiterationState = 0_pInt
@ -3105,7 +3103,7 @@ subroutine crystallite_integrateStateFPI()
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) & if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) &
call constitutive_microstructure(crystallite_temperature(i,e), crystallite_Fe(1:3,1:3,g,i,e), & call constitutive_microstructure(crystallite_temperature(i,e), crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states
call constitutive_damage_microstructure(crystallite_Tstar_v(1:6,g,i,e), & call constitutive_damage_microstructure(crystallite_Tstar_v(1:6,g,i,e), & !ToDo: should this be independent of todo/converged?
crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fe(1:3,1:3,g,i,e), &
g,i,e) ! update dependent state variables to be consistent with basic states g,i,e) ! update dependent state variables to be consistent with basic states
call constitutive_thermal_microstructure(crystallite_Tstar_v(1:6,g,i,e), & call constitutive_thermal_microstructure(crystallite_Tstar_v(1:6,g,i,e), &
@ -3122,7 +3120,6 @@ subroutine crystallite_integrateStateFPI()
enddo; enddo; enddo enddo; enddo; enddo
!$OMP ENDDO !$OMP ENDDO
! --- STRESS INTEGRATION --- ! --- STRESS INTEGRATION ---
!$OMP DO !$OMP DO
@ -3140,7 +3137,6 @@ subroutine crystallite_integrateStateFPI()
enddo; enddo; enddo enddo; enddo; enddo
!$OMP ENDDO !$OMP ENDDO
!$OMP SINGLE !$OMP SINGLE
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) & if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) &
@ -3184,11 +3180,9 @@ subroutine crystallite_integrateStateFPI()
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMP ENDDO !$OMP ENDDO
! --- UPDATE STATE --- ! --- UPDATE STATE ---
!$OMP DO PRIVATE(dot_prod12,dot_prod22, & !$OMP DO PRIVATE(dot_prod12,dot_prod22, &
@ -3206,10 +3200,8 @@ subroutine crystallite_integrateStateFPI()
mySizeThermalDotState = thermalState(p)%sizeDotState mySizeThermalDotState = thermalState(p)%sizeDotState
dot_prod12 = dot_product(plasticState(p)%dotState(:,c)- plasticState(p)%previousDotState(:,c), & dot_prod12 = dot_product(plasticState(p)%dotState(:,c)- plasticState(p)%previousDotState(:,c), &
plasticState(p)%previousDotState(:,c) - plasticState(p)%previousDotState2(:,c)) plasticState(p)%previousDotState(:,c) - plasticState(p)%previousDotState2(:,c))
dot_prod22 = dot_product(plasticState(p)%previousDotState(:,c) - plasticState(p)%previousDotState2(:,c), & dot_prod22 = dot_product(plasticState(p)%previousDotState(:,c) - plasticState(p)%previousDotState2(:,c), &
plasticState(p)%previousDotState(:,c) - plasticState(p)%previousDotState2(:,c)) plasticState(p)%previousDotState(:,c) - plasticState(p)%previousDotState2(:,c))
if ( dot_prod22 > 0.0_pReal & if ( dot_prod22 > 0.0_pReal &
.and. ( dot_prod12 < 0.0_pReal & .and. ( dot_prod12 < 0.0_pReal &
.or. dot_product(plasticState(p)%dotState(:,c), & .or. dot_product(plasticState(p)%dotState(:,c), &
@ -3218,13 +3210,10 @@ subroutine crystallite_integrateStateFPI()
else else
stateDamper = 1.0_pReal stateDamper = 1.0_pReal
endif endif
dot_prod12 = dot_product(damageState(p)%dotState(:,c)- damageState(p)%previousDotState(:,c), & dot_prod12 = dot_product(damageState(p)%dotState(:,c)- damageState(p)%previousDotState(:,c), &
damageState(p)%previousDotState(:,c) - damageState(p)%previousDotState2(:,c)) damageState(p)%previousDotState(:,c) - damageState(p)%previousDotState2(:,c))
dot_prod22 = dot_product(damageState(p)%previousDotState(:,c) - damageState(p)%previousDotState2(:,c), & dot_prod22 = dot_product(damageState(p)%previousDotState(:,c) - damageState(p)%previousDotState2(:,c), &
damageState(p)%previousDotState(:,c) - damageState(p)%previousDotState2(:,c)) damageState(p)%previousDotState(:,c) - damageState(p)%previousDotState2(:,c))
if ( dot_prod22 > 0.0_pReal & if ( dot_prod22 > 0.0_pReal &
.and. ( dot_prod12 < 0.0_pReal & .and. ( dot_prod12 < 0.0_pReal &
.or. dot_product(damageState(p)%dotState(:,c), & .or. dot_product(damageState(p)%dotState(:,c), &
@ -3233,13 +3222,10 @@ subroutine crystallite_integrateStateFPI()
else else
damageStateDamper = 1.0_pReal damageStateDamper = 1.0_pReal
endif endif
dot_prod12 = dot_product(thermalState(p)%dotState(:,c)- thermalState(p)%previousDotState(:,c), & dot_prod12 = dot_product(thermalState(p)%dotState(:,c)- thermalState(p)%previousDotState(:,c), &
thermalState(p)%previousDotState(:,c) - thermalState(p)%previousDotState2(:,c)) thermalState(p)%previousDotState(:,c) - thermalState(p)%previousDotState2(:,c))
dot_prod22 = dot_product(thermalState(p)%previousDotState(:,c) - thermalState(p)%previousDotState2(:,c), & dot_prod22 = dot_product(thermalState(p)%previousDotState(:,c) - thermalState(p)%previousDotState2(:,c), &
thermalState(p)%previousDotState(:,c) - thermalState(p)%previousDotState2(:,c)) thermalState(p)%previousDotState(:,c) - thermalState(p)%previousDotState2(:,c))
if ( dot_prod22 > 0.0_pReal & if ( dot_prod22 > 0.0_pReal &
.and. ( dot_prod12 < 0.0_pReal & .and. ( dot_prod12 < 0.0_pReal &
.or. dot_product(thermalState(p)%dotState(:,c), & .or. dot_product(thermalState(p)%dotState(:,c), &
@ -3267,9 +3253,7 @@ subroutine crystallite_integrateStateFPI()
- (thermalState(p)%dotState(1:mySizeThermalDotState,c) * thermalStateDamper & - (thermalState(p)%dotState(1:mySizeThermalDotState,c) * thermalStateDamper &
+ thermalState(p)%previousDotState(1:mySizeThermalDotState,c) & + thermalState(p)%previousDotState(1:mySizeThermalDotState,c) &
* (1.0_pReal - thermalStateDamper)) * crystallite_subdt(g,i,e) * (1.0_pReal - thermalStateDamper)) * crystallite_subdt(g,i,e)
! --- correct state with residuum --- ! --- correct state with residuum ---
tempState(1:mySizePlasticDotState) = plasticState(p)%state(1:mySizePlasticDotState,c) & tempState(1:mySizePlasticDotState) = plasticState(p)%state(1:mySizePlasticDotState,c) &
- stateResiduum(1:mySizePlasticDotState) ! need to copy to local variable, since we cant flush a pointer in openmp - stateResiduum(1:mySizePlasticDotState) ! need to copy to local variable, since we cant flush a pointer in openmp
tempDamageState(1:mySizeDamageDotState) = damageState(p)%state(1:mySizeDamageDotState,c) & tempDamageState(1:mySizeDamageDotState) = damageState(p)%state(1:mySizeDamageDotState,c) &
@ -3282,11 +3266,10 @@ subroutine crystallite_integrateStateFPI()
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,i8,1x,i2,1x,i3,/)') '<< CRYST >> update state at el ip g ',e,i,g write(6,'(a,i8,1x,i2,1x,i3,/)') '<< CRYST >> update state at el ip g ',e,i,g
write(6,'(a,f6.1,/)') '<< CRYST >> statedamper ',statedamper write(6,'(a,f6.1,/)') '<< CRYST >> statedamper ',statedamper
write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> state residuum',stateResiduum(1:mySizeDotState) write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> state residuum',stateResiduum(1:mySizePlasticDotState)
write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> new state',tempState(1:mySizeDotState) write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> new state',tempState(1:mySizePlasticDotState)
endif endif
#endif #endif
! --- store corrected dotState --- (cannot do this before state update, because not sure how to flush pointers in openmp) ! --- store corrected dotState --- (cannot do this before state update, because not sure how to flush pointers in openmp)
plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * stateDamper & plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * stateDamper &
@ -3317,13 +3300,12 @@ subroutine crystallite_integrateStateFPI()
!$OMP END CRITICAL (distributionState) !$OMP END CRITICAL (distributionState)
endif endif
endif endif
plasticState(p)%state(1:mySizePlasticDotState,c) = tempState(1:mySizePlasticDotState) ! copy local backup to global pointer plasticState(p)%dotState(1:mySizePlasticDotState,c) = tempState(1:mySizePlasticDotState) ! copy local backup to global pointer
damageState(p)%state (1:mySizeDamageDotState, c) = tempDamageState(1:mySizeDamageDotState) ! copy local backup to global pointer damageState(p)%dotState (1:mySizeDamageDotState, c) = tempDamageState(1:mySizeDamageDotState) ! copy local backup to global pointer
thermalState(p)%state(1:mySizeThermalDotState,c) = tempThermalState(1:mySizeThermalDotState) ! copy local backup to global pointer thermalState(p)%dotState(1:mySizeThermalDotState,c) = tempThermalState(1:mySizeThermalDotState) ! copy local backup to global pointer
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMP ENDDO !$OMP ENDDO
! --- STATE JUMP --- ! --- STATE JUMP ---
!$OMP DO !$OMP DO
@ -3343,7 +3325,6 @@ subroutine crystallite_integrateStateFPI()
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMP ENDDO !$OMP ENDDO
!$OMP END PARALLEL !$OMP END PARALLEL
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) & if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) &
@ -3364,7 +3345,6 @@ subroutine crystallite_integrateStateFPI()
write(6,'(a,i8,a,i2,/)') '<< CRYST >> ', count(crystallite_todo(:,:,:)), & write(6,'(a,i8,a,i2,/)') '<< CRYST >> ', count(crystallite_todo(:,:,:)), &
' grains todo after state integration no. ', NiterationState ' grains todo after state integration no. ', NiterationState
endif endif
! --- CHECK IF DONE WITH INTEGRATION --- ! --- CHECK IF DONE WITH INTEGRATION ---
doneWithIntegration = .true. doneWithIntegration = .true.