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

Fixed a physics issue: sense change with respect to difference between resolved stress and backstress

This commit is contained in:
Zhuowen Zhao 2018-02-14 22:13:10 -05:00
parent 7f487bb77b
commit f26fd1d1dc
7 changed files with 18 additions and 179 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

0
lib/damask/util.py Executable file → Normal file
View File

51
src/DAMASK_spectral.f90 Executable file → Normal file
View File

@ -522,13 +522,9 @@ program DAMASK_spectral
real(loadCases(currentLoadCase)%incs ,pReal)))
endif
endif
<<<<<<< HEAD
timeinc = timeinc / 2.0_pReal**real(cutBackLevel,pReal) ! depending on cut back level, decrease time step
! QUESTION: what happens to inc-counter when cutbacklevel is not zero? not clear where half an inc gets incremented..?
=======
timeinc = timeinc / real(subStepFactor,pReal)**real(cutBackLevel,pReal) ! depending on cut back level, decrease time step
>>>>>>> spectralSolver-cutbackfix
skipping: if (totalIncsCounter < restartInc) then ! not yet at restart inc?
time = time + timeinc ! just advance time, skip already performed calculation
guess = .true. ! QUESTION:why forced guessing instead of inheriting loadcase preference
@ -633,38 +629,7 @@ program DAMASK_spectral
stagIter = stagIter + 1_pInt
stagIterate = stagIter < stagItMax &
.and. all(solres(:)%converged) &
<<<<<<< HEAD
.and. .not. all(solres(:)%stagConverged)
enddo
!--------------------------------------------------------------------------------------------------
! check solution
cutBack = .False.
if (solres(1)%termIll &
.or. .not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found
! QUESTION: why termIll checked only for first field? only one that can be mechanic?
if (cutBackLevel < maxCutBack) then ! further cutbacking tolerated?
write(6,'(/,a)') ' cutting back '
cutBack = .true.
stepFraction = (stepFraction - 1_pInt) * subStepFactor ! adjust to new denominator
cutBackLevel = cutBackLevel + 1_pInt
time = time - timeinc ! rewind time
timeinc = timeinc/2.0_pReal
elseif (continueCalculation == 1_pInt .and. .not. solres(1)%termIll) then
guess = .true. ! accept non converged BVP solution
else ! material point model cannot find a solution
call IO_warning(850_pInt)
call MPI_file_close(resUnit,ierr)
close(statUnit)
call quit(-1_pInt*(lastRestartWritten+1_pInt)) ! quit and provide information about last restart inc written
endif
else
guess = .true. ! start guessing after first converged (sub)inc
endif
if (.not. cutBack) then
=======
.and. .not. all(solres(:)%stagConverged) ! stationary with respect to staggered iteration
enddo
@ -676,7 +641,6 @@ program DAMASK_spectral
timeIncOld = timeinc
cutBack = .false.
guess = .true. ! start guessing after first converged (sub)inc
>>>>>>> spectralSolver-cutbackfix
if (worldrank == 0) then
write(statUnit,*) totalIncsCounter, time, cutBackLevel, &
solres%converged, solres%iterationsNeeded
@ -695,10 +659,7 @@ program DAMASK_spectral
close(statUnit)
call quit(-1_pInt*(lastRestartWritten+1_pInt)) ! quit and provide information about last restart inc written
endif
<<<<<<< HEAD
=======
>>>>>>> spectralSolver-cutbackfix
enddo subStepLooping
cutBackLevel = max(0_pInt, cutBackLevel - 1_pInt) ! try half number of subincs next inc
@ -719,11 +680,8 @@ program DAMASK_spectral
flush(6)
call materialpoint_postResults()
call MPI_file_seek (resUnit,fileOffset,MPI_SEEK_SET,ierr)
<<<<<<< HEAD
if (ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_file_seek')
=======
if (ierr /= 0_pInt) call IO_error(894_pInt, ext_msg='MPI_file_seek')
>>>>>>> spectralSolver-cutbackfix
do i=1, size(materialpoint_results,3)/(maxByteOut/(materialpoint_sizeResults*pReal))+1 ! slice the output of my process in chunks not exceeding the limit for one output
outputIndex=int([(i-1_pInt)*((maxRealOut)/materialpoint_sizeResults)+1_pInt, &
min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt)
@ -740,9 +698,7 @@ program DAMASK_spectral
restartWrite = .true. ! set restart parameter for FEsolving
lastRestartWritten = inc ! QUESTION: first call to CPFEM_general will write?
endif
<<<<<<< HEAD
endif skipping
=======
else forwarding
time = time + timeinc
guess = .true.
@ -812,7 +768,6 @@ program DAMASK_spectral
call quit(0_pInt)
endif
endif
>>>>>>> development
enddo incLooping
enddo loadCaseLooping

0
src/crystallite.f90 Executable file → Normal file
View File

0
src/math.f90 Executable file → Normal file
View File

18
src/plastic_kinematichardening.f90
View File

@ -1,5 +1,5 @@
!--------------------------------------------------------------------------------------------------
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Michigan State University
!> @author Zhuowen Zhao, Michigan State University
!> @brief Introducing Voce-type kinematic hardening rule into crystal plasticity
!! formulation using a power law fitting
@ -611,13 +611,13 @@ subroutine plastic_kinehardening_shearRates(gdot_pos,gdot_neg,tau_pos,tau_neg, &
enddo slipFamilies
gdot_pos = 0.5_pReal * param(instance)%gdot0 * &
(abs(tau_pos-state(instance)%sense(:,of)*state(instance)%crss_back(:,of))/ &
(abs(tau_pos-state(instance)%crss_back(:,of))/ &
state(instance)%crss(:,of))**param(instance)%n_slip &
*sign(1.0_pReal,tau_pos)
*sign(1.0_pReal,tau_pos-state(instance)%crss_back(:,of))
gdot_neg = 0.5_pReal * param(instance)%gdot0 * &
(abs(tau_neg-state(instance)%sense(:,of)*state(instance)%crss_back(:,of))/ &
(abs(tau_neg-state(instance)%crss_back(:,of))/ &
state(instance)%crss(:,of))**param(instance)%n_slip &
*sign(1.0_pReal,tau_neg)
*sign(1.0_pReal,tau_neg-state(instance)%crss_back(:,of))
! gdot_pos = 0.5_pReal * param(instance)%gdot0 * &
! exp(-param(instance)%F0/(1.38e-23*298.15)* &
@ -765,7 +765,8 @@ end subroutine plastic_kinehardening_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
subroutine plastic_kinehardening_deltaState(Tstar_v,ipc,ip,el)
use prec, only: &
dNeq
dNeq, &
dEq0
use debug, only: &
debug_level, &
debug_constitutive, &
@ -804,8 +805,9 @@ subroutine plastic_kinehardening_deltaState(Tstar_v,ipc,ip,el)
call plastic_kinehardening_shearRates(gdot_pos,gdot_neg,tau_pos,tau_neg, &
Tstar_v,ph,instance,of)
sense = sign(1.0_pReal,gdot_pos+gdot_neg) ! current sense of shear direction
sense = merge(state(instance)%sense(:,of), & ! keep existing...
sign(1.0_pReal,gdot_pos+gdot_neg), & ! ...or have a defined
dEq0(gdot_pos+gdot_neg,1e-10_pReal)) ! current sense of shear direction
#ifdef DEBUG
if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0_pInt &

114
src/spectral_mech_Basic.f90
View File

@ -280,13 +280,8 @@ type(tSolutionState) function basicPETSc_solution(incInfoIn,timeinc,timeinc_old,
basicPETSC_solution%iterationsNeeded = totalIter
basicPETSc_solution%termIll = terminallyIll
terminallyIll = .false.
<<<<<<< HEAD
if (reason == -4) call IO_error(893_pInt)
BasicPETSc_solution%converged = reason > 0
basicPETSC_solution%iterationsNeeded = totalIter
=======
if (reason == -4) call IO_error(893_pInt) ! MPI error
>>>>>>> spectralSolver-cutbackfix
end function BasicPETSc_solution
@ -343,10 +338,6 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr)
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment
<<<<<<< HEAD
newIteration: if (totalIter <= PETScIter) then
=======
>>>>>>> spectralSolver-cutbackfix
!--------------------------------------------------------------------------------------------------
! begin of new iteration
newIteration: if (totalIter <= PETScIter) then
@ -449,106 +440,6 @@ end subroutine BasicPETSc_converged
!> possibly writing restart information, triggering of state increment in DAMASK, and updating of IPcoordinates
!--------------------------------------------------------------------------------------------------
subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC)
<<<<<<< HEAD
use math, only: &
math_mul33x33 ,&
math_rotate_backward33
use numerics, only: &
worldrank
use mesh, only: &
grid, &
grid3
use spectral_utilities, only: &
Utilities_calculateRate, &
Utilities_forwardField, &
Utilities_updateIPcoords, &
tBoundaryCondition, &
cutBack
use IO, only: &
IO_write_JobRealFile
use FEsolving, only: &
restartWrite
implicit none
real(pReal), intent(in) :: &
timeinc_old, &
timeinc, &
loadCaseTime !< remaining time of current load case
type(tBoundaryCondition), intent(in) :: &
stress_BC, &
deformation_BC
real(pReal), dimension(3,3), intent(in) :: rotation_BC
logical, intent(in) :: &
guess
PetscErrorCode :: ierr
PetscScalar, pointer :: F(:,:,:,:)
character(len=1024) :: rankStr
call DMDAVecGetArrayF90(da,solution_vec,F,ierr) ! get F from PETSc data structure
!--------------------------------------------------------------------------------------------------
! restart information for spectral solver
if (restartWrite) then ! QUESTION: where is this logical properly set?
write(6,'(/,a)') ' writing converged results for restart'
flush(6)
write(rankStr,'(a1,i0)')'_',worldrank
call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file
write (777,rec=1) F
close (777)
call IO_write_jobRealFile(777,'F_lastInc'//trim(rankStr),size(F_lastInc)) ! writing F_lastInc field to file
write (777,rec=1) F_lastInc
close (777)
if (worldrank == 0_pInt) then
call IO_write_jobRealFile(777,'F_aimDot',size(F_aimDot))
write (777,rec=1) F_aimDot
close(777)
call IO_write_jobRealFile(777,'C_volAvg',size(C_volAvg))
write (777,rec=1) C_volAvg
close(777)
call IO_write_jobRealFile(777,'C_volAvgLastInc',size(C_volAvgLastInc))
write (777,rec=1) C_volAvgLastInc
close(777)
endif
endif
call utilities_updateIPcoords(F) ! QUESTION: why do this even when cutback happened??
if (cutBack) then ! reset to former inc's values
F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) ! QUESTION: purpose of resetting this when updating in line 541?
F_aim = F_aim_lastInc
C_volAvg = C_volAvgLastInc
else
ForwardData = .true. ! QUESTION: who is resetting this?
C_volAvgLastInc = C_volAvg
!--------------------------------------------------------------------------------------------------
! calculate rate for aim
if (deformation_BC%myType=='l') then ! calculate f_aimDot from given L and current F
f_aimDot = deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim)
elseif(deformation_BC%myType=='fdot') then ! f_aimDot is prescribed
f_aimDot = deformation_BC%maskFloat * deformation_BC%values
elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed
f_aimDot = deformation_BC%maskFloat * (deformation_BC%values - F_aim)/loadCaseTime
endif
if (guess) f_aimDot = f_aimDot + stress_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old
F_aim_lastInc = F_aim
!--------------------------------------------------------------------------------------------------
! update coordinates and rate and forward last inc
call utilities_updateIPcoords(F)
Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3])) ! QUESTION: what do we need Fdot for and why is it not restored at cutback?
F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid3])
endif
F_aim = F_aim + f_aimDot * timeinc
!--------------------------------------------------------------------------------------------------
! update local deformation gradient
F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! ensure that it matches rotated F_aim
math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid3])
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
=======
use math, only: &
math_mul33x33 ,&
math_rotate_backward33
@ -660,8 +551,7 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation
F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average
math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid3])
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
>>>>>>> spectralSolver-cutbackfix
end subroutine BasicPETSc_forward
!--------------------------------------------------------------------------------------------------

14
src/spectral_utilities.f90 Executable file → Normal file
View File

@ -823,11 +823,8 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
.or. (m/=n .and. abs(sTimesC(m,n)) > 1.0e-12_pReal) ! off-diagonal elements of S*C should be 0
enddo
enddo
<<<<<<< HEAD
if(debugGeneral .or. errmatinv) then
=======
if (debugGeneral .or. errmatinv) then
>>>>>>> spectralSolver-cutbackfix
write(formatString, '(i2)') size_reduced
formatString = '(/,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))'
write(6,trim(formatString),advance='no') ' C * S (load) ', &
@ -841,18 +838,13 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
else
temp99_real = 0.0_pReal
endif
<<<<<<< HEAD
if(debugGeneral) &
write(6,'(/,a,/,9(9(2x,f10.5,1x)/),/)',advance='no') ' Masked Compliance (load) / GPa =', &
transpose(temp99_Real*1.e9_pReal)
flush(6)
=======
if(debugGeneral) then
write(6,'(/,a,/,9(9(2x,f10.5,1x)/),/)',advance='no') &
' Masked Compliance (load) / GPa =', transpose(temp99_Real*1.e-9_pReal)
flush(6)
endif
>>>>>>> spectralSolver-cutbackfix
utilities_maskedCompliance = math_Plain99to3333(temp99_Real)
end function utilities_maskedCompliance