From dcf9e139d0d603d3245b9b95c518c6c703cc8d56 Mon Sep 17 00:00:00 2001 From: Zhuowen Zhao Date: Wed, 13 Dec 2017 19:18:45 -0500 Subject: [PATCH] question marks on those files --- src/DAMASK_spectral.f90 | 82 +++++++++++++++++++------------------ src/crystallite.f90 | 70 +++++++++++++++++++++++-------- src/spectral_mech_Basic.f90 | 39 +++++++++--------- 3 files changed, 114 insertions(+), 77 deletions(-) diff --git a/src/DAMASK_spectral.f90 b/src/DAMASK_spectral.f90 index f32bfb7b3..39eb77bc7 100644 --- a/src/DAMASK_spectral.f90 +++ b/src/DAMASK_spectral.f90 @@ -442,8 +442,9 @@ program DAMASK_spectral 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 + 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 - 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) call MPI_file_write(resUnit, & 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') enddo fileOffset = fileOffset + sum(outputSize) ! forward to current file position - write(6,'(1/,a)') ' ... writing initial configuration to file ........................' endif !-------------------------------------------------------------------------------------------------- ! loopping over loadcases @@ -487,19 +487,22 @@ program DAMASK_spectral endif endif timeinc = timeinc / 2.0_pReal**real(cutBackLevel,pReal) ! depending on cut back level, decrease time step - - forwarding: if (totalIncsCounter >= restartInc) then - stepFraction = 0_pInt + ! QUESTION: what happens to inc-counter when cutbacklevel is not zero? not clear where half an inc gets incremented..? + 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 + else skipping + stepFraction = 0_pInt ! fraction scaled by stepFactor**cutLevel !-------------------------------------------------------------------------------------------------- -! loop over sub incs - subIncLooping: do while (stepFraction/subStepFactor**cutBackLevel <1_pInt) - time = time + timeinc ! forward time - stepFraction = stepFraction + 1_pInt - remainingLoadCaseTime = time0 - time + loadCases(currentLoadCase)%time + timeInc +! loop over sub step + subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel) + remainingLoadCaseTime = loadCases(currentLoadCase)%time+time0 - time + time = time + timeinc ! forward target time + stepFraction = stepFraction + 1_pInt ! count step !-------------------------------------------------------------------------------------------------- -! report begin of new increment +! report begin of new step write(6,'(/,a)') ' ###########################################################################' write(6,'(1x,a,es12.5'//& ',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) end select + if (.not. solres(field)%converged) exit ! no solution found + enddo stagIter = stagIter + 1_pInt - stagIterate = stagIter < stagItMax .and. & - all(solres(:)%converged) .and. & - .not. all(solres(:)%stagConverged) + stagIterate = stagIter < stagItMax & + .and. all(solres(:)%converged) & + .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 - if (cutBackLevel < maxCutBack) then ! do cut back - write(6,'(/,a)') ' cut back detected' - cutBack = .True. + + 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 (solres(1)%termIll) then ! material point model cannot find a solution, exit in any casy - 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 + elseif (continueCalculation == 1_pInt .and. .not. solres(1)%termIll) then 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 MPI_file_close(resUnit,ierr) close(statUnit) @@ -617,6 +620,7 @@ program DAMASK_spectral else guess = .true. ! start guessing after first converged (sub)inc endif + if (.not. cutBack) then if (worldrank == 0) then write(statUnit,*) totalIncsCounter, time, cutBackLevel, & @@ -624,23 +628,26 @@ program DAMASK_spectral flush(statUnit) endif endif - enddo subIncLooping + enddo subStepLooping + 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 - write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & + write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report converged inc ' increment ', totalIncsCounter, ' converged' else - write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report non-converged inc - ' increment ', totalIncsCounter, ' NOT converged' notConvergedCounter = notConvergedCounter + 1_pInt + write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report non-converged inc + ' increment ', totalIncsCounter, ' NOT converged' endif; flush(6) + if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency if (worldrank == 0) & write(6,'(1/,a)') ' ... writing results to file ......................................' call materialpoint_postResults() call MPI_file_seek (resUnit,fileOffset,MPI_SEEK_SET,ierr) - 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') 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) @@ -652,15 +659,12 @@ program DAMASK_spectral enddo fileOffset = fileOffset + sum(outputSize) ! forward to current file position endif - if( loadCases(currentLoadCase)%restartFrequency > 0_pInt .and. & ! at frequency of writing restart information set restart parameter for FEsolving - mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0_pInt) then ! first call to CPFEM_general will write? - restartWrite = .true. - lastRestartWritten = inc + if ( loadCases(currentLoadCase)%restartFrequency > 0_pInt & ! writing of restart info requested ... + .and. mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0_pInt) then ! ... and at frequency of writing restart information + restartWrite = .true. ! set restart parameter for FEsolving + lastRestartWritten = inc ! QUESTION: first call to CPFEM_general will write? endif - else forwarding - time = time + timeinc - guess = .true. - endif forwarding + endif skipping enddo incLooping enddo loadCaseLooping diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 2f451d953..5a8919b79 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -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) !$OMP FLUSH(crystallite_todo) #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 write(6,'(a,f12.8,a,i8,1x,i2,1x,i3,/)') '<< CRYST >> cutback step in crystallite_stressAndItsTangent & &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 >> max(subFrac) ',maxval(crystallite_subFrac) 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 ! --- integrate --- requires fully defined state array (basic + dependent state) @@ -2752,7 +2759,7 @@ subroutine crystallite_integrateStateFPI() enddo if (NaN) then ! NaN occured in any dotState 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... !$OMP CRITICAL (checkTodo) 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 !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then - - - crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e) !$OMP FLUSH(crystallite_todo) 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 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,/,(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 >> reltol dotstate',rTol_crystalliteState* & + abs(tempPlasticState(1:mySizePlasticDotState)) write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> new state',tempPlasticState(1:mySizePlasticDotState) endif #endif @@ -3202,9 +3209,9 @@ end function crystallite_push33ToRef !> intermediate acceleration of the Newton-Raphson correction !-------------------------------------------------------------------------------------------------- logical function crystallite_integrateStress(& - ipc,& ! grain number - ip,& ! integration point number - el,& ! element number + ipc,& ! grain number + ip,& ! integration point number + el,& ! element number timeFraction & ) use, intrinsic :: & @@ -3255,10 +3262,10 @@ logical function crystallite_integrateStress(& use mesh, only: mesh_element implicit none - integer(pInt), intent(in):: el, & ! element index - ip, & ! integration point index - ipc ! grain index - real(pReal), optional, intent(in) :: timeFraction ! fraction of timestep + integer(pInt), intent(in):: el, & ! element index + ip, & ! integration point index + ipc ! grain index + real(pReal), optional, intent(in) :: timeFraction ! fraction of timestep !*** local variables ***! real(pReal), dimension(3,3):: Fg_new, & ! deformation gradient at end of timestep @@ -3419,7 +3426,7 @@ logical function crystallite_integrateStress(& #ifdef DEBUG 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, & - ' 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 return endif loopsExeced @@ -3428,7 +3435,8 @@ logical function crystallite_integrateStress(& B = math_I3 - dt*Lpguess 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) !* 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) & 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, & 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 & .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.7,1x)/))') '<< CRYST >> Lp_constitutive', math_transpose33(Lp_constitutive) - write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lpguess', math_transpose33(Lpguess) + write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST >> Lp_constitutive', math_transpose33(Lp_constitutive) endif #endif @@ -3485,6 +3502,13 @@ logical function crystallite_integrateStress(& else ! not converged and residuum not improved... steplengthLp = subStepSizeLp * steplengthLp ! ...try with smaller step length in same direction 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 endif @@ -3498,6 +3522,16 @@ logical function crystallite_integrateStress(& dFe_dLp3333 = - dt * dFe_dLp3333 dRLp_dLp = math_identity2nd(9_pInt) & - 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 work = math_plain33to9(residuumLp) call dgesv(9,1,dRLp_dLp2,9,ipiv,work,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp diff --git a/src/spectral_mech_Basic.f90 b/src/spectral_mech_Basic.f90 index 55403ee7c..ea6526091 100644 --- a/src/spectral_mech_Basic.f90 +++ b/src/spectral_mech_Basic.f90 @@ -196,8 +196,9 @@ subroutine basicPETSc_init .false., & math_I3) 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) & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & '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 !-------------------------------------------------------------------------------------------------- -type(tSolutionState) function & - basicPETSc_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) +type(tSolutionState) function basicPETSc_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) use IO, only: & IO_error use numerics, only: & @@ -283,9 +283,8 @@ type(tSolutionState) function & CHKERRQ(ierr) basicPETSc_solution%termIll = terminallyIll terminallyIll = .false. - BasicPETSc_solution%converged =.true. if (reason == -4) call IO_error(893_pInt) - if (reason < 1) basicPETSC_solution%converged = .false. + BasicPETSc_solution%converged = reason > 0 basicPETSC_solution%iterationsNeeded = totalIter end function BasicPETSc_solution @@ -343,8 +342,8 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) - if(nfuncs== 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment - newIteration: if(totalIter <= PETScIter) then + if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment + newIteration: if (totalIter <= PETScIter) then !-------------------------------------------------------------------------------------------------- ! report begin of new iteration totalIter = totalIter + 1_pInt @@ -480,10 +479,10 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation 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 - if (restartWrite) then + 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 @@ -506,23 +505,23 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation endif endif - call utilities_updateIPcoords(F) + call utilities_updateIPcoords(F) ! QUESTION: why do this even when cutback happened?? - if (cutBack) then - F_aim = F_aim_lastInc - F = reshape(F_lastInc, [9,grid(1),grid(2),grid3]) + 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. + 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 + 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 + 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 + 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 @@ -531,8 +530,8 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation ! 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])) - 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]) endif F_aim = F_aim + f_aimDot * timeinc