diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index b3848a9eb..e34b5baa8 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -162,6 +162,7 @@ subroutine CPFEM_init write(6,'(/,a)') ' <<<+- CPFEM init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" + flush(6) endif mainProcess ! initialize stress and jacobian to zero @@ -242,8 +243,8 @@ subroutine CPFEM_init write(6,'(a32,1x,6(i8,1x))') 'CPFEM_dcsdE: ', shape(CPFEM_dcsdE) write(6,'(a32,1x,6(i8,1x),/)') 'CPFEM_dcsdE_knownGood: ', shape(CPFEM_dcsdE_knownGood) write(6,'(a32,l1)') 'symmetricSolver: ', symmetricSolver + flush(6) endif - flush(6) end subroutine CPFEM_init diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index 0ac916046..a16aee54f 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -9,7 +9,7 @@ module CPFEM2 private public :: & - CPFEM_general, & + CPFEM_age, & CPFEM_initAll contains @@ -127,6 +127,7 @@ subroutine CPFEM_init write(6,'(/,a)') ' <<<+- CPFEM init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" + flush(6) endif mainProcess ! *** restore the last converged values of each essential variable from the binary file @@ -194,7 +195,6 @@ subroutine CPFEM_init restartRead = .false. endif - flush(6) end subroutine CPFEM_init @@ -202,7 +202,7 @@ end subroutine CPFEM_init !-------------------------------------------------------------------------------------------------- !> @brief perform initialization at first call, update variables and call the actual material model !-------------------------------------------------------------------------------------------------- -subroutine CPFEM_general(age, dt) +subroutine CPFEM_age() use prec, only: & pReal, & pInt @@ -215,7 +215,6 @@ subroutine CPFEM_general(age, dt) debug_levelExtensive, & debug_levelSelective use FEsolving, only: & - terminallyIll, & restartWrite use math, only: & math_identity2nd, & @@ -254,114 +253,99 @@ subroutine CPFEM_general(age, dt) crystallite_dPdF, & crystallite_Tstar0_v, & crystallite_Tstar_v - use homogenization, only: & - materialpoint_stressAndItsTangent, & - materialpoint_postResults use IO, only: & IO_write_jobRealFile, & IO_warning use DAMASK_interface implicit none - real(pReal), intent(in) :: dt !< time increment - logical, intent(in) :: age !< age results integer(pInt) :: i, k, l, m, ph, homog, mySource - character(len=1024) :: rankStr + character(len=32) :: rankStr - !*** age results and write restart data if requested - if (age) then - crystallite_F0 = crystallite_partionedF ! crystallite deformation (_subF is perturbed...) - crystallite_Fp0 = crystallite_Fp ! crystallite plastic deformation - crystallite_Lp0 = crystallite_Lp ! crystallite plastic velocity - crystallite_Fi0 = crystallite_Fi ! crystallite intermediate deformation - crystallite_Li0 = crystallite_Li ! crystallite intermediate velocity - crystallite_dPdF0 = crystallite_dPdF ! crystallite stiffness - crystallite_Tstar0_v = crystallite_Tstar_v ! crystallite 2nd Piola Kirchhoff stress +if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) & + write(6,'(a)') '<< CPFEM >> aging states' - forall ( i = 1:size(plasticState )) plasticState(i)%state0 = plasticState(i)%state ! copy state in this lenghty way because: A component cannot be an array if the encompassing structure is an array - do i = 1, size(sourceState) - do mySource = 1,phase_Nsources(i) - sourceState(i)%p(mySource)%state0 = sourceState(i)%p(mySource)%state ! copy state in this lenghty way because: A component cannot be an array if the encompassing structure is an array - enddo; enddo - if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) & - write(6,'(a)') '<< CPFEM >> aging states' +crystallite_F0 = crystallite_partionedF ! crystallite deformation (_subF is perturbed...) +crystallite_Fp0 = crystallite_Fp ! crystallite plastic deformation +crystallite_Lp0 = crystallite_Lp ! crystallite plastic velocity +crystallite_Fi0 = crystallite_Fi ! crystallite intermediate deformation +crystallite_Li0 = crystallite_Li ! crystallite intermediate velocity +crystallite_dPdF0 = crystallite_dPdF ! crystallite stiffness +crystallite_Tstar0_v = crystallite_Tstar_v ! crystallite 2nd Piola Kirchhoff stress - do homog = 1_pInt, material_Nhomogenization - homogState (homog)%state0 = homogState (homog)%state - thermalState (homog)%state0 = thermalState (homog)%state - damageState (homog)%state0 = damageState (homog)%state - vacancyfluxState (homog)%state0 = vacancyfluxState (homog)%state - hydrogenfluxState(homog)%state0 = hydrogenfluxState(homog)%state - enddo +forall (i = 1:size(plasticState)) plasticState(i)%state0 = plasticState(i)%state ! copy state in this lengthy way because: A component cannot be an array if the encompassing structure is an array +do i = 1, size(sourceState) + do mySource = 1,phase_Nsources(i) + sourceState(i)%p(mySource)%state0 = sourceState(i)%p(mySource)%state ! copy state in this lengthy way because: A component cannot be an array if the encompassing structure is an array +enddo; enddo - if (restartWrite) then - if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) & - write(6,'(a)') '<< CPFEM >> writing state variables of last converged step to binary files' - - write(rankStr,'(a1,i0)')'_',worldrank +do homog = 1_pInt, material_Nhomogenization + homogState (homog)%state0 = homogState (homog)%state + thermalState (homog)%state0 = thermalState (homog)%state + damageState (homog)%state0 = damageState (homog)%state + vacancyfluxState (homog)%state0 = vacancyfluxState (homog)%state + hydrogenfluxState(homog)%state0 = hydrogenfluxState(homog)%state +enddo - call IO_write_jobRealFile(777,'recordedPhase'//trim(rankStr),size(material_phase)) - write (777,rec=1) material_phase - close (777) +if (restartWrite) then + if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) & + write(6,'(a)') '<< CPFEM >> writing state variables of last converged step to binary files' - call IO_write_jobRealFile(777,'convergedF'//trim(rankStr),size(crystallite_F0)) - write (777,rec=1) crystallite_F0 - close (777) + write(rankStr,'(a1,i0)')'_',worldrank - call IO_write_jobRealFile(777,'convergedFp'//trim(rankStr),size(crystallite_Fp0)) - write (777,rec=1) crystallite_Fp0 - close (777) + call IO_write_jobRealFile(777,'recordedPhase'//trim(rankStr),size(material_phase)) + write (777,rec=1) material_phase; close (777) - call IO_write_jobRealFile(777,'convergedFi'//trim(rankStr),size(crystallite_Fi0)) - write (777,rec=1) crystallite_Fi0 - close (777) + call IO_write_jobRealFile(777,'convergedF'//trim(rankStr),size(crystallite_F0)) + write (777,rec=1) crystallite_F0; close (777) - call IO_write_jobRealFile(777,'convergedLp'//trim(rankStr),size(crystallite_Lp0)) - write (777,rec=1) crystallite_Lp0 - close (777) + call IO_write_jobRealFile(777,'convergedFp'//trim(rankStr),size(crystallite_Fp0)) + write (777,rec=1) crystallite_Fp0; close (777) - call IO_write_jobRealFile(777,'convergedLi'//trim(rankStr),size(crystallite_Li0)) - write (777,rec=1) crystallite_Li0 - close (777) + call IO_write_jobRealFile(777,'convergedFi'//trim(rankStr),size(crystallite_Fi0)) + write (777,rec=1) crystallite_Fi0; close (777) - call IO_write_jobRealFile(777,'convergeddPdF'//trim(rankStr),size(crystallite_dPdF0)) - write (777,rec=1) crystallite_dPdF0 - close (777) + call IO_write_jobRealFile(777,'convergedLp'//trim(rankStr),size(crystallite_Lp0)) + write (777,rec=1) crystallite_Lp0; close (777) - call IO_write_jobRealFile(777,'convergedTstar'//trim(rankStr),size(crystallite_Tstar0_v)) - write (777,rec=1) crystallite_Tstar0_v - close (777) + call IO_write_jobRealFile(777,'convergedLi'//trim(rankStr),size(crystallite_Li0)) + write (777,rec=1) crystallite_Li0; close (777) - call IO_write_jobRealFile(777,'convergedStateConst'//trim(rankStr)) - m = 0_pInt - writePlasticityInstances: do ph = 1_pInt, size(phase_plasticity) - do k = 1_pInt, plasticState(ph)%sizeState - do l = 1, size(plasticState(ph)%state0(1,:)) - m = m+1_pInt - write(777,rec=m) plasticState(ph)%state0(k,l) - enddo; enddo - enddo writePlasticityInstances - close (777) + call IO_write_jobRealFile(777,'convergeddPdF'//trim(rankStr),size(crystallite_dPdF0)) + write (777,rec=1) crystallite_dPdF0; close (777) - call IO_write_jobRealFile(777,'convergedStateHomog'//trim(rankStr)) - m = 0_pInt - writeHomogInstances: do homog = 1_pInt, material_Nhomogenization - do k = 1_pInt, homogState(homog)%sizeState - do l = 1, size(homogState(homog)%state0(1,:)) - m = m+1_pInt - write(777,rec=m) homogState(homog)%state0(k,l) - enddo; enddo - enddo writeHomogInstances - close (777) + call IO_write_jobRealFile(777,'convergedTstar'//trim(rankStr),size(crystallite_Tstar0_v)) + write (777,rec=1) crystallite_Tstar0_v; close (777) - endif - endif + call IO_write_jobRealFile(777,'convergedStateConst'//trim(rankStr)) + m = 0_pInt + writePlasticityInstances: do ph = 1_pInt, size(phase_plasticity) + do k = 1_pInt, plasticState(ph)%sizeState + do l = 1, size(plasticState(ph)%state0(1,:)) + m = m+1_pInt + write(777,rec=m) plasticState(ph)%state0(k,l) + enddo; enddo + enddo writePlasticityInstances + close (777) - if (.not. terminallyIll) & - call materialpoint_stressAndItsTangent(.True., dt) + call IO_write_jobRealFile(777,'convergedStateHomog'//trim(rankStr)) + m = 0_pInt + writeHomogInstances: do homog = 1_pInt, material_Nhomogenization + do k = 1_pInt, homogState(homog)%sizeState + do l = 1, size(homogState(homog)%state0(1,:)) + m = m+1_pInt + write(777,rec=m) homogState(homog)%state0(k,l) + enddo; enddo + enddo writeHomogInstances + close (777) -end subroutine CPFEM_general +endif + +if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) & + write(6,'(a)') '<< CPFEM >> done aging states' + +end subroutine CPFEM_age end module CPFEM2 diff --git a/src/DAMASK_spectral.f90 b/src/DAMASK_spectral.f90 index f32bfb7b3..ac3fbf5a2 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,24 +454,23 @@ 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 +! looping over loadcases loadCaseLooping: do currentLoadCase = 1_pInt, size(loadCases) time0 = time ! currentLoadCase start time guess = loadCases(currentLoadCase)%followFormerTrajectory ! change of load case? homogeneous guess for the first inc !-------------------------------------------------------------------------------------------------- -! loop oper incs defined in input file for current currentLoadCase +! loop over incs defined in input file for current currentLoadCase incLooping: do inc = 1_pInt, loadCases(currentLoadCase)%incs totalIncsCounter = totalIncsCounter + 1_pInt !-------------------------------------------------------------------------------------------------- ! forwarding time - timeIncOld = timeinc + timeIncOld = timeinc ! last timeinc that brought former inc to an end if (loadCases(currentLoadCase)%logscale == 0_pInt) then ! linear scale - timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal) ! only valid for given linear time scale. will be overwritten later in case loglinear scale is used + timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal) else if (currentLoadCase == 1_pInt) then ! 1st currentLoadCase of logarithmic scale if (inc == 1_pInt) then ! 1st inc of 1st currentLoadCase of logarithmic scale @@ -486,20 +486,23 @@ program DAMASK_spectral real(loadCases(currentLoadCase)%incs ,pReal))) endif endif - timeinc = timeinc / 2.0_pReal**real(cutBackLevel,pReal) ! depending on cut back level, decrease time step + timeinc = timeinc / real(subStepFactor,pReal)**real(cutBackLevel,pReal) ! depending on cut back level, decrease time step - forwarding: if (totalIncsCounter >= restartInc) then - stepFraction = 0_pInt + 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)//& @@ -509,11 +512,11 @@ program DAMASK_spectral 's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,& '-', stepFraction, '/', subStepFactor**cutBackLevel,& ' of load case ', currentLoadCase,'/',size(loadCases) - flush(6) write(incInfo,'(a,'//IO_intOut(totalIncsCounter)//',a,'//IO_intOut(sum(loadCases%incs))//& ',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//')') & 'Increment ',totalIncsCounter,'/',sum(loadCases%incs),& '-',stepFraction, '/', subStepFactor**cutBackLevel + flush(6) !-------------------------------------------------------------------------------------------------- ! forward fields @@ -542,7 +545,7 @@ program DAMASK_spectral end select case(FIELD_THERMAL_ID); call spectral_thermal_forward() - case(FIELD_DAMAGE_ID); call spectral_damage_forward() + case(FIELD_DAMAGE_ID); call spectral_damage_forward() end select enddo @@ -582,65 +585,64 @@ 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) ! stationary with respect to staggered iteration 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. - 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 - guess = .true. ! accept non converged BVP solution - else ! default behavior, exit if spectral solver does not converge - 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 +! check solution for either advance or retry + + if ( (continueCalculation .or. all(solres(:)%converged .and. solres(:)%stagConverged)) & ! don't care or did converge + .and. .not. solres(1)%termIll) then ! and acceptable solution found + timeIncOld = timeinc + cutBack = .false. guess = .true. ! start guessing after first converged (sub)inc - endif - if (.not. cutBack) then if (worldrank == 0) then write(statUnit,*) totalIncsCounter, time, cutBackLevel, & - solres%converged, solres%iterationsNeeded ! write statistics about accepted solution + solres%converged, solres%iterationsNeeded flush(statUnit) endif + elseif (cutBackLevel < maxCutBack) then ! further cutbacking tolerated? + cutBack = .true. + stepFraction = (stepFraction - 1_pInt) * subStepFactor ! adjust to new denominator + cutBackLevel = cutBackLevel + 1_pInt + time = time - timeinc ! rewind time + timeinc = timeinc/real(subStepFactor,pReal) ! cut timestep + write(6,'(/,a)') ' cutting back ' + else ! no more options to continue + 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 - 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 ......................................' + flush(6) 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 +654,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 @@ -673,6 +672,7 @@ program DAMASK_spectral real(convergedCounter, pReal)/& real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, & ' %) increments converged!' + flush(6) call MPI_file_close(resUnit,ierr) close(statUnit) diff --git a/src/numerics.f90 b/src/numerics.f90 index 70c7f3c30..8392ac61c 100644 --- a/src/numerics.f90 +++ b/src/numerics.f90 @@ -120,9 +120,9 @@ module numerics petsc_options = '' integer(pInt), protected, public :: & fftw_planner_flag = 32_pInt, & !< conversion of fftw_plan_mode to integer, basically what is usually done in the include file of fftw - continueCalculation = 0_pInt, & !< 0: exit if BVP solver does not converge, 1: continue calculation if BVP solver does not converge divergence_correction = 2_pInt !< correct divergence calculation in fourier space 0: no correction, 1: size scaled to 1, 2: size scaled to Npoints logical, protected, public :: & + continueCalculation = .false., & !< false:exit if BVP solver does not converge, true: continue calculation despite BVP solver not converging memory_efficient = .true., & !< for fast execution (pre calculation of gamma_hat), Default .true.: do not precalculate update_gamma = .false. !< update gamma operator with current stiffness, Default .false.: use initial stiffness #endif @@ -424,9 +424,9 @@ subroutine numerics_init case ('err_stress_tolabs') err_stress_tolabs = IO_floatValue(line,chunkPos,2_pInt) case ('continuecalculation') - continueCalculation = IO_intValue(line,chunkPos,2_pInt) + continueCalculation = IO_intValue(line,chunkPos,2_pInt) > 0_pInt case ('memory_efficient') - memory_efficient = IO_intValue(line,chunkPos,2_pInt) > 0_pInt + memory_efficient = IO_intValue(line,chunkPos,2_pInt) > 0_pInt case ('fftw_timelimit') fftw_timelimit = IO_floatValue(line,chunkPos,2_pInt) case ('fftw_plan_mode') @@ -436,7 +436,7 @@ subroutine numerics_init case ('divergence_correction') divergence_correction = IO_intValue(line,chunkPos,2_pInt) case ('update_gamma') - update_gamma = IO_intValue(line,chunkPos,2_pInt) > 0_pInt + update_gamma = IO_intValue(line,chunkPos,2_pInt) > 0_pInt case ('petsc_options') petsc_options = trim(line(chunkPos(4):)) case ('spectralsolver','myspectralsolver') @@ -599,7 +599,7 @@ subroutine numerics_init !-------------------------------------------------------------------------------------------------- ! spectral parameters #ifdef Spectral - write(6,'(a24,1x,i8)') ' continueCalculation: ',continueCalculation + write(6,'(a24,1x,L8)') ' continueCalculation: ',continueCalculation write(6,'(a24,1x,L8)') ' memory_efficient: ',memory_efficient write(6,'(a24,1x,i8)') ' divergence_correction: ',divergence_correction write(6,'(a24,1x,a)') ' spectral_derivative: ',trim(spectral_derivative) @@ -698,8 +698,6 @@ subroutine numerics_init if (err_hydrogenflux_tolabs <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_hydrogenflux_tolabs') if (err_hydrogenflux_tolrel <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_hydrogenflux_tolrel') #ifdef Spectral - if (continueCalculation /= 0_pInt .and. & - continueCalculation /= 1_pInt) call IO_error(301_pInt,ext_msg='continueCalculation') if (divergence_correction < 0_pInt .or. & divergence_correction > 2_pInt) call IO_error(301_pInt,ext_msg='divergence_correction') if (update_gamma .and. & @@ -713,7 +711,7 @@ subroutine numerics_init if (polarAlpha <= 0.0_pReal .or. & polarAlpha > 2.0_pReal) call IO_error(301_pInt,ext_msg='polarAlpha') if (polarBeta < 0.0_pReal .or. & - polarBeta > 2.0_pReal) call IO_error(301_pInt,ext_msg='polarBeta') + polarBeta > 2.0_pReal) call IO_error(301_pInt,ext_msg='polarBeta') #endif end subroutine numerics_init diff --git a/src/spectral_mech_AL.f90 b/src/spectral_mech_AL.f90 index 6d0fff286..4695d4faa 100644 --- a/src/spectral_mech_AL.f90 +++ b/src/spectral_mech_AL.f90 @@ -213,8 +213,9 @@ subroutine AL_init endif restart call Utilities_updateIPcoords(reshape(F,shape(F_lastInc))) - call Utilities_constitutiveResponse(F_lastInc, reshape(F,shape(F_lastInc)), & - 0.0_pReal,P,C_volAvg,C_minMaxAvg,temp33_Real,.false.,math_I3) + call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & + reshape(F,shape(F_lastInc)), 0.0_pReal, math_I3) + nullify(F) nullify(F_lambda) call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! write data back to PETSc @@ -364,12 +365,10 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) DMDALocalInfo, dimension(& DMDA_LOCAL_INFO_SIZE) :: & in - PetscScalar, target, dimension(3,3,2, & - XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: & - x_scal - PetscScalar, target, dimension(3,3,2, & - X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: & - f_scal + PetscScalar, & + target, dimension(3,3,2, XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: x_scal + PetscScalar, & + target, dimension(3,3,2, X_RANGE, Y_RANGE, Z_RANGE), intent(out) :: f_scal PetscScalar, pointer, dimension(:,:,:,:,:) :: & F, & F_lambda, & @@ -441,8 +440,9 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response P_avLastEval = P_av - call Utilities_constitutiveResponse(F_lastInc,F - residual_F_lambda/polarBeta,params%timeinc, & - residual_F,C_volAvg,C_minMaxAvg,P_av,ForwardData,params%rotation_BC) + + call Utilities_constitutiveResponse(residual_F,P_av,C_volAvg,C_minMaxAvg, & + F - residual_F_lambda/polarBeta,params%timeinc, params%rotation_BC) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) ForwardData = .False. @@ -655,10 +655,12 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stre !-------------------------------------------------------------------------------------------------- ! 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_lambdaDot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & - timeinc_old,guess,F_lambda_lastInc,reshape(F_lambda,[3,3,grid(1),grid(2),grid3])) + Fdot = Utilities_calculateRate(guess, & + F_lastInc, reshape(F, [3,3,grid(1),grid(2),grid3]), timeinc_old, & + math_rotate_backward33(f_aimDot,rotation_BC)) + F_lambdaDot = Utilities_calculateRate(guess, & + F_lambda_lastInc,reshape(F_lambda,[3,3,grid(1),grid(2),grid3]), timeinc_old, & + math_rotate_backward33(f_aimDot,rotation_BC)) F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) F_lambda_lastInc = reshape(F_lambda,[3,3,grid(1),grid(2),grid3]) endif diff --git a/src/spectral_mech_Basic.f90 b/src/spectral_mech_Basic.f90 index 55403ee7c..acdcbee3a 100644 --- a/src/spectral_mech_Basic.f90 +++ b/src/spectral_mech_Basic.f90 @@ -39,16 +39,16 @@ module spectral_mech_basic ! stress, stiffness and compliance average etc. real(pReal), private, dimension(3,3) :: & F_aim = math_I3, & - F_aim_lastIter = math_I3, & F_aim_lastInc = math_I3, & P_av = 0.0_pReal, & - F_aimDot=0.0_pReal + F_aimDot = 0.0_pReal character(len=1024), private :: incInfo real(pReal), private, dimension(3,3,3,3) :: & C_volAvg = 0.0_pReal, & !< current volume average stiffness C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 stiffness - S = 0.0_pReal !< current compliance (filled up with zeros) + C_minMaxAvgLastInc = 0.0_pReal, & !< previous (min+max)/2 stiffness + S = 0.0_pReal !< current compliance (filled up with zeros) real(pReal), private :: err_stress, err_div logical, private :: ForwardData integer(pInt), private :: & @@ -69,7 +69,7 @@ module spectral_mech_basic contains !-------------------------------------------------------------------------------------------------- -!> @brief allocates all neccessary fields and fills them with data, potentially from restart info +!> @brief allocates all necessary fields and fills them with data, potentially from restart info !-------------------------------------------------------------------------------------------------- subroutine basicPETSc_init #ifdef __GFORTRAN__ @@ -90,6 +90,8 @@ subroutine basicPETSc_init use numerics, only: & worldrank, & worldsize + use homogenization, only: & + materialpoint_F0 use DAMASK_interface, only: & getSolverJobName use spectral_utilities, only: & @@ -172,14 +174,11 @@ subroutine basicPETSc_init flush(6) write(rankStr,'(a1,i0)')'_',worldrank call IO_read_realFile(777,'F'//trim(rankStr),trim(getSolverJobName()),size(F)) - read (777,rec=1) F - close (777) + read (777,rec=1) F; close (777) call IO_read_realFile(777,'F_lastInc'//trim(rankStr),trim(getSolverJobName()),size(F_lastInc)) - read (777,rec=1) F_lastInc - close (777) + read (777,rec=1) F_lastInc; close (777) call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot)) - read (777,rec=1) f_aimDot - close (777) + read (777,rec=1) f_aimDot; close (777) F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc elseif (restartInc == 1_pInt) then restart @@ -187,41 +186,36 @@ subroutine basicPETSc_init F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) endif restart + materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent call Utilities_updateIPcoords(reshape(F,shape(F_lastInc))) - call Utilities_constitutiveResponse(F_lastInc, reshape(F,shape(F_lastInc)), & - 0.0_pReal, & - P, & - C_volAvg,C_minMaxAvg, & ! global average of stiffness and (min+max)/2 - temp33_Real, & - .false., & - math_I3) + call Utilities_constitutiveResponse(P, temp33_Real, C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 + reshape(F,shape(F_lastInc)), & ! target F + 0.0_pReal, & ! time increment + math_I3) ! no rotation of boundary condition 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' + 'reading more values of increment', restartInc-1_pInt, 'from file' flush(6) call IO_read_realFile(777,'C_volAvg',trim(getSolverJobName()),size(C_volAvg)) - read (777,rec=1) C_volAvg - close (777) + read (777,rec=1) C_volAvg; close (777) call IO_read_realFile(777,'C_volAvgLastInc',trim(getSolverJobName()),size(C_volAvgLastInc)) - read (777,rec=1) C_volAvgLastInc - close (777) + read (777,rec=1) C_volAvgLastInc; close (777) call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg)) - read (777,rec=1) C_minMaxAvg - close (777) + read (777,rec=1) C_minMaxAvg; close (777) endif restartRead - call Utilities_updateGamma(C_minmaxAvg,.True.) + call Utilities_updateGamma(C_minmaxAvg,.true.) 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: & @@ -238,13 +232,13 @@ type(tSolutionState) function & !-------------------------------------------------------------------------------------------------- ! input data for solution - real(pReal), intent(in) :: & - timeinc, & !< increment in time for current solution - timeinc_old !< increment in time of last increment - type(tBoundaryCondition), intent(in) :: & - stress_BC character(len=*), intent(in) :: & incInfoIn + real(pReal), intent(in) :: & + timeinc, & !< increment time for current solution + timeinc_old !< increment time of last successful increment + type(tBoundaryCondition), intent(in) :: & + stress_BC real(pReal), dimension(3,3), intent(in) :: rotation_BC !-------------------------------------------------------------------------------------------------- @@ -279,14 +273,13 @@ type(tSolutionState) function & !-------------------------------------------------------------------------------------------------- ! check convergence - call SNESGetConvergedReason(snes,reason,ierr) - CHKERRQ(ierr) + call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr) + + BasicPETSc_solution%converged = reason > 0 + basicPETSC_solution%iterationsNeeded = totalIter 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%iterationsNeeded = totalIter + if (reason == -4) call IO_error(893_pInt) ! MPI error end function BasicPETSc_solution @@ -322,19 +315,18 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) terminallyIll implicit none - DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & - in - PetscScalar, dimension(3,3, & - XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: & - x_scal - PetscScalar, dimension(3,3, & - X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: & - f_scal + DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in + PetscScalar, & + dimension(3,3, XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: x_scal !< what is this? + PetscScalar, & + dimension(3,3, X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: f_scal !< what is this? PetscInt :: & PETScIter, & nfuncs PetscObject :: dummy PetscErrorCode :: ierr + real(pReal), dimension(3,3) :: & + deltaF_aim external :: & SNESGetNumberFunctionEvals, & @@ -343,46 +335,45 @@ 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 !-------------------------------------------------------------------------------------------------- -! report begin of new iteration +! begin of new iteration + newIteration: if (totalIter <= PETScIter) then totalIter = totalIter + 1_pInt - write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & - ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax + write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') & + trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', & - math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', & - math_transpose33(F_aim) + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & + ' deformation gradient aim (lab) =', math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & + ' deformation gradient aim =', math_transpose33(F_aim) flush(6) endif newIteration !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - call Utilities_constitutiveResponse(F_lastInc,x_scal,params%timeinc, & - f_scal,C_volAvg,C_minmaxAvg,P_av,ForwardData,params%rotation_BC) + call Utilities_constitutiveResponse(f_scal,P_av,C_volAvg,C_minmaxAvg, & + x_scal,params%timeinc, params%rotation_BC) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) - ForwardData = .false. !-------------------------------------------------------------------------------------------------- ! stress BC handling - F_aim_lastIter = F_aim - F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc - err_stress = maxval(abs(mask_stress * (P_av - params%stress_BC))) ! mask = 0.0 for no bc + deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC) + F_aim = F_aim - deltaF_aim + err_stress = maxval(abs(mask_stress * (P_av - params%stress_BC))) ! mask = 0.0 when no stress bc !-------------------------------------------------------------------------------------------------- ! updated deformation gradient using fix point algorithm of basic scheme tensorField_real = 0.0_pReal tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = f_scal - call utilities_FFTtensorForward() - err_div = Utilities_divergenceRMS() - call utilities_fourierGammaConvolution(math_rotate_backward33(F_aim_lastIter-F_aim,params%rotation_BC)) - call utilities_FFTtensorBackward() + call utilities_FFTtensorForward() ! FFT forward of global "tensorField_real" + err_div = Utilities_divergenceRMS() ! divRMS of tensorField_fourier + call utilities_fourierGammaConvolution(math_rotate_backward33(deltaF_aim,params%rotation_BC)) ! convolution of Gamma and tensorField_fourier, with arg + call utilities_FFTtensorBackward() ! FFT backward of global tensorField_fourier !-------------------------------------------------------------------------------------------------- ! constructing residual - f_scal = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) + f_scal = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too end subroutine BasicPETSc_formResidual @@ -443,106 +434,120 @@ end subroutine BasicPETSc_converged !-------------------------------------------------------------------------------------------------- !> @brief forwarding routine +!> @details find new boundary conditions and best F estimate for end of current timestep +!> 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) - 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 + use math, only: & + math_mul33x33 ,& + math_rotate_backward33 + use numerics, only: & + worldrank + use homogenization, only: & + materialpoint_F0 + use mesh, only: & + grid, & + grid3 + use CPFEM2, only: & + CPFEM_age + 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(:,:,:,:) + implicit none + logical, intent(in) :: & + guess + 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 + PetscErrorCode :: ierr + PetscScalar, pointer :: F(:,:,:,:) - character(len=1024) :: rankStr + character(len=32) :: rankStr - call DMDAVecGetArrayF90(da,solution_vec,F,ierr) -!-------------------------------------------------------------------------------------------------- -! restart information for spectral solver - if (restartWrite) then - 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 DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) + + if (cutBack) then + C_volAvg = C_volAvgLastInc ! QUESTION: where is this required? + C_minMaxAvg = C_minMaxAvgLastInc ! QUESTION: where is this required? + else + !-------------------------------------------------------------------------------------------------- + ! 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) - call utilities_updateIPcoords(F) + if (worldrank == 0_pInt) then + 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) + call IO_write_jobRealFile(777,'C_minMaxAvg',size(C_volAvg)) + write (777,rec=1) C_minMaxAvg; close(777) + call IO_write_jobRealFile(777,'C_minMaxAvgLastInc',size(C_volAvgLastInc)) + write (777,rec=1) C_minMaxAvgLastInc; close(777) + endif - if (cutBack) then - F_aim = F_aim_lastInc - F = reshape(F_lastInc, [9,grid(1),grid(2),grid3]) - C_volAvg = C_volAvgLastInc - else - ForwardData = .True. - 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 + 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) + endif + + call CPFEM_age() ! age state and kinematics + call utilities_updateIPcoords(F) + + C_volAvgLastInc = C_volAvg + C_minMaxAvgLastInc = C_minMaxAvg + + if (guess) then ! QUESTION: better with a = L ? x:y + F_aimDot = stress_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old ! initialize with correction based on last inc + else + F_aimDot = 0.0_pReal + endif + F_aim_lastInc = F_aim + !-------------------------------------------------------------------------------------------------- + ! calculate rate for aim + if (deformation_BC%myType=='l') then ! calculate f_aimDot from given L and current F + F_aimDot = & + F_aimDot + deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim_lastInc) + elseif(deformation_BC%myType=='fdot') then ! f_aimDot is prescribed + F_aimDot = & + F_aimDot + deformation_BC%maskFloat * deformation_BC%values + elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed + F_aimDot = & + F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime + endif + + + Fdot = Utilities_calculateRate(guess, & + F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, & + math_rotate_backward33(f_aimDot,rotation_BC)) + F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) ! winding F forward + materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + endif !-------------------------------------------------------------------------------------------------- -! 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]) - 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) - +! update average and local deformation gradients + F_aim = F_aim_lastInc + f_aimDot * timeinc + 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) + end subroutine BasicPETSc_forward !-------------------------------------------------------------------------------------------------- diff --git a/src/spectral_mech_Polarisation.f90 b/src/spectral_mech_Polarisation.f90 index ecf707d46..fc65f14cf 100644 --- a/src/spectral_mech_Polarisation.f90 +++ b/src/spectral_mech_Polarisation.f90 @@ -213,8 +213,8 @@ subroutine Polarisation_init endif restart call Utilities_updateIPcoords(reshape(F,shape(F_lastInc))) - call Utilities_constitutiveResponse(F_lastInc, reshape(F,shape(F_lastInc)), & - 0.0_pReal,P,C_volAvg,C_minMaxAvg,temp33_Real,.false.,math_I3) + call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & + reshape(F,shape(F_lastInc)),0.0_pReal,math_I3) nullify(F) nullify(F_tau) call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! write data back to PETSc @@ -364,12 +364,10 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) DMDALocalInfo, dimension(& DMDA_LOCAL_INFO_SIZE) :: & in - PetscScalar, target, dimension(3,3,2, & - XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: & - x_scal - PetscScalar, target, dimension(3,3,2, & - X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: & - f_scal + PetscScalar, & + target, dimension(3,3,2, XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: x_scal + PetscScalar, & + target, dimension(3,3,2, X_RANGE, Y_RANGE, Z_RANGE), intent(out) :: f_scal PetscScalar, pointer, dimension(:,:,:,:,:) :: & F, & F_tau, & @@ -440,8 +438,8 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response P_avLastEval = P_av - call Utilities_constitutiveResponse(F_lastInc,F - residual_F_tau/polarBeta,params%timeinc, & - residual_F,C_volAvg,C_minMaxAvg,P_av,ForwardData,params%rotation_BC) + call Utilities_constitutiveResponse(residual_F,P_av,C_volAvg,C_minMaxAvg, & + F - residual_F_tau/polarBeta,params%timeinc,params%rotation_BC) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) ForwardData = .False. @@ -654,13 +652,13 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,deformati !-------------------------------------------------------------------------------------------------- ! 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_tauDot = Utilities_calculateRate(math_rotate_backward33(2.0_pReal*f_aimDot,rotation_BC), & - timeinc_old,guess,F_tau_lastInc, & - reshape(F_tau,[3,3,grid(1),grid(2),grid3])) - F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) + Fdot = Utilities_calculateRate(guess, & + F_lastInc, reshape(F, [3,3,grid(1),grid(2),grid3]), timeinc_old, & + math_rotate_backward33( f_aimDot,rotation_BC)) + F_tauDot = Utilities_calculateRate(guess, & + F_tau_lastInc, reshape(F_tau,[3,3,grid(1),grid(2),grid3]), timeinc_old, & + math_rotate_backward33(2.0_pReal*f_aimDot,rotation_BC)) + F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3]) endif diff --git a/src/spectral_utilities.f90 b/src/spectral_utilities.f90 index 1bbf2e608..3295aa2bd 100644 --- a/src/spectral_utilities.f90 +++ b/src/spectral_utilities.f90 @@ -16,7 +16,7 @@ module spectral_utilities #include include 'fftw3-mpi.f03' - logical, public :: cutBack =.false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill + logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill integer(pInt), public, parameter :: maxPhaseFields = 2_pInt integer(pInt), public :: nActiveFields = 0_pInt @@ -799,7 +799,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) call math_invert(size_reduced, c_reduced, s_reduced, errmatinv) ! invert reduced stiffness if (any(IEEE_is_NaN(s_reduced))) errmatinv = .true. - if(errmatinv) call IO_error(error_ID=400_pInt,ext_msg='utilities_maskedCompliance') + if (errmatinv) call IO_error(error_ID=400_pInt,ext_msg='utilities_maskedCompliance') temp99_Real = 0.0_pReal ! fill up compliance with zeros k = 0_pInt do n = 1_pInt,9_pInt @@ -817,28 +817,30 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) sTimesC = matmul(c_reduced,s_reduced) do m=1_pInt, size_reduced do n=1_pInt, size_reduced - if(m==n .and. abs(sTimesC(m,n)) > (1.0_pReal + 10.0e-12_pReal)) errmatinv = .true. ! diagonal elements of S*C should be 1 - if(m/=n .and. abs(sTimesC(m,n)) > (0.0_pReal + 10.0e-12_pReal)) errmatinv = .true. ! off diagonal elements of S*C should be 0 + errmatinv = errmatinv & + .or. (m==n .and. abs(sTimesC(m,n)-1.0_pReal) > 1.0e-12_pReal) & ! diagonal elements of S*C should be 1 + .or. (m/=n .and. abs(sTimesC(m,n)) > 1.0e-12_pReal) ! off-diagonal elements of S*C should be 0 enddo enddo - if(debugGeneral .or. errmatinv) then - write(formatString, '(I16.16)') size_reduced + if (debugGeneral .or. errmatinv) then + write(formatString, '(i2)') size_reduced formatString = '(/,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))' write(6,trim(formatString),advance='no') ' C * S (load) ', & transpose(matmul(c_reduced,s_reduced)) write(6,trim(formatString),advance='no') ' S (load) ', transpose(s_reduced) + if(errmatinv) call IO_error(error_ID=400_pInt,ext_msg='utilities_maskedCompliance') endif - if(errmatinv) call IO_error(error_ID=400_pInt,ext_msg='utilities_maskedCompliance') deallocate(c_reduced) deallocate(s_reduced) deallocate(sTimesC) else temp99_real = 0.0_pReal endif - if(debugGeneral) & - write(6,'(/,a,/,9(9(2x,f12.7,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 utilities_maskedCompliance = math_Plain99to3333(temp99_Real) end function utilities_maskedCompliance @@ -924,10 +926,10 @@ end subroutine utilities_fourierTensorDivergence !-------------------------------------------------------------------------------------------------- -!> @brief calculates constitutive response +!> @brief calculate constitutive response from materialpoint_F0 to F during timeinc !-------------------------------------------------------------------------------------------------- -subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, & - P,C_volAvg,C_minmaxAvg,P_av,forwardData,rotation_BC) +subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& + F,timeinc,rotation_BC) use IO, only: & IO_error use debug, only: & @@ -940,31 +942,22 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, & use mesh, only: & grid,& grid3 - use FEsolving, only: & - restartWrite - use CPFEM2, only: & - CPFEM_general use homogenization, only: & - materialpoint_F0, & materialpoint_F, & materialpoint_P, & - materialpoint_dPdF + materialpoint_dPdF, & + materialpoint_stressAndItsTangent implicit none - real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: & - F_lastInc, & !< target deformation gradient - F !< previous deformation gradient - real(pReal), intent(in) :: timeinc !< loading time - logical, intent(in) :: forwardData !< age results - real(pReal), intent(in), dimension(3,3) :: rotation_BC !< rotation of load frame - real(pReal),intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress real(pReal),intent(out), dimension(3,3,grid(1),grid(2),grid3) :: P !< PK stress - logical :: & - age + real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: F !< deformation gradient target !< previous deformation gradient + real(pReal), intent(in) :: timeinc !< loading time + real(pReal), intent(in), dimension(3,3) :: rotation_BC !< rotation of load frame + integer(pInt) :: & j,k,ierr real(pReal), dimension(3,3,3,3) :: max_dPdF, min_dPdF @@ -975,17 +968,9 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, & write(6,'(/,a)') ' ... evaluating constitutive response ......................................' flush(6) - age = .False. - - if (forwardData) then ! aging results - age = .True. - materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) - endif - if (cutBack) age = .False. ! restore saved variables - - materialpoint_F = reshape(F,[3,3,1,product(grid(1:2))*grid3]) - call debug_reset() ! this has no effect on rank >0 - + + materialpoint_F = reshape(F,[3,3,1,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field + !-------------------------------------------------------------------------------------------------- ! calculate bounds of det(F) and report if(debugGeneral) then @@ -1002,7 +987,19 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, & flush(6) endif - call CPFEM_general(age,timeinc) + call debug_reset() ! this has no effect on rank >0 + call materialpoint_stressAndItsTangent(.true.,timeinc) ! calculate P field + + P = reshape(materialpoint_P, [3,3,grid(1),grid(2),grid3]) + P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P + call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + if (debugRotation) & + write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress (lab) / MPa =',& + math_transpose33(P_av)*1.e-6_pReal + P_av = math_rotate_forward33(P_av,rotation_BC) + write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',& + math_transpose33(P_av)*1.e-6_pReal + flush(6) max_dPdF = 0.0_pReal max_dPdF_norm = 0.0_pReal @@ -1020,38 +1017,24 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, & end do call MPI_Allreduce(MPI_IN_PLACE,max_dPdF,81,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) - if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce max') + if (ierr /= 0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce max') call MPI_Allreduce(MPI_IN_PLACE,min_dPdF,81,MPI_DOUBLE,MPI_MIN,PETSC_COMM_WORLD,ierr) - if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce min') + if (ierr /= 0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce min') C_minmaxAvg = 0.5_pReal*(max_dPdF + min_dPdF) - C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) * wgt + C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) * wgt call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) call debug_info() ! this has no effect on rank >0 - restartWrite = .false. ! reset restartWrite status - cutBack = .false. ! reset cutBack status - - P = reshape(materialpoint_P, [3,3,grid(1),grid(2),grid3]) - P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P - call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - if (debugRotation) & - write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress (lab) / MPa =',& - math_transpose33(P_av)*1.e-6_pReal - P_av = math_rotate_forward33(P_av,rotation_BC) - write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',& - math_transpose33(P_av)*1.e-6_pReal - flush(6) - end subroutine utilities_constitutiveResponse !-------------------------------------------------------------------------------------------------- !> @brief calculates forward rate, either guessing or just add delta/timeinc !-------------------------------------------------------------------------------------------------- -pure function utilities_calculateRate(avRate,timeinc_old,guess,field_lastInc,field) +pure function utilities_calculateRate(heterogeneous,field0,field,dt,avRate) use mesh, only: & grid3, & grid @@ -1059,17 +1042,17 @@ pure function utilities_calculateRate(avRate,timeinc_old,guess,field_lastInc,fie implicit none real(pReal), intent(in), dimension(3,3) :: avRate !< homogeneous addon real(pReal), intent(in) :: & - timeinc_old !< timeinc of last step + dt !< timeinc between field0 and field logical, intent(in) :: & - guess !< guess along former trajectory + heterogeneous !< calculate field of rates real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: & - field_lastInc, & !< data of previous step + field0, & !< data of previous step field !< data of current step real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: & utilities_calculateRate - if (guess) then - utilities_calculateRate = (field-field_lastInc) / timeinc_old + if (heterogeneous) then + utilities_calculateRate = (field-field0) / dt else utilities_calculateRate = spread(spread(spread(avRate,3,grid(1)),4,grid(2)),5,grid3) endif