From f26fd1d1dc4f5567fcc9269c3598232710c0b649 Mon Sep 17 00:00:00 2001 From: Zhuowen Zhao Date: Wed, 14 Feb 2018 22:13:10 -0500 Subject: [PATCH] Fixed a physics issue: sense change with respect to difference between resolved stress and backstress --- lib/damask/util.py | 0 src/DAMASK_spectral.f90 | 51 +------------ src/crystallite.f90 | 0 src/math.f90 | 0 src/plastic_kinematichardening.f90 | 18 +++-- src/spectral_mech_Basic.f90 | 114 +---------------------------- src/spectral_utilities.f90 | 14 +--- 7 files changed, 18 insertions(+), 179 deletions(-) mode change 100755 => 100644 lib/damask/util.py mode change 100755 => 100644 src/DAMASK_spectral.f90 mode change 100755 => 100644 src/crystallite.f90 mode change 100755 => 100644 src/math.f90 mode change 100755 => 100644 src/spectral_utilities.f90 diff --git a/lib/damask/util.py b/lib/damask/util.py old mode 100755 new mode 100644 diff --git a/src/DAMASK_spectral.f90 b/src/DAMASK_spectral.f90 old mode 100755 new mode 100644 index 3b5fedf18..5cda40249 --- a/src/DAMASK_spectral.f90 +++ b/src/DAMASK_spectral.f90 @@ -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 diff --git a/src/crystallite.f90 b/src/crystallite.f90 old mode 100755 new mode 100644 diff --git a/src/math.f90 b/src/math.f90 old mode 100755 new mode 100644 diff --git a/src/plastic_kinematichardening.f90 b/src/plastic_kinematichardening.f90 index c457c1344..6d7812a74 100644 --- a/src/plastic_kinematichardening.f90 +++ b/src/plastic_kinematichardening.f90 @@ -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 & diff --git a/src/spectral_mech_Basic.f90 b/src/spectral_mech_Basic.f90 index 0ee977f6a..0b3409b52 100644 --- a/src/spectral_mech_Basic.f90 +++ b/src/spectral_mech_Basic.f90 @@ -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 !-------------------------------------------------------------------------------------------------- diff --git a/src/spectral_utilities.f90 b/src/spectral_utilities.f90 old mode 100755 new mode 100644 index 3bcc914a4..a66fa558e --- a/src/spectral_utilities.f90 +++ b/src/spectral_utilities.f90 @@ -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