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

question marks on those files

This commit is contained in:
Zhuowen Zhao 2017-12-13 19:18:45 -05:00
parent 2b8baa2f01
commit dcf9e139d0
3 changed files with 114 additions and 77 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

78
src/DAMASK_spectral.f90
View File

@ -442,8 +442,9 @@ program DAMASK_spectral
if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_seek') if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_seek')
if (.not. appendToOutFile) then ! if not restarting, write 0th increment if (.not. appendToOutFile) then ! if not restarting, write 0th increment
write(6,'(1/,a)') ' ... writing initial configuration to file ........................'
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 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, & outputIndex = int([(i-1_pInt)*((maxRealOut)/materialpoint_sizeResults)+1_pInt, & ! QUESTION: why not starting i at 0 instead of murky 1?
min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt) min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt)
call MPI_file_write(resUnit, & call MPI_file_write(resUnit, &
reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)), & reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)), &
@ -453,7 +454,6 @@ program DAMASK_spectral
if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_write') if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_write')
enddo enddo
fileOffset = fileOffset + sum(outputSize) ! forward to current file position fileOffset = fileOffset + sum(outputSize) ! forward to current file position
write(6,'(1/,a)') ' ... writing initial configuration to file ........................'
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! loopping over loadcases ! loopping over loadcases
@ -487,19 +487,22 @@ program DAMASK_spectral
endif endif
endif endif
timeinc = timeinc / 2.0_pReal**real(cutBackLevel,pReal) ! depending on cut back level, decrease time step 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..?
forwarding: if (totalIncsCounter >= restartInc) then skipping: if (totalIncsCounter < restartInc) then ! not yet at restart inc?
stepFraction = 0_pInt time = time + timeinc ! just advance time, skip already performed calculation
guess = .true. ! QUESTION:why forced guessing instead of inheriting loadcase preference
else skipping
stepFraction = 0_pInt ! fraction scaled by stepFactor**cutLevel
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! loop over sub incs ! loop over sub step
subIncLooping: do while (stepFraction/subStepFactor**cutBackLevel <1_pInt) subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel)
time = time + timeinc ! forward time remainingLoadCaseTime = loadCases(currentLoadCase)%time+time0 - time
stepFraction = stepFraction + 1_pInt time = time + timeinc ! forward target time
remainingLoadCaseTime = time0 - time + loadCases(currentLoadCase)%time + timeInc stepFraction = stepFraction + 1_pInt ! count step
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report begin of new increment ! report begin of new step
write(6,'(/,a)') ' ###########################################################################' write(6,'(/,a)') ' ###########################################################################'
write(6,'(1x,a,es12.5'//& write(6,'(1x,a,es12.5'//&
',a,'//IO_intOut(inc)//',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//& ',a,'//IO_intOut(inc)//',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//&
@ -582,33 +585,33 @@ program DAMASK_spectral
solres(field) = spectral_damage_solution(timeinc,timeIncOld,remainingLoadCaseTime) solres(field) = spectral_damage_solution(timeinc,timeIncOld,remainingLoadCaseTime)
end select end select
if (.not. solres(field)%converged) exit ! no solution found if (.not. solres(field)%converged) exit ! no solution found
enddo enddo
stagIter = stagIter + 1_pInt stagIter = stagIter + 1_pInt
stagIterate = stagIter < stagItMax .and. & stagIterate = stagIter < stagItMax &
all(solres(:)%converged) .and. & .and. all(solres(:)%converged) &
.not. all(solres(:)%stagConverged) .and. .not. all(solres(:)%stagConverged)
enddo enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! check solution ! check solution
cutBack = .False. cutBack = .False.
if(solres(1)%termIll .or. .not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found
if (cutBackLevel < maxCutBack) then ! do cut back if (solres(1)%termIll &
write(6,'(/,a)') ' cut back detected' .or. .not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found
cutBack = .True. ! 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 stepFraction = (stepFraction - 1_pInt) * subStepFactor ! adjust to new denominator
cutBackLevel = cutBackLevel + 1_pInt cutBackLevel = cutBackLevel + 1_pInt
time = time - timeinc ! rewind time time = time - timeinc ! rewind time
timeinc = timeinc/2.0_pReal timeinc = timeinc/2.0_pReal
elseif (solres(1)%termIll) then ! material point model cannot find a solution, exit in any casy elseif (continueCalculation == 1_pInt .and. .not. solres(1)%termIll) then
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
elseif (continueCalculation == 1_pInt) then
guess = .true. ! accept non converged BVP solution guess = .true. ! accept non converged BVP solution
else ! default behavior, exit if spectral solver does not converge else ! material point model cannot find a solution
call IO_warning(850_pInt) call IO_warning(850_pInt)
call MPI_file_close(resUnit,ierr) call MPI_file_close(resUnit,ierr)
close(statUnit) close(statUnit)
@ -617,6 +620,7 @@ program DAMASK_spectral
else else
guess = .true. ! start guessing after first converged (sub)inc guess = .true. ! start guessing after first converged (sub)inc
endif endif
if (.not. cutBack) then if (.not. cutBack) then
if (worldrank == 0) then if (worldrank == 0) then
write(statUnit,*) totalIncsCounter, time, cutBackLevel, & write(statUnit,*) totalIncsCounter, time, cutBackLevel, &
@ -624,17 +628,20 @@ program DAMASK_spectral
flush(statUnit) flush(statUnit)
endif endif
endif endif
enddo subIncLooping enddo subStepLooping
cutBackLevel = max(0_pInt, cutBackLevel - 1_pInt) ! try half number of subincs next inc cutBackLevel = max(0_pInt, cutBackLevel - 1_pInt) ! try half number of subincs next inc
if(all(solres(:)%converged)) then ! report converged inc
if (all(solres(:)%converged)) then
convergedCounter = convergedCounter + 1_pInt convergedCounter = convergedCounter + 1_pInt
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report converged inc
' increment ', totalIncsCounter, ' converged' ' increment ', totalIncsCounter, ' converged'
else else
notConvergedCounter = notConvergedCounter + 1_pInt
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report non-converged inc write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report non-converged inc
' increment ', totalIncsCounter, ' NOT converged' ' increment ', totalIncsCounter, ' NOT converged'
notConvergedCounter = notConvergedCounter + 1_pInt
endif; flush(6) endif; flush(6)
if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency
if (worldrank == 0) & if (worldrank == 0) &
write(6,'(1/,a)') ' ... writing results to file ......................................' write(6,'(1/,a)') ' ... writing results to file ......................................'
@ -652,15 +659,12 @@ program DAMASK_spectral
enddo enddo
fileOffset = fileOffset + sum(outputSize) ! forward to current file position fileOffset = fileOffset + sum(outputSize) ! forward to current file position
endif endif
if( loadCases(currentLoadCase)%restartFrequency > 0_pInt .and. & ! at frequency of writing restart information set restart parameter for FEsolving if ( loadCases(currentLoadCase)%restartFrequency > 0_pInt & ! writing of restart info requested ...
mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0_pInt) then ! first call to CPFEM_general will write? .and. mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0_pInt) then ! ... and at frequency of writing restart information
restartWrite = .true. restartWrite = .true. ! set restart parameter for FEsolving
lastRestartWritten = inc lastRestartWritten = inc ! QUESTION: first call to CPFEM_general will write?
endif endif
else forwarding endif skipping
time = time + timeinc
guess = .true.
endif forwarding
enddo incLooping enddo incLooping
enddo loadCaseLooping enddo loadCaseLooping

56
src/crystallite.f90
View File

@ -986,7 +986,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
crystallite_todo(c,i,e) = crystallite_subStep(c,i,e) > subStepMinCryst ! still on track or already done (beyond repair) crystallite_todo(c,i,e) = crystallite_subStep(c,i,e) > subStepMinCryst ! still on track or already done (beyond repair)
!$OMP FLUSH(crystallite_todo) !$OMP FLUSH(crystallite_todo)
#ifdef DEBUG #ifdef DEBUG
if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) then 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)) then
if (crystallite_todo(c,i,e)) then if (crystallite_todo(c,i,e)) then
write(6,'(a,f12.8,a,i8,1x,i2,1x,i3,/)') '<< CRYST >> cutback step in crystallite_stressAndItsTangent & write(6,'(a,f12.8,a,i8,1x,i2,1x,i3,/)') '<< CRYST >> cutback step in crystallite_stressAndItsTangent &
&with new crystallite_subStep: ',& &with new crystallite_subStep: ',&
@ -1047,6 +1049,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
write(6,'(a,f8.5)') '<< CRYST >> min(subFrac) ',minval(crystallite_subFrac) write(6,'(a,f8.5)') '<< CRYST >> min(subFrac) ',minval(crystallite_subFrac)
write(6,'(a,f8.5,/)') '<< CRYST >> max(subFrac) ',maxval(crystallite_subFrac) write(6,'(a,f8.5,/)') '<< CRYST >> max(subFrac) ',maxval(crystallite_subFrac)
flush(6) flush(6)
if (iand(debug_level(debug_crystallite),debug_levelSelective) /= 0_pInt) then
write(6,'(/,a,f8.5,1x,a,1x,f8.5,1x,a)') '<< CRYST >> subFrac + subStep = ',&
crystallite_subFrac(debug_g,debug_i,debug_e),'+',crystallite_subStep(debug_g,debug_i,debug_e),'@selective'
flush(6)
endif
endif endif
! --- integrate --- requires fully defined state array (basic + dependent state) ! --- integrate --- requires fully defined state array (basic + dependent state)
@ -2752,7 +2759,7 @@ subroutine crystallite_integrateStateFPI()
enddo enddo
if (NaN) then ! NaN occured in any dotState if (NaN) then ! NaN occured in any dotState
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) & if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) &
write(6,*) '<< CRYST >> ',plasticState(p)%dotState(:,c) write(6,*) '<< CRYST >> dotstate ',plasticState(p)%dotState(:,c)
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)
@ -2831,9 +2838,6 @@ subroutine crystallite_integrateStateFPI()
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
!$OMP FLUSH(crystallite_todo) !$OMP FLUSH(crystallite_todo)
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e) crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e)
!$OMP FLUSH(crystallite_todo) !$OMP FLUSH(crystallite_todo)
if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! broken non-local... if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! broken non-local...
@ -2984,8 +2988,11 @@ 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 >> plasticstatedamper ',plasticStatedamper write(6,'(a,f6.1,/)') '<< CRYST >> plasticstatedamper ',plasticStatedamper
write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> plastic state residuum',plasticStateResiduum(1:mySizePlasticDotState) write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> plastic state residuum',&
abs(plasticStateResiduum(1:mySizePlasticDotState))
write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> abstol dotstate',plasticState(p)%aTolState(1:mySizePlasticDotState) write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> abstol dotstate',plasticState(p)%aTolState(1:mySizePlasticDotState)
write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> reltol dotstate',rTol_crystalliteState* &
abs(tempPlasticState(1:mySizePlasticDotState))
write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> new state',tempPlasticState(1:mySizePlasticDotState) write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> new state',tempPlasticState(1:mySizePlasticDotState)
endif endif
#endif #endif
@ -3419,7 +3426,7 @@ logical function crystallite_integrateStress(&
#ifdef DEBUG #ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) & if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) &
write(6,'(a,i3,a,i8,1x,a,i8,a,1x,i2,1x,i3,/)') '<< CRYST >> integrateStress reached loop limit',nStress, & write(6,'(a,i3,a,i8,1x,a,i8,a,1x,i2,1x,i3,/)') '<< CRYST >> integrateStress reached loop limit',nStress, &
' at el (elFE) ip ipc ', el,mesh_element(1,el),ip,ipc ' at el (elFE) ip ipc ', el,'(',mesh_element(1,el),')',ip,ipc
#endif #endif
return return
endif loopsExeced endif loopsExeced
@ -3428,7 +3435,8 @@ logical function crystallite_integrateStress(&
B = math_I3 - dt*Lpguess B = math_I3 - dt*Lpguess
Fe = math_mul33x33(math_mul33x33(A,B), invFi_new) ! current elastic deformation tensor Fe = math_mul33x33(math_mul33x33(A,B), invFi_new) ! current elastic deformation tensor
call constitutive_TandItsTangent(Tstar, dT_dFe3333, dT_dFi3333, Fe, Fi_new, ipc, ip, el) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative in unloaded configuration call constitutive_TandItsTangent(Tstar, dT_dFe3333, dT_dFi3333, &
Fe, Fi_new, ipc, ip, el) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative in unloaded configuration
Tstar_v = math_Mandel33to6(Tstar) Tstar_v = math_Mandel33to6(Tstar)
!* calculate plastic velocity gradient and its tangent from constitutive law !* calculate plastic velocity gradient and its tangent from constitutive law
@ -3436,6 +3444,17 @@ logical function crystallite_integrateStress(&
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) & if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) &
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,i3,/)') '<< CRYST >> stress iteration ', NiterationStressLp
write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST >> Lpguess', math_transpose33(Lpguess)
write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST >> Fi', math_transpose33(Fi_new)
write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST >> Fe', math_transpose33(Fe)
write(6,'(a,/,6(e20.10,1x))') '<< CRYST >> Tstar', Tstar_v
endif
#endif
call constitutive_LpAndItsTangent(Lp_constitutive, dLp_dT3333, dLp_dFi3333, & call constitutive_LpAndItsTangent(Lp_constitutive, dLp_dT3333, dLp_dFi3333, &
Tstar_v, Fi_new, ipc, ip, el) Tstar_v, Fi_new, ipc, ip, el)
@ -3453,9 +3472,7 @@ logical function crystallite_integrateStress(&
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
.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,i3,/)') '<< CRYST >> stress iteration ', NiterationStressLp write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST >> Lp_constitutive', math_transpose33(Lp_constitutive)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lp_constitutive', math_transpose33(Lp_constitutive)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lpguess', math_transpose33(Lpguess)
endif endif
#endif #endif
@ -3485,6 +3502,13 @@ logical function crystallite_integrateStress(&
else ! not converged and residuum not improved... else ! not converged and residuum not improved...
steplengthLp = subStepSizeLp * steplengthLp ! ...try with smaller step length in same direction steplengthLp = subStepSizeLp * steplengthLp ! ...try with smaller step length in same direction
Lpguess = Lpguess_old + steplengthLp * deltaLp Lpguess = Lpguess_old + steplengthLp * deltaLp
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,1x,f7.4)') '<< CRYST >> linear search for Lpguess with step', steplengthLp
endif
#endif
cycle LpLoop cycle LpLoop
endif endif
@ -3498,6 +3522,16 @@ logical function crystallite_integrateStress(&
dFe_dLp3333 = - dt * dFe_dLp3333 dFe_dLp3333 = - dt * dFe_dLp3333
dRLp_dLp = math_identity2nd(9_pInt) & dRLp_dLp = math_identity2nd(9_pInt) &
- math_Plain3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dT3333,dT_dFe3333),dFe_dLp3333)) - math_Plain3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dT3333,dT_dFe3333),dFe_dLp3333))
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,/,9(12x,9(e12.4,1x)/))') '<< CRYST >> dLp_dT', math_Plain3333to99(dLp_dT3333)
write(6,'(a,1x,e20.10)') '<< CRYST >> dLp_dT norm', norm2(math_Plain3333to99(dLp_dT3333))
write(6,'(a,/,9(12x,9(e12.4,1x)/))') '<< CRYST >> dRLp_dLp', dRLp_dLp - math_identity2nd(9_pInt)
write(6,'(a,1x,e20.10)') '<< CRYST >> dRLp_dLp norm', norm2(dRLp_dLp - math_identity2nd(9_pInt))
endif
#endif
dRLp_dLp2 = dRLp_dLp ! will be overwritten in first call to LAPACK routine dRLp_dLp2 = dRLp_dLp ! will be overwritten in first call to LAPACK routine
work = math_plain33to9(residuumLp) work = math_plain33to9(residuumLp)
call dgesv(9,1,dRLp_dLp2,9,ipiv,work,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp call dgesv(9,1,dRLp_dLp2,9,ipiv,work,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp

23
src/spectral_mech_Basic.f90
View File

@ -196,8 +196,9 @@ subroutine basicPETSc_init
.false., & .false., &
math_I3) math_I3)
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! write data back to PETSc call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! write data back to PETSc
! QUESTION: why not writing back right after reading (l.189)?
restartRead: if (restartInc > 1_pInt) then restartRead: if (restartInc > 1_pInt) then ! QUESTION: are those values not calc'ed by constitutiveResponse? why reading from file?
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) & if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) &
write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'reading more values of increment', restartInc - 1_pInt, 'from file' 'reading more values of increment', restartInc - 1_pInt, 'from file'
@ -220,8 +221,7 @@ end subroutine basicPETSc_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief solution for the Basic PETSC scheme with internal iterations !> @brief solution for the Basic PETSC scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
type(tSolutionState) function & type(tSolutionState) function basicPETSc_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC)
basicPETSc_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC)
use IO, only: & use IO, only: &
IO_error IO_error
use numerics, only: & use numerics, only: &
@ -283,9 +283,8 @@ type(tSolutionState) function &
CHKERRQ(ierr) CHKERRQ(ierr)
basicPETSc_solution%termIll = terminallyIll basicPETSc_solution%termIll = terminallyIll
terminallyIll = .false. terminallyIll = .false.
BasicPETSc_solution%converged =.true.
if (reason == -4) call IO_error(893_pInt) if (reason == -4) call IO_error(893_pInt)
if (reason < 1) basicPETSC_solution%converged = .false. BasicPETSc_solution%converged = reason > 0
basicPETSC_solution%iterationsNeeded = totalIter basicPETSC_solution%iterationsNeeded = totalIter
end function BasicPETSc_solution end function BasicPETSc_solution
@ -480,10 +479,10 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation
character(len=1024) :: rankStr character(len=1024) :: rankStr
call DMDAVecGetArrayF90(da,solution_vec,F,ierr) call DMDAVecGetArrayF90(da,solution_vec,F,ierr) ! get F from PETSc data structure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! restart information for spectral solver ! restart information for spectral solver
if (restartWrite) then if (restartWrite) then ! QUESTION: where is this logical properly set?
write(6,'(/,a)') ' writing converged results for restart' write(6,'(/,a)') ' writing converged results for restart'
flush(6) flush(6)
write(rankStr,'(a1,i0)')'_',worldrank write(rankStr,'(a1,i0)')'_',worldrank
@ -506,14 +505,14 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation
endif endif
endif endif
call utilities_updateIPcoords(F) call utilities_updateIPcoords(F) ! QUESTION: why do this even when cutback happened??
if (cutBack) then 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 F_aim = F_aim_lastInc
F = reshape(F_lastInc, [9,grid(1),grid(2),grid3])
C_volAvg = C_volAvgLastInc C_volAvg = C_volAvgLastInc
else else
ForwardData = .True. ForwardData = .true. ! QUESTION: who is resetting this?
C_volAvgLastInc = C_volAvg C_volAvgLastInc = C_volAvg
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate rate for aim ! calculate rate for aim
@ -531,7 +530,7 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation
! update coordinates and rate and forward last inc ! update coordinates and rate and forward last inc
call utilities_updateIPcoords(F) call utilities_updateIPcoords(F)
Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3])) 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]) F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid3])
endif endif