diff --git a/examples/SpectralMethod/20grains.seeds b/examples/SpectralMethod/Polycrystal/20grains.seeds similarity index 100% rename from examples/SpectralMethod/20grains.seeds rename to examples/SpectralMethod/Polycrystal/20grains.seeds diff --git a/examples/SpectralMethod/20grains16x16x16.geom b/examples/SpectralMethod/Polycrystal/20grains16x16x16.geom similarity index 100% rename from examples/SpectralMethod/20grains16x16x16.geom rename to examples/SpectralMethod/Polycrystal/20grains16x16x16.geom diff --git a/examples/SpectralMethod/20grains32x32x32.geom b/examples/SpectralMethod/Polycrystal/20grains32x32x32.geom similarity index 100% rename from examples/SpectralMethod/20grains32x32x32.geom rename to examples/SpectralMethod/Polycrystal/20grains32x32x32.geom diff --git a/examples/SpectralMethod/20grains64x64x64.geom b/examples/SpectralMethod/Polycrystal/20grains64x64x64.geom similarity index 100% rename from examples/SpectralMethod/20grains64x64x64.geom rename to examples/SpectralMethod/Polycrystal/20grains64x64x64.geom diff --git a/examples/SpectralMethod/material.config b/examples/SpectralMethod/Polycrystal/material.config similarity index 100% rename from examples/SpectralMethod/material.config rename to examples/SpectralMethod/Polycrystal/material.config diff --git a/examples/SpectralMethod/numerics.config b/examples/SpectralMethod/Polycrystal/numerics.config similarity index 100% rename from examples/SpectralMethod/numerics.config rename to examples/SpectralMethod/Polycrystal/numerics.config diff --git a/examples/SpectralMethod/shearXY.load b/examples/SpectralMethod/Polycrystal/shearXY.load similarity index 100% rename from examples/SpectralMethod/shearXY.load rename to examples/SpectralMethod/Polycrystal/shearXY.load diff --git a/examples/SpectralMethod/shearZX.load b/examples/SpectralMethod/Polycrystal/shearZX.load similarity index 100% rename from examples/SpectralMethod/shearZX.load rename to examples/SpectralMethod/Polycrystal/shearZX.load diff --git a/examples/SpectralMethod/tensionX.load b/examples/SpectralMethod/Polycrystal/tensionX.load similarity index 100% rename from examples/SpectralMethod/tensionX.load rename to examples/SpectralMethod/Polycrystal/tensionX.load diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index aa2615b5f..a89bfc294 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 @@ -139,36 +140,28 @@ subroutine CPFEM_init write(rankStr,'(a1,i0)')'_',worldrank call IO_read_intFile(777,'recordedPhase'//trim(rankStr),modelName,size(material_phase)) - read (777,rec=1) material_phase - close (777) + read (777,rec=1) material_phase; close (777) call IO_read_realFile(777,'convergedF'//trim(rankStr),modelName,size(crystallite_F0)) - read (777,rec=1) crystallite_F0 - close (777) + read (777,rec=1) crystallite_F0; close (777) call IO_read_realFile(777,'convergedFp'//trim(rankStr),modelName,size(crystallite_Fp0)) - read (777,rec=1) crystallite_Fp0 - close (777) + read (777,rec=1) crystallite_Fp0; close (777) call IO_read_realFile(777,'convergedFi'//trim(rankStr),modelName,size(crystallite_Fi0)) - read (777,rec=1) crystallite_Fi0 - close (777) + read (777,rec=1) crystallite_Fi0; close (777) call IO_read_realFile(777,'convergedLp'//trim(rankStr),modelName,size(crystallite_Lp0)) - read (777,rec=1) crystallite_Lp0 - close (777) + read (777,rec=1) crystallite_Lp0; close (777) call IO_read_realFile(777,'convergedLi'//trim(rankStr),modelName,size(crystallite_Li0)) - read (777,rec=1) crystallite_Li0 - close (777) + read (777,rec=1) crystallite_Li0; close (777) call IO_read_realFile(777,'convergeddPdF'//trim(rankStr),modelName,size(crystallite_dPdF0)) - read (777,rec=1) crystallite_dPdF0 - close (777) + read (777,rec=1) crystallite_dPdF0; close (777) call IO_read_realFile(777,'convergedTstar'//trim(rankStr),modelName,size(crystallite_Tstar0_v)) - read (777,rec=1) crystallite_Tstar0_v - close (777) + read (777,rec=1) crystallite_Tstar0_v; close (777) call IO_read_realFile(777,'convergedStateConst'//trim(rankStr),modelName) m = 0_pInt @@ -194,7 +187,6 @@ subroutine CPFEM_init restartRead = .false. endif - flush(6) end subroutine CPFEM_init @@ -202,7 +194,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 +207,6 @@ subroutine CPFEM_general(age, dt) debug_levelExtensive, & debug_levelSelective use FEsolving, only: & - terminallyIll, & restartWrite use math, only: & math_identity2nd, & @@ -254,114 +245,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 old mode 100644 new mode 100755 index fe7d78621..7ec2dde4c --- a/src/DAMASK_spectral.f90 +++ b/src/DAMASK_spectral.f90 @@ -78,8 +78,7 @@ program DAMASK_spectral FIELD_UNDEFINED_ID, & FIELD_MECH_ID, & FIELD_THERMAL_ID, & - FIELD_DAMAGE_ID, & - utilities_calcPlasticity + FIELD_DAMAGE_ID use spectral_mech_Basic use spectral_mech_AL use spectral_mech_Polarisation @@ -157,19 +156,6 @@ program DAMASK_spectral MPI_finalize, & MPI_allreduce, & PETScFinalize -!-------------------------------------------------------------------------------------------------- -! variables related to stop criterion for yielding - real(pReal) :: plasticWorkOld, plasticWorkNew, & ! plastic work - eqTotalStrainOld, eqTotalStrainNew, & ! total equivalent strain - eqPlasticStrainOld, eqPlasticStrainNew, & ! total equivalent plastic strain - eqStressOld, eqStressNew , & ! equivalent stress - yieldStopValue - real(pReal), dimension(3,3) :: yieldStress,yieldStressOld,yieldStressNew, & - plasticStrainOld, plasticStrainNew, plasticStrainRate - integer(pInt) :: yieldResUnit = 0_pInt - integer(pInt) :: stressstrainUnit = 0_pInt - character(len=13) :: stopFlag - logical :: yieldStop, yieldStopSatisfied !-------------------------------------------------------------------------------------------------- ! init DAMASK (all modules) @@ -227,8 +213,6 @@ program DAMASK_spectral !-------------------------------------------------------------------------------------------------- ! reading the load case and assign values to the allocated data structure - yieldStop = .False. - yieldStopSatisfied = .False. rewind(FILEUNIT) do line = IO_read(FILEUNIT) @@ -303,30 +287,10 @@ program DAMASK_spectral temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) enddo loadCases(currentLoadCase)%rotation = math_plain9to33(temp_valueVector) - case('totalstrain') - yieldStop = .True. - stopFlag = 'totalStrain' - yieldStopValue = IO_floatValue(line,chunkPos,i+1_pInt) - case('plasticstrain') - yieldStop = .True. - stopFlag = 'plasticStrain' - yieldStopValue = IO_floatValue(line,chunkPos,i+1_pInt) - case('plasticwork') - yieldStop = .True. - stopFlag = 'plasticWork' - yieldStopValue = IO_floatValue(line,chunkPos,i+1_pInt) end select enddo; enddo close(FILEUNIT) - if(yieldStop) then ! initialize variables related to yield stop - yieldStressNew = 0.0_pReal - plasticStrainNew = 0.0_pReal - eqStressNew = 0.0_pReal - eqTotalStrainNew = 0.0_pReal - eqPlasticStrainNew = 0.0_pReal - plasticWorkNew = 0.0_pReal - endif !-------------------------------------------------------------------------------------------------- ! consistency checks and output of load case loadCases(1)%followFormerTrajectory = .false. ! cannot guess along trajectory for first inc of first currentLoadCase @@ -387,7 +351,7 @@ program DAMASK_spectral if (loadCases(currentLoadCase)%outputfrequency < 1_pInt) errorID = 836_pInt ! non-positive result frequency write(6,'(2x,a,i5)') 'output frequency: ', & loadCases(currentLoadCase)%outputfrequency - write(6,'(2x,a,i5,/)') 'restart frequency: ', & + write(6,'(2x,a,i5,/)') 'restart frequency: ', & loadCases(currentLoadCase)%restartfrequency if (errorID > 0_pInt) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message enddo checkLoadcases @@ -478,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)), & @@ -489,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 @@ -519,37 +483,43 @@ program DAMASK_spectral ( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc,pReal)/& real(loadCases(currentLoadCase)%incs ,pReal))& -(1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( (inc-1_pInt),pReal)/& - real(loadCases(currentLoadCase)%incs ,pReal))) + 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)//& - ',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//& + ',a,'//IO_intOut(inc) //',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//& + ',a,'//IO_intOut(stepFraction) //',a,'//IO_intOut(subStepFactor**cutBackLevel)//& ',a,'//IO_intOut(currentLoadCase)//',a,'//IO_intOut(size(loadCases))//')') & 'Time', time, & 's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,& '-', stepFraction, '/', subStepFactor**cutBackLevel,& ' of load case ', currentLoadCase,'/',size(loadCases) + 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) - 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 !-------------------------------------------------------------------------------------------------- ! forward fields @@ -578,7 +548,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 @@ -618,65 +588,63 @@ 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 ......................................' + 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) @@ -688,93 +656,29 @@ 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 - - yieldCheck: if(yieldStop) then ! check if it yields or satisfies the certain stop condition - yieldStressOld = yieldStressNew - plasticStrainOld = plasticStrainNew - eqStressOld = eqStressNew - eqTotalStrainOld = eqTotalStrainNew - eqPlasticStrainOld = eqPlasticStrainNew - plasticWorkOld = plasticWorkNew - - call utilities_calcPlasticity(yieldStressNew, plasticStrainNew, eqStressNew, eqTotalStrainNew, & - eqPlasticStrainNew, plasticWorkNew, loadCases(currentLoadCase)%rotation) - - if (worldrank == 0) then ! output the stress-strain curve to file if yield stop criterion is used - if ((currentLoadCase == 1_pInt) .and. (inc == 1_pInt)) then - open(newunit=stressstrainUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& - '.stressstrain',form='FORMATTED',status='REPLACE') - write(stressstrainUnit,*) 0.0_pReal, 0.0_pReal - write(stressstrainUnit,*) eqTotalStrainNew, eqStressNew - close(stressstrainUnit) - else - open(newunit=stressstrainUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& - '.stressstrain',form='FORMATTED', position='APPEND', status='OLD') - write(stressstrainUnit,*) eqTotalStrainNew, eqStressNew - close(stressstrainUnit) - endif - endif - - if(stopFlag == 'totalStrain') then - if(eqTotalStrainNew > yieldStopValue) then - yieldStress = yieldStressOld * (eqTotalStrainNew - yieldStopValue)/(eqTotalStrainNew - eqTotalStrainOld) & ! linear interpolation of stress values - + yieldStressNew * (yieldStopValue - eqTotalStrainOld)/(eqTotalStrainNew - eqTotalStrainOld) - plasticStrainRate = (plasticStrainNew - plasticStrainOld)/(time - time0) ! calculate plastic strain rate - yieldStopSatisfied = .True. - endif - elseif(stopFlag == 'plasticStrain') then - if(eqPlasticStrainNew > yieldStopValue) then - yieldStress = yieldStressOld * (eqPlasticStrainNew - yieldStopValue)/(eqPlasticStrainNew - eqPlasticStrainOld) & - + yieldStressNew * (yieldStopValue - eqPlasticStrainOld)/(eqPlasticStrainNew - eqPlasticStrainOld) - plasticStrainRate = (plasticStrainNew - plasticStrainOld)/(time - time0) - yieldStopSatisfied = .True. - endif - elseif(stopFlag == 'plasticWork') then - if(plasticWorkNew > yieldStopValue) then - yieldStress = yieldStressOld * (plasticWorkNew - yieldStopValue)/(plasticWorkNew - plasticWorkOld) & - + yieldStressNew * (yieldStopValue - plasticWorkOld)/(plasticWorkNew - plasticWorkOld) - plasticStrainRate = (plasticStrainNew - plasticStrainOld)/(time - time0) - yieldStopSatisfied = .True. - endif - endif - endif yieldCheck - if (yieldStopSatisfied) then ! when yield, write the yield stress and strain rate to file and quit the job - if (worldrank == 0) then - open(newunit=yieldResUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& - '.yield',form='FORMATTED',status='REPLACE') - do i = 1_pInt,3_pInt - write(yieldResUnit,*) (yieldStress(i,j), j=1,3) - enddo - do i = 1_pInt,3_pInt - write(yieldResUnit,*) (plasticStrainRate(i,j), j=1,3) - enddo - close(yieldResUnit) - call quit(0_pInt) - endif - endif + endif skipping enddo incLooping + enddo loadCaseLooping !-------------------------------------------------------------------------------------------------- ! report summary of whole calculation write(6,'(/,a)') ' ###########################################################################' - write(6,'(1x,i6.6,a,i6.6,a,f5.1,a)') convergedCounter, ' out of ', & - notConvergedCounter + convergedCounter, ' (', & - real(convergedCounter, pReal)/& - real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, & - ' %) increments converged!' + write(6,'(1x,'//IO_intOut(convergedCounter)//',a,'//IO_intOut(notConvergedCounter + convergedCounter)//',a,f5.1,a)') & + convergedCounter, ' out of ', & + notConvergedCounter + convergedCounter, ' (', & + real(convergedCounter, pReal)/& + real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, & + ' %) increments converged!' + flush(6) call MPI_file_close(resUnit,ierr) close(statUnit) @@ -831,8 +735,6 @@ subroutine quit(stop_id) call PETScFinalize(ierr) if (ierr /= 0) write(6,'(a)') ' Error in PETScFinalize' #ifdef _OPENMP - ! If openMP is enabled, MPI is initialized before and independently of PETSc. Hence, also - ! take care of the finalization call MPI_finalize(error) if (error /= 0) write(6,'(a)') ' Error in MPI_finalize' #endif diff --git a/src/IO.f90 b/src/IO.f90 index d024d04ff..9e8033f73 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -1268,11 +1268,11 @@ integer(pInt) function IO_countNumericalDataLines(fileUnit) line = IO_read(fileUnit) chunkPos = IO_stringPos(line) tmp = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) - if (verify(trim(tmp) ,"0123456789")/=0) then ! found keyword + if (verify(trim(tmp),'0123456789') == 0) then ! numerical values + IO_countNumericalDataLines = IO_countNumericalDataLines + 1_pInt + else line = IO_read(fileUnit, .true.) ! reset IO_read exit - else - IO_countNumericalDataLines = IO_countNumericalDataLines + 1_pInt endif enddo backspace(fileUnit) diff --git a/src/crystallite.f90 b/src/crystallite.f90 old mode 100644 new mode 100755 index dd6a01786..3ff4417fe --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -55,14 +55,14 @@ module crystallite crystallite_Li0, & !< intermediate velocitiy grad at start of FE inc crystallite_partionedLi0,& !< intermediate velocity grad at start of homog inc crystallite_Fe, & !< current "elastic" def grad (end of converged time step) - crystallite_P, & !< 1st Piola-Kirchhoff stress per grain - crystallite_subF !< def grad to be reached at end of crystallite inc + crystallite_P !< 1st Piola-Kirchhoff stress per grain real(pReal), dimension(:,:,:,:,:), allocatable, private :: & crystallite_subFe0,& !< "elastic" def grad at start of crystallite inc crystallite_invFp, & !< inverse of current plastic def grad (end of converged time step) crystallite_subFp0,& !< plastic def grad at start of crystallite inc crystallite_invFi, & !< inverse of current intermediate def grad (end of converged time step) crystallite_subFi0,& !< intermediate def grad at start of crystallite inc + crystallite_subF, & !< def grad to be reached at end of crystallite inc crystallite_subF0, & !< def grad at start of crystallite inc crystallite_subLp0,& !< plastic velocity grad at start of crystallite inc crystallite_subLi0,& !< intermediate velocity grad at start of crystallite inc @@ -193,7 +193,7 @@ subroutine crystallite_init c, & !< counter in integration point component loop i, & !< counter in integration point loop e, & !< counter in element loop - o, & !< counter in output loop + o = 0_pInt, & !< counter in output loop r, & !< counter in crystallite loop cMax, & !< maximum number of integration point components iMax, & !< maximum number of integration points @@ -1239,14 +1239,16 @@ subroutine crystallite_integrateStateRK4() use numerics, only: & numerics_integrationMode use debug, only: & +#ifdef DEBUG + debug_e, & + debug_i, & + debug_g, & +#endif debug_level, & debug_crystallite, & debug_levelBasic, & debug_levelExtensive, & debug_levelSelective, & - debug_e, & - debug_i, & - debug_g, & debug_StateLoopDistribution use FEsolving, only: & FEsolving_execElem, & @@ -1533,14 +1535,16 @@ subroutine crystallite_integrateStateRKCK45() use, intrinsic :: & IEEE_arithmetic use debug, only: & +#ifdef DEBUG + debug_e, & + debug_i, & + debug_g, & +#endif debug_level, & debug_crystallite, & debug_levelBasic, & debug_levelExtensive, & debug_levelSelective, & - debug_e, & - debug_i, & - debug_g, & debug_StateLoopDistribution use numerics, only: & rTol_crystalliteState, & @@ -2036,14 +2040,16 @@ subroutine crystallite_integrateStateAdaptiveEuler() use, intrinsic :: & IEEE_arithmetic use debug, only: & +#ifdef DEBUG + debug_e, & + debug_i, & + debug_g, & +#endif debug_level, & debug_crystallite, & debug_levelBasic, & debug_levelExtensive, & debug_levelSelective, & - debug_e, & - debug_i, & - debug_g, & debug_StateLoopDistribution use numerics, only: & rTol_crystalliteState, & @@ -2396,14 +2402,16 @@ subroutine crystallite_integrateStateEuler() use, intrinsic :: & IEEE_arithmetic use debug, only: & +#ifdef DEBUG + debug_e, & + debug_i, & + debug_g, & +#endif debug_level, & debug_crystallite, & debug_levelBasic, & debug_levelExtensive, & debug_levelSelective, & - debug_e, & - debug_i, & - debug_g, & debug_StateLoopDistribution use numerics, only: & numerics_integrationMode, & @@ -2619,9 +2627,11 @@ subroutine crystallite_integrateStateFPI() use, intrinsic :: & IEEE_arithmetic use debug, only: & +#ifdef DEBUG debug_e, & debug_i, & debug_g, & +#endif debug_level,& debug_crystallite, & debug_levelBasic, & @@ -3068,14 +3078,16 @@ logical function crystallite_stateJump(ipc,ip,el) IEEE_arithmetic use prec, only: & dNeq0 +#ifdef DEBUG use debug, only: & + debug_e, & + debug_i, & + debug_g, & debug_level, & debug_crystallite, & debug_levelExtensive, & - debug_levelSelective, & - debug_e, & - debug_i, & - debug_g + debug_levelSelective +#endif use material, only: & plasticState, & sourceState, & @@ -3205,9 +3217,11 @@ logical function crystallite_integrateStress(& debug_levelBasic, & debug_levelExtensive, & debug_levelSelective, & +#ifdef DEBUG debug_e, & debug_i, & debug_g, & +#endif debug_cumLpCalls, & debug_cumLpTicks, & debug_StressLoopLpDistribution, & @@ -3233,7 +3247,9 @@ logical function crystallite_integrateStress(& math_Plain33to9, & math_Plain9to33, & math_Plain99to3333 +#ifdef DEBUG use mesh, only: mesh_element +#endif implicit none integer(pInt), intent(in):: el, & ! element index @@ -3297,8 +3313,8 @@ logical function crystallite_integrateStress(& p, & jacoCounterLp, & jacoCounterLi ! counters to check for Jacobian update - integer(pLongInt) tick, & - tock, & + integer(pLongInt) :: tick = 0_pLongInt, & + tock = 0_pLongInt, & tickrate, & maxticks diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 7cc5fbd04..7dbea41d5 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -16,7 +16,7 @@ module homogenization ! General variables for the homogenization at a material point implicit none private - real(pReal), dimension(:,:,:,:), allocatable, public :: & + real(pReal), dimension(:,:,:,:), allocatable, public :: & materialpoint_F0, & !< def grad of IP at start of FE increment materialpoint_F, & !< def grad of IP to be reached at end of FE increment materialpoint_P !< first P--K stress of IP @@ -128,7 +128,7 @@ subroutine homogenization_init integer(pInt), dimension(:) , pointer :: thisNoutput character(len=64), dimension(:,:), pointer :: thisOutput character(len=32) :: outputName !< name of output, intermediate fix until HDF5 output is ready - logical :: knownHomogenization, knownThermal, knownDamage, knownVacancyflux, knownPorosity, knownHydrogenflux + logical :: valid !-------------------------------------------------------------------------------------------------- @@ -199,7 +199,7 @@ subroutine homogenization_init do p = 1,material_Nhomogenization if (any(material_homog == p)) then i = homogenization_typeInstance(p) ! which instance of this homogenization type - knownHomogenization = .true. ! assume valid + valid = .true. ! assume valid select case(homogenization_type(p)) ! split per homogenization type case (HOMOGENIZATION_NONE_ID) outputName = HOMOGENIZATION_NONE_label @@ -217,10 +217,10 @@ subroutine homogenization_init thisOutput => homogenization_RGC_output thisSize => homogenization_RGC_sizePostResult case default - knownHomogenization = .false. + valid = .false. end select write(FILEUNIT,'(/,a,/)') '['//trim(homogenization_name(p))//']' - if (knownHomogenization) then + if (valid) then write(FILEUNIT,'(a)') '(type)'//char(9)//trim(outputName) write(FILEUNIT,'(a,i4)') '(ngrains)'//char(9),homogenization_Ngrains(p) if (homogenization_type(p) /= HOMOGENIZATION_NONE_ID) then @@ -230,7 +230,7 @@ subroutine homogenization_init endif endif i = thermal_typeInstance(p) ! which instance of this thermal type - knownThermal = .true. ! assume valid + valid = .true. ! assume valid select case(thermal_type(p)) ! split per thermal type case (THERMAL_isothermal_ID) outputName = THERMAL_isothermal_label @@ -248,9 +248,9 @@ subroutine homogenization_init thisOutput => thermal_conduction_output thisSize => thermal_conduction_sizePostResult case default - knownThermal = .false. + valid = .false. end select - if (knownThermal) then + if (valid) then write(FILEUNIT,'(a)') '(thermal)'//char(9)//trim(outputName) if (thermal_type(p) /= THERMAL_isothermal_ID) then do e = 1,thisNoutput(i) @@ -259,7 +259,7 @@ subroutine homogenization_init endif endif i = damage_typeInstance(p) ! which instance of this damage type - knownDamage = .true. ! assume valid + valid = .true. ! assume valid select case(damage_type(p)) ! split per damage type case (DAMAGE_none_ID) outputName = DAMAGE_none_label @@ -277,9 +277,9 @@ subroutine homogenization_init thisOutput => damage_nonlocal_output thisSize => damage_nonlocal_sizePostResult case default - knownDamage = .false. + valid = .false. end select - if (knownDamage) then + if (valid) then write(FILEUNIT,'(a)') '(damage)'//char(9)//trim(outputName) if (damage_type(p) /= DAMAGE_none_ID) then do e = 1,thisNoutput(i) @@ -288,7 +288,7 @@ subroutine homogenization_init endif endif i = vacancyflux_typeInstance(p) ! which instance of this vacancy flux type - knownVacancyflux = .true. ! assume valid + valid = .true. ! assume valid select case(vacancyflux_type(p)) ! split per vacancy flux type case (VACANCYFLUX_isoconc_ID) outputName = VACANCYFLUX_isoconc_label @@ -306,9 +306,9 @@ subroutine homogenization_init thisOutput => vacancyflux_cahnhilliard_output thisSize => vacancyflux_cahnhilliard_sizePostResult case default - knownVacancyflux = .false. + valid = .false. end select - if (knownVacancyflux) then + if (valid) then write(FILEUNIT,'(a)') '(vacancyflux)'//char(9)//trim(outputName) if (vacancyflux_type(p) /= VACANCYFLUX_isoconc_ID) then do e = 1,thisNoutput(i) @@ -317,7 +317,7 @@ subroutine homogenization_init endif endif i = porosity_typeInstance(p) ! which instance of this porosity type - knownPorosity = .true. ! assume valid + valid = .true. ! assume valid select case(porosity_type(p)) ! split per porosity type case (POROSITY_none_ID) outputName = POROSITY_none_label @@ -330,9 +330,9 @@ subroutine homogenization_init thisOutput => porosity_phasefield_output thisSize => porosity_phasefield_sizePostResult case default - knownPorosity = .false. + valid = .false. end select - if (knownPorosity) then + if (valid) then write(FILEUNIT,'(a)') '(porosity)'//char(9)//trim(outputName) if (porosity_type(p) /= POROSITY_none_ID) then do e = 1,thisNoutput(i) @@ -341,7 +341,7 @@ subroutine homogenization_init endif endif i = hydrogenflux_typeInstance(p) ! which instance of this hydrogen flux type - knownHydrogenflux = .true. ! assume valid + valid = .true. ! assume valid select case(hydrogenflux_type(p)) ! split per hydrogen flux type case (HYDROGENFLUX_isoconc_ID) outputName = HYDROGENFLUX_isoconc_label @@ -354,9 +354,9 @@ subroutine homogenization_init thisOutput => hydrogenflux_cahnhilliard_output thisSize => hydrogenflux_cahnhilliard_sizePostResult case default - knownHydrogenflux = .false. + valid = .false. end select - if (knownHydrogenflux) then + if (valid) then write(FILEUNIT,'(a)') '(hydrogenflux)'//char(9)//trim(outputName) if (hydrogenflux_type(p) /= HYDROGENFLUX_isoconc_ID) then do e = 1,thisNoutput(i) diff --git a/src/numerics.f90 b/src/numerics.f90 index 00d97297e..c854d9d2b 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 f2b9d3f7f..176ac7194 100644 --- a/src/spectral_mech_AL.f90 +++ b/src/spectral_mech_AL.f90 @@ -43,17 +43,20 @@ module spectral_mech_AL !-------------------------------------------------------------------------------------------------- ! stress, stiffness and compliance average etc. real(pReal), private, dimension(3,3) :: & - F_aimDot, & !< assumed rate of average deformation gradient + F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient F_aim = math_I3, & !< current prescribed deformation gradient F_aim_lastInc = math_I3, & !< previous average deformation gradient F_av = 0.0_pReal, & !< average incompatible def grad field P_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress P_avLastEval = 0.0_pReal !< average 1st Piola--Kirchhoff stress last call of CPFEM_general - character(len=1024), private :: incInfo !< time and increment information + + character(len=1024), private :: incInfo !< time and increment information + 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 + C_minMaxAvgLastInc = 0.0_pReal, & !< previous (min+max)/2 stiffness S = 0.0_pReal, & !< current compliance (filled up with zeros) C_scale = 0.0_pReal, & S_scale = 0.0_pReal @@ -62,7 +65,7 @@ module spectral_mech_AL err_BC, & !< deviation from stress BC err_curl, & !< RMS of curl of F err_div !< RMS of div of P - logical, private :: ForwardData + integer(pInt), private :: & totalIter = 0_pInt !< total iteration in current increment @@ -80,7 +83,7 @@ module spectral_mech_AL 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 !> @todo use sourced allocation, e.g. allocate(Fdot,source = F_lastInc) !-------------------------------------------------------------------------------------------------- subroutine AL_init @@ -102,12 +105,15 @@ subroutine AL_init use numerics, only: & worldrank, & worldsize + use homogenization, only: & + materialpoint_F0 use DAMASK_interface, only: & getSolverJobName use spectral_utilities, only: & Utilities_constitutiveResponse, & Utilities_updateGamma, & - Utilities_updateIPcoords + Utilities_updateIPcoords, & + wgt use mesh, only: & grid, & grid3 @@ -120,7 +126,11 @@ subroutine AL_init temp33_Real = 0.0_pReal PetscErrorCode :: ierr - PetscScalar, pointer, dimension(:,:,:,:) :: xx_psc, F, F_lambda + PetscScalar, pointer, dimension(:,:,:,:) :: & + FandF_lambda, & ! overall pointer to solution data + F, & ! specific (sub)pointer + F_lambda ! specific (sub)pointer + integer(pInt), dimension(:), allocatable :: localK integer(pInt) :: proc character(len=1024) :: rankStr @@ -160,51 +170,43 @@ subroutine AL_init DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point grid(1),grid(2),grid(3), & ! global grid 1 , 1, worldsize, & - 18, 0, & ! #dof (F tensor), ghost boundary width (domain overlap) + 18, 0, & ! #dof (F tensor), ghost boundary width (domain overlap) grid(1),grid(2),localK, & ! local grid da,ierr) ! handle, error CHKERRQ(ierr) - call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) - call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) - call DMDASNESSetFunctionLocal(da,INSERT_VALUES,AL_formResidual,PETSC_NULL_OBJECT,ierr) + call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) ! connect snes to da + call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 9, i.e. every def grad tensor) + call DMDASNESSetFunctionLocal(da,INSERT_VALUES,AL_formResidual,PETSC_NULL_OBJECT,ierr) ! residual vector of same shape as solution vector + CHKERRQ(ierr) + call SNESSetConvergenceTest(snes,AL_converged,PETSC_NULL_OBJECT,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged" CHKERRQ(ierr) - call SNESSetConvergenceTest(snes,AL_converged,PETSC_NULL_OBJECT,PETSC_NULL_FUNCTION,ierr) - CHKERRQ(ierr) - call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr) + call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional cli arguments !-------------------------------------------------------------------------------------------------- ! init fields - call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! places pointer xx_psc on PETSc data - F => xx_psc(0:8,:,:,:) - F_lambda => xx_psc(9:17,:,:,:) + call DMDAVecGetArrayF90(da,solution_vec,FandF_lambda,ierr); CHKERRQ(ierr) ! places pointer on PETSc data + F => FandF_lambda( 0: 8,:,:,:) + F_lambda => FandF_lambda( 9:17,:,:,:) + restart: if (restartInc > 1_pInt) then - if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) & + if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0) then write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & - 'reading values of increment ', restartInc - 1_pInt, ' from file' - flush(6) + 'reading values of increment ', restartInc-1_pInt, ' from file' + flush(6) + endif 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_lambda'//trim(rankStr),trim(getSolverJobName()),size(F_lambda)) - read (777,rec=1) F_lambda - close (777) - call IO_read_realFile(777,'F_lambda_lastInc'//trim(rankStr),& - trim(getSolverJobName()),size(F_lambda_lastInc)) - read (777,rec=1) F_lambda_lastInc - close (777) - call IO_read_realFile(777,'F_aim', trim(getSolverJobName()),size(F_aim)) - read (777,rec=1) F_aim - close (777) - call IO_read_realFile(777,'F_aim_lastInc', trim(getSolverJobName()),size(F_aim_lastInc)) - read (777,rec=1) F_aim_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_lambda; close (777) + call IO_read_realFile(777,'F_lambda_lastInc'//trim(rankStr),trim(getSolverJobName()),size(F_lambda_lastInc)) + read (777,rec=1) F_lambda_lastInc; close (777) + call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(F_aimDot)) + 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 F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) @@ -212,30 +214,30 @@ subroutine AL_init F_lambda_lastInc = F_lastInc 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,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 nullify(F) nullify(F_lambda) - call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! write data back to PETSc + call DMDAVecRestoreArrayF90(da,solution_vec,FandF_lambda,ierr); CHKERRQ(ierr) ! write data back to PETSc restartRead: if (restartInc > 1_pInt) then 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.) C_scale = C_minMaxAvg S_scale = math_invSym3333(C_minMaxAvg) @@ -245,8 +247,7 @@ end subroutine AL_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the AL scheme with internal iterations !-------------------------------------------------------------------------------------------------- -type(tSolutionState) function & - AL_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) +type(tSolutionState) function AL_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) use IO, only: & IO_error use numerics, only: & @@ -265,13 +266,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 !-------------------------------------------------------------------------------------------------- @@ -304,18 +305,17 @@ type(tSolutionState) function & !-------------------------------------------------------------------------------------------------- ! solve BVP - call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr) - CHKERRQ(ierr) + call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr); CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! check convergence - call SNESGetConvergedReason(snes,reason,ierr) - CHKERRQ(ierr) + call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr) + + AL_solution%converged = reason > 0 + AL_solution%iterationsNeeded = totalIter AL_solution%termIll = terminallyIll terminallyIll = .false. - if (reason == -4) call IO_error(893_pInt) - if (reason < 1) AL_solution%converged = .false. - AL_solution%iterationsNeeded = totalIter + if (reason == -4) call IO_error(893_pInt) ! MPI error end function AL_solution @@ -330,8 +330,8 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) polarAlpha, & polarBeta use mesh, only: & - grid3, & - grid + grid, & + grid3 use IO, only: & IO_intOut use math, only: & @@ -340,6 +340,10 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) math_mul3333xx33, & math_invSym3333, & math_mul33x33 + use debug, only: & + debug_level, & + debug_spectral, & + debug_spectralRotation use spectral_utilities, only: & wgt, & tensorField_real, & @@ -349,27 +353,17 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) Utilities_constitutiveResponse, & Utilities_divergenceRMS, & Utilities_curlRMS - use debug, only: & - debug_level, & - debug_spectral, & - debug_spectralRotation use homogenization, only: & materialpoint_dPdF use FEsolving, only: & terminallyIll implicit none -!-------------------------------------------------------------------------------------------------- -! strange syntax in the next line because otherwise macros expand beyond 132 character limit - 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 + DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in + PetscScalar, & + target, dimension(3,3,2, XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: x_scal !< what is this? + PetscScalar, & + target, dimension(3,3,2, X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: f_scal !< what is this? PetscScalar, pointer, dimension(:,:,:,:,:) :: & F, & F_lambda, & @@ -387,33 +381,33 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) SNESGetNumberFunctionEvals, & SNESGetIterationNumber - F => x_scal(1:3,1:3,1,& - XG_RANGE,YG_RANGE,ZG_RANGE) - F_lambda => x_scal(1:3,1:3,2,& - XG_RANGE,YG_RANGE,ZG_RANGE) - residual_F => f_scal(1:3,1:3,1,& - X_RANGE,Y_RANGE,Z_RANGE) + F => x_scal(1:3,1:3,1,& + XG_RANGE,YG_RANGE,ZG_RANGE) + F_lambda => x_scal(1:3,1:3,2,& + XG_RANGE,YG_RANGE,ZG_RANGE) + residual_F => f_scal(1:3,1:3,1,& + X_RANGE, Y_RANGE, Z_RANGE) residual_F_lambda => f_scal(1:3,1:3,2,& - X_RANGE,Y_RANGE,Z_RANGE) - - call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) - call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) + X_RANGE, Y_RANGE, Z_RANGE) F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - if(nfuncs== 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment - newIteration: if(totalIter <= PETScIter) then + 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 !-------------------------------------------------------------------------------------------------- -! 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 @@ -425,7 +419,6 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& polarAlpha*math_mul33x33(F(1:3,1:3,i,j,k), & math_mul3333xx33(C_scale,F_lambda(1:3,1:3,i,j,k) - math_I3)) - enddo; enddo; enddo !-------------------------------------------------------------------------------------------------- @@ -435,24 +428,23 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) call utilities_FFTtensorBackward() !-------------------------------------------------------------------------------------------------- -! constructing residual - residual_F_lambda = polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) +! constructing F_lambda residual + residual_F_lambda = polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) !< eq (16) in doi: 10.1016/j.ijplas.2014.02.006 !-------------------------------------------------------------------------------------------------- ! 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. !-------------------------------------------------------------------------------------------------- ! calculate divergence tensorField_real = 0.0_pReal - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F + tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F !< stress field in disguise call utilities_FFTtensorForward() - err_div = Utilities_divergenceRMS() - call utilities_FFTtensorBackward() + err_div = Utilities_divergenceRMS() !< root mean squared error in divergence of stress !-------------------------------------------------------------------------------------------------- ! constructing residual @@ -463,7 +455,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) residual_F(1:3,1:3,i,j,k) - & math_mul33x33(F(1:3,1:3,i,j,k), & math_mul3333xx33(C_scale,F_lambda(1:3,1:3,i,j,k) - math_I3))) & - + residual_F_lambda(1:3,1:3,i,j,k) + + residual_F_lambda(1:3,1:3,i,j,k) !< eq (16) in doi: 10.1016/j.ijplas.2014.02.006 enddo; enddo; enddo !-------------------------------------------------------------------------------------------------- @@ -472,8 +464,11 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F call utilities_FFTtensorForward() err_curl = Utilities_curlRMS() - call utilities_FFTtensorBackward() - + + nullify(F) + nullify(F_lambda) + nullify(residual_F) + nullify(residual_F_lambda) end subroutine AL_formResidual @@ -488,8 +483,8 @@ subroutine AL_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr err_div_tolAbs, & err_curl_tolRel, & err_curl_tolAbs, & - err_stress_tolAbs, & - err_stress_tolRel + err_stress_tolRel, & + err_stress_tolAbs use math, only: & math_mul3333xx33 use FEsolving, only: & @@ -508,24 +503,24 @@ subroutine AL_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr real(pReal) :: & curlTol, & divTol, & - BC_tol + BCTol !-------------------------------------------------------------------------------------------------- ! stress BC handling F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc - err_BC = maxval(abs((-mask_stress+1.0_pReal)*math_mul3333xx33(C_scale,F_aim-F_av) + & - mask_stress *(P_av - params%stress_BC))) ! mask = 0.0 for no bc + err_BC = maxval(abs((1.0_pReal-mask_stress) * math_mul3333xx33(C_scale,F_aim-F_av) + & + mask_stress * (P_av-params%stress_BC))) ! mask = 0.0 for no bc !-------------------------------------------------------------------------------------------------- ! error calculation - curlTol = max(maxval(abs(F_aim-math_I3))*err_curl_tolRel,err_curl_tolAbs) - divTol = max(maxval(abs(P_av)) *err_div_tolRel,err_div_tolAbs) - BC_tol = max(maxval(abs(P_av)) *err_stress_tolrel,err_stress_tolabs) + curlTol = max(maxval(abs(F_aim-math_I3))*err_curl_tolRel ,err_curl_tolAbs) + divTol = max(maxval(abs(P_av)) *err_div_tolRel ,err_div_tolAbs) + BCTol = max(maxval(abs(P_av)) *err_stress_tolRel,err_stress_tolAbs) converged: if ((totalIter >= itmin .and. & - all([ err_div/divTol, & + all([ err_div /divTol, & err_curl/curlTol, & - err_BC/BC_tol ] < 1.0_pReal)) & + err_BC /BCTol ] < 1.0_pReal)) & .or. terminallyIll) then reason = 1 elseif (totalIter >= itmax) then converged @@ -537,12 +532,12 @@ subroutine AL_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr !-------------------------------------------------------------------------------------------------- ! report write(6,'(1/,a)') ' ... reporting .............................................................' - write(6,'(/,a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', & - err_curl/curlTol,' (',err_curl,' -, tol =',curlTol,')' - write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & - err_div/divTol, ' (',err_div, ' / m, tol =',divTol,')' - write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', & - err_BC/BC_tol, ' (',err_BC, ' Pa, tol =',BC_tol,')' + write(6,'(/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & + err_div/divTol, ' (',err_div, ' / m, tol = ',divTol,')' + write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', & + err_curl/curlTol,' (',err_curl,' -, tol = ',curlTol,')' + write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', & + err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' write(6,'(/,a)') ' ===========================================================================' flush(6) @@ -550,130 +545,142 @@ end subroutine AL_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 AL_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC) - use math, only: & - math_mul33x33, & - math_mul3333xx33, & - math_transpose33, & - math_rotate_backward33 - use numerics, only: & - worldrank - use mesh, only: & - grid3, & - grid - 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_mul3333xx33, & + math_transpose33, & + 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, dimension(:,:,:,:), pointer :: xx_psc, F, F_lambda - integer(pInt) :: i, j, k - real(pReal), dimension(3,3) :: F_lambda33 - character(len=1024) :: rankStr + 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, dimension(:,:,:,:), pointer :: FandF_lambda, F, F_lambda + integer(pInt) :: i, j, k + real(pReal), dimension(3,3) :: F_lambda33 + character(len=32) :: rankStr !-------------------------------------------------------------------------------------------------- ! update coordinates and rate and forward last inc - call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr) - F => xx_psc(0:8,:,:,:) - F_lambda => xx_psc(9:17,:,:,:) - 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) - call IO_write_jobRealFile(777,'F_lambda'//trim(rankStr),size(F_lambda)) ! writing deformation gradient field to file - write (777,rec=1) F_lambda - close (777) - call IO_write_jobRealFile(777,'F_lambda_lastInc'//trim(rankStr),size(F_lambda_lastInc)) ! writing F_lastInc field to file - write (777,rec=1) F_lambda_lastInc - close (777) - if (worldrank == 0_pInt) then - call IO_write_jobRealFile(777,'F_aim',size(F_aim)) - write (777,rec=1) F_aim - close(777) - call IO_write_jobRealFile(777,'F_aim_lastInc',size(F_aim_lastInc)) - write (777,rec=1) F_aim_lastInc - close(777) - 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,FandF_lambda,ierr); CHKERRQ(ierr) + F => FandF_lambda( 0: 8,:,:,:) + F_lambda => FandF_lambda( 9:17,:,:,:) - call utilities_updateIPcoords(F) + 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) - if (cutBack) then - F_aim = F_aim_lastInc - F_lambda = reshape(F_lambda_lastInc,[9,grid(1),grid(2),grid3]) - 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 + 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) + call IO_write_jobRealFile(777,'F_aimDot',size(F_aimDot)) + write (777,rec=1) F_aimDot; close(777) + endif + + 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) + call IO_write_jobRealFile(777,'F_lambda'//trim(rankStr),size(F_lambda)) ! writing deformation gradient field to file + write (777,rec=1) F_lambda; close (777) + call IO_write_jobRealFile(777,'F_lambda_lastInc'//trim(rankStr),size(F_lambda_lastInc)) ! writing F_lastInc field to file + write (777,rec=1) F_lambda_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_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]) ! winding F forward + F_lambda_lastInc = reshape(F_lambda, [3,3,grid(1),grid(2),grid3]) ! winding F_lambda 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_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])) - 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 +! update average and local deformation gradients + F_aim = F_aim_lastInc + F_aimDot * timeinc - 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]) - F_lambda = reshape(Utilities_forwardField(timeinc,F_lambda_lastInc,F_lambdadot), & - [9,grid(1),grid(2),grid3]) ! does not have any average value as boundary condition - if (.not. guess) then ! large strain forwarding - do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) + 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]) + if (guess) then + F_lambda = reshape(Utilities_forwardField(timeinc,F_lambda_lastInc,F_lambdadot), & + [9,grid(1),grid(2),grid3]) ! does not have any average value as boundary condition + else + do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) F_lambda33 = reshape(F_lambda(1:9,i,j,k),[3,3]) F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, & math_mul3333xx33(C_scale,& @@ -681,9 +688,12 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stre F_lambda33) -math_I3))*0.5_pReal)& + math_I3 F_lambda(1:9,i,j,k) = reshape(F_lambda33,[9]) - enddo; enddo; enddo - endif - call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) + enddo; enddo; enddo + endif + + nullify(F) + nullify(F_lambda) + call DMDAVecRestoreArrayF90(da,solution_vec,FandF_lambda,ierr); CHKERRQ(ierr) end subroutine AL_forward diff --git a/src/spectral_mech_Basic.f90 b/src/spectral_mech_Basic.f90 index 00502b6b5..44e043575 100644 --- a/src/spectral_mech_Basic.f90 +++ b/src/spectral_mech_Basic.f90 @@ -38,21 +38,27 @@ 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 - character(len=1024), private :: incInfo + F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient + F_aim = math_I3, & !< current prescribed deformation gradient + F_aim_lastInc = math_I3, & !< previous average deformation gradient + P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress + + character(len=1024), private :: incInfo !< time and increment information + 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) - real(pReal), private :: err_stress, err_div - logical, private :: ForwardData + C_minMaxAvgLastInc = 0.0_pReal, & !< previous (min+max)/2 stiffness + S = 0.0_pReal !< current compliance (filled up with zeros) + + real(pReal), private :: & + err_BC, & !< deviation from stress BC + err_div !< RMS of div of P + integer(pInt), private :: & totalIter = 0_pInt !< total iteration in current increment + real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal public :: & @@ -69,7 +75,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 #if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 @@ -90,6 +96,8 @@ subroutine basicPETSc_init use numerics, only: & worldrank, & worldsize + use homogenization, only: & + materialpoint_F0 use DAMASK_interface, only: & getSolverJobName use spectral_utilities, only: & @@ -132,8 +140,8 @@ subroutine basicPETSc_init !-------------------------------------------------------------------------------------------------- ! allocate global fields - allocate (F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - allocate (Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate (F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate (Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc @@ -152,11 +160,10 @@ subroutine basicPETSc_init grid(1),grid(2),localK, & ! local grid da,ierr) ! handle, error CHKERRQ(ierr) - call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) + call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) ! connect snes to da call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 9, i.e. every def grad tensor) call DMDASNESSetFunctionLocal(da,INSERT_VALUES,BasicPETSC_formResidual,PETSC_NULL_OBJECT,ierr) ! residual vector of same shape as solution vector CHKERRQ(ierr) - call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) ! connect snes to da call SNESSetConvergenceTest(snes,BasicPETSC_converged,PETSC_NULL_OBJECT,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged" CHKERRQ(ierr) call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional cli arguments @@ -166,62 +173,55 @@ subroutine basicPETSc_init call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with restart: if (restartInc > 1_pInt) then - if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) & + if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0) then write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & - 'reading values of increment ', restartInc - 1_pInt, ' from file' - flush(6) + 'reading values of increment ', restartInc-1_pInt, ' from file' + flush(6) + endif 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) - 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_lastInc; close (777) + call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(F_aimDot)) + 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 + elseif (restartInc == 1_pInt) then restart F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity 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 +238,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 !-------------------------------------------------------------------------------------------------- @@ -261,7 +261,7 @@ type(tSolutionState) function & !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) - if (update_gamma) call Utilities_updateGamma(C_minmaxAvg,restartWrite) + if (update_gamma) call Utilities_updateGamma(C_minMaxAvg,restartWrite) !-------------------------------------------------------------------------------------------------- @@ -274,19 +274,17 @@ type(tSolutionState) function & !-------------------------------------------------------------------------------------------------- ! solve BVP - call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr) - CHKERRQ(ierr) + call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr); CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! 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 +320,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 +340,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_BC = 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 @@ -413,14 +409,14 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du PetscErrorCode :: ierr real(pReal) :: & divTol, & - stressTol + BCTol - divTol = max(maxval(abs(P_av))*err_div_tolRel,err_div_tolAbs) - stressTol = max(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs) + divTol = max(maxval(abs(P_av))*err_div_tolRel ,err_div_tolAbs) + BCTol = max(maxval(abs(P_av))*err_stress_tolRel,err_stress_tolAbs) converged: if ((totalIter >= itmin .and. & all([ err_div/divTol, & - err_stress/stressTol ] < 1.0_pReal)) & + err_BC /BCTol ] < 1.0_pReal)) & .or. terminallyIll) then reason = 1 elseif (totalIter >= itmax) then converged @@ -433,9 +429,9 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du ! report write(6,'(1/,a)') ' ... reporting .............................................................' write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & - err_div/divTol, ' (',err_div,' / m, tol =',divTol,')' - write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', & - err_stress/stressTol, ' (',err_stress, ' Pa, tol =',stressTol,')' + err_div/divTol, ' (',err_div,' / m, tol = ',divTol,')' + write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', & + err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' write(6,'(/,a)') ' ===========================================================================' flush(6) @@ -443,106 +439,118 @@ 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, dimension(:,:,:,:), 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,'F_aimDot',size(F_aimDot)) + write (777,rec=1) F_aimDot; 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 ! components of deformation_BC%maskFloat always start out with zero + 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 2d261ff1b..1a89718da 100644 --- a/src/spectral_mech_Polarisation.f90 +++ b/src/spectral_mech_Polarisation.f90 @@ -43,7 +43,7 @@ module spectral_mech_Polarisation !-------------------------------------------------------------------------------------------------- ! stress, stiffness and compliance average etc. real(pReal), private, dimension(3,3) :: & - F_aimDot, & !< assumed rate of average deformation gradient + F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient F_aim = math_I3, & !< current prescribed deformation gradient F_aim_lastInc = math_I3, & !< previous average deformation gradient F_av = 0.0_pReal, & !< average incompatible def grad field @@ -54,6 +54,7 @@ module spectral_mech_Polarisation 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 + C_minMaxAvgLastInc = 0.0_pReal, & !< previous (min+max)/2 stiffness S = 0.0_pReal, & !< current compliance (filled up with zeros) C_scale = 0.0_pReal, & S_scale = 0.0_pReal @@ -62,7 +63,7 @@ module spectral_mech_Polarisation err_BC, & !< deviation from stress BC err_curl, & !< RMS of curl of F err_div !< RMS of div of P - logical, private :: ForwardData + integer(pInt), private :: & totalIter = 0_pInt !< total iteration in current increment @@ -80,7 +81,7 @@ module spectral_mech_Polarisation 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 !> @todo use sourced allocation, e.g. allocate(Fdot,source = F_lastInc) !-------------------------------------------------------------------------------------------------- subroutine Polarisation_init @@ -102,12 +103,15 @@ subroutine Polarisation_init use numerics, only: & worldrank, & worldsize + use homogenization, only: & + materialpoint_F0 use DAMASK_interface, only: & getSolverJobName use spectral_utilities, only: & Utilities_constitutiveResponse, & Utilities_updateGamma, & - Utilities_updateIPcoords + Utilities_updateIPcoords, & + wgt use mesh, only: & grid, & grid3 @@ -120,7 +124,11 @@ subroutine Polarisation_init temp33_Real = 0.0_pReal PetscErrorCode :: ierr - PetscScalar, pointer, dimension(:,:,:,:) :: xx_psc, F, F_tau + PetscScalar, pointer, dimension(:,:,:,:) :: & + FandF_tau, & ! overall pointer to solution data + F, & ! specific (sub)pointer + F_tau ! specific (sub)pointer + integer(pInt), dimension(:), allocatable :: localK integer(pInt) :: proc character(len=1024) :: rankStr @@ -164,78 +172,69 @@ subroutine Polarisation_init grid(1),grid(2),localK, & ! local grid da,ierr) ! handle, error CHKERRQ(ierr) - call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) - call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) - call DMDASNESSetFunctionLocal(da,INSERT_VALUES,Polarisation_formResidual,PETSC_NULL_OBJECT,ierr) + call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) ! connect snes to da + call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 9, i.e. every def grad tensor) + call DMDASNESSetFunctionLocal(da,INSERT_VALUES,Polarisation_formResidual,PETSC_NULL_OBJECT,ierr) ! residual vector of same shape as solution vector + CHKERRQ(ierr) + call SNESSetConvergenceTest(snes,Polarisation_converged,PETSC_NULL_OBJECT,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged" CHKERRQ(ierr) - call SNESSetConvergenceTest(snes,Polarisation_converged,PETSC_NULL_OBJECT,PETSC_NULL_FUNCTION,ierr) - CHKERRQ(ierr) - call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr) + call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional cli arguments !-------------------------------------------------------------------------------------------------- ! init fields - call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! places pointer xx_psc on PETSc data - F => xx_psc(0:8,:,:,:) - F_tau => xx_psc(9:17,:,:,:) - restart: if (restartInc > 1_pInt) then - if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) & + call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! places pointer on PETSc data + F => FandF_tau( 0: 8,:,:,:) + F_tau => FandF_tau( 9:17,:,:,:) + restart: if (restartInc > 1_pInt) then + if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0) then write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & 'reading values of increment ', restartInc - 1_pInt, ' from file' - flush(6) + flush(6) + endif 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_tau'//trim(rankStr),trim(getSolverJobName()),size(F_tau)) - read (777,rec=1) F_tau - close (777) - call IO_read_realFile(777,'F_tau_lastInc'//trim(rankStr),& - trim(getSolverJobName()),size(F_tau_lastInc)) - read (777,rec=1) F_tau_lastInc - close (777) - call IO_read_realFile(777,'F_aim', trim(getSolverJobName()),size(F_aim)) - read (777,rec=1) F_aim - close (777) - call IO_read_realFile(777,'F_aim_lastInc', trim(getSolverJobName()),size(F_aim_lastInc)) - read (777,rec=1) F_aim_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_tau; close (777) + call IO_read_realFile(777,'F_tau_lastInc'//trim(rankStr),trim(getSolverJobName()),size(F_tau_lastInc)) + read (777,rec=1) F_tau_lastInc; close (777) + call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(F_aimDot)) + 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 - F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity + F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) - F_tau = 2.0_pReal* F + F_tau = 2.0_pReal*F F_tau_lastInc = 2.0_pReal*F_lastInc 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,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 nullify(F) nullify(F_tau) - call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! write data back to PETSc + call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! write data back to PETSc restartRead: if (restartInc > 1_pInt) then 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.) C_scale = C_minMaxAvg S_scale = math_invSym3333(C_minMaxAvg) @@ -245,8 +244,7 @@ end subroutine Polarisation_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the Polarisation scheme with internal iterations !-------------------------------------------------------------------------------------------------- -type(tSolutionState) function & - Polarisation_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) +type(tSolutionState) function Polarisation_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) use IO, only: & IO_error use numerics, only: & @@ -265,13 +263,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 !-------------------------------------------------------------------------------------------------- @@ -304,18 +302,17 @@ type(tSolutionState) function & !-------------------------------------------------------------------------------------------------- ! solve BVP - call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr) - CHKERRQ(ierr) + call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr); CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! check convergence - call SNESGetConvergedReason(snes,reason,ierr) - CHKERRQ(ierr) + call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr) + + Polarisation_solution%converged = reason > 0 + Polarisation_solution%iterationsNeeded = totalIter Polarisation_solution%termIll = terminallyIll terminallyIll = .false. - if (reason == -4) call IO_error(893_pInt) - if (reason < 1) Polarisation_solution%converged = .false. - Polarisation_solution%iterationsNeeded = totalIter + if (reason == -4) call IO_error(893_pInt) ! MPI error end function Polarisation_solution @@ -330,8 +327,8 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) polarAlpha, & polarBeta use mesh, only: & - grid3, & - grid + grid, & + grid3 use IO, only: & IO_intOut use math, only: & @@ -340,6 +337,10 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) math_mul3333xx33, & math_invSym3333, & math_mul33x33 + use debug, only: & + debug_level, & + debug_spectral, & + debug_spectralRotation use spectral_utilities, only: & wgt, & tensorField_real, & @@ -349,27 +350,17 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) Utilities_constitutiveResponse, & Utilities_divergenceRMS, & Utilities_curlRMS - use debug, only: & - debug_level, & - debug_spectral, & - debug_spectralRotation use homogenization, only: & materialpoint_dPdF use FEsolving, only: & terminallyIll implicit none -!-------------------------------------------------------------------------------------------------- -! strange syntax in the next line because otherwise macros expand beyond 132 character limit - 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 + 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, pointer, dimension(:,:,:,:,:) :: & F, & F_tau, & @@ -387,33 +378,33 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) SNESGetNumberFunctionEvals, & SNESGetIterationNumber - F => x_scal(1:3,1:3,1,& - XG_RANGE,YG_RANGE,ZG_RANGE) - F_tau => x_scal(1:3,1:3,2,& - XG_RANGE,YG_RANGE,ZG_RANGE) - residual_F => f_scal(1:3,1:3,1,& - X_RANGE,Y_RANGE,Z_RANGE) - residual_F_tau => f_scal(1:3,1:3,2,& - X_RANGE,Y_RANGE,Z_RANGE) - - call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) - call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) + F => x_scal(1:3,1:3,1,& + XG_RANGE,YG_RANGE,ZG_RANGE) + F_tau => x_scal(1:3,1:3,2,& + XG_RANGE,YG_RANGE,ZG_RANGE) + residual_F => f_scal(1:3,1:3,1,& + X_RANGE, Y_RANGE, Z_RANGE) + residual_F_tau => f_scal(1:3,1:3,2,& + X_RANGE, Y_RANGE, Z_RANGE) F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - if(nfuncs== 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment - newIteration: if(totalIter <= PETScIter) then + 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 !-------------------------------------------------------------------------------------------------- -! 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 @@ -440,18 +431,16 @@ 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. !-------------------------------------------------------------------------------------------------- ! calculate divergence tensorField_real = 0.0_pReal - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F + tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F !< stress field in disguise call utilities_FFTtensorForward() - err_div = Utilities_divergenceRMS() - call utilities_FFTtensorBackward() + err_div = Utilities_divergenceRMS() !< root mean squared error in divergence of stress !-------------------------------------------------------------------------------------------------- ! constructing residual @@ -471,8 +460,11 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F call utilities_FFTtensorForward() err_curl = Utilities_curlRMS() - call utilities_FFTtensorBackward() - + + nullify(F) + nullify(F_tau) + nullify(residual_F) + nullify(residual_F_tau) end subroutine Polarisation_formResidual @@ -487,8 +479,8 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason, err_div_tolAbs, & err_curl_tolRel, & err_curl_tolAbs, & - err_stress_tolAbs, & - err_stress_tolRel + err_stress_tolRel, & + err_stress_tolAbs use math, only: & math_mul3333xx33 use FEsolving, only: & @@ -507,24 +499,24 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason, real(pReal) :: & curlTol, & divTol, & - BC_tol + BCTol !-------------------------------------------------------------------------------------------------- ! stress BC handling F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc - err_BC = maxval(abs((-mask_stress+1.0_pReal)*math_mul3333xx33(C_scale,F_aim-F_av) + & - mask_stress *(P_av - params%stress_BC))) ! mask = 0.0 for no bc + err_BC = maxval(abs((1.0_pReal-mask_stress) * math_mul3333xx33(C_scale,F_aim-F_av) + & + mask_stress * (P_av-params%stress_BC))) ! mask = 0.0 for no bc !-------------------------------------------------------------------------------------------------- ! error calculation - curlTol = max(maxval(abs(F_aim-math_I3))*err_curl_tolRel,err_curl_tolAbs) - divTol = max(maxval(abs(P_av)) *err_div_tolRel,err_div_tolAbs) - BC_tol = max(maxval(abs(P_av)) *err_stress_tolrel,err_stress_tolabs) + curlTol = max(maxval(abs(F_aim-math_I3))*err_curl_tolRel ,err_curl_tolAbs) + divTol = max(maxval(abs(P_av)) *err_div_tolRel ,err_div_tolAbs) + BCTol = max(maxval(abs(P_av)) *err_stress_tolRel,err_stress_tolAbs) converged: if ((totalIter >= itmin .and. & - all([ err_div/divTol, & + all([ err_div /divTol, & err_curl/curlTol, & - err_BC/BC_tol ] < 1.0_pReal)) & + err_BC /BCTol ] < 1.0_pReal)) & .or. terminallyIll) then reason = 1 elseif (totalIter >= itmax) then converged @@ -536,12 +528,12 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason, !-------------------------------------------------------------------------------------------------- ! report write(6,'(1/,a)') ' ... reporting .............................................................' - write(6,'(/,a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', & - err_curl/curlTol,' (',err_curl,' -, tol =',curlTol,')' - write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & - err_div/divTol, ' (',err_div, ' / m, tol =',divTol,')' - write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', & - err_BC/BC_tol, ' (',err_BC, ' Pa, tol =',BC_tol,')' + write(6,'(/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & + err_div/divTol, ' (',err_div, ' / m, tol = ',divTol,')' + write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', & + err_curl/curlTol,' (',err_curl,' -, tol = ',curlTol,')' + write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', & + err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' write(6,'(/,a)') ' ===========================================================================' flush(6) @@ -549,6 +541,8 @@ end subroutine Polarisation_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 Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC) use math, only: & @@ -558,133 +552,140 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,deformati math_rotate_backward33 use numerics, only: & worldrank + use homogenization, only: & + materialpoint_F0 use mesh, only: & - grid3, & - grid - use spectral_utilities, only: & - Utilities_calculateRate, & - Utilities_forwardField, & - Utilities_updateIPcoords, & - tBoundaryCondition, & - cutBack - use IO, only: & - IO_write_JobRealFile - use FEsolving, only: & - restartWrite + 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, dimension(:,:,:,:), pointer :: xx_psc, F, F_tau - integer(pInt) :: i, j, k - real(pReal), dimension(3,3) :: F_lambda33 - character(len=1024) :: rankStr + 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, dimension(:,:,:,:), pointer :: FandF_tau, F, F_tau + integer(pInt) :: i, j, k + real(pReal), dimension(3,3) :: F_lambda33 + character(len=32) :: rankStr !-------------------------------------------------------------------------------------------------- ! update coordinates and rate and forward last inc - call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr) - F => xx_psc(0:8,:,:,:) - F_tau => xx_psc(9:17,:,:,:) - 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) - call IO_write_jobRealFile(777,'F_tau'//trim(rankStr),size(F_tau)) ! writing deformation gradient field to file - write (777,rec=1) F_tau - close (777) - call IO_write_jobRealFile(777,'F_tau_lastInc'//trim(rankStr),size(F_tau_lastInc)) ! writing F_lastInc field to file - write (777,rec=1) F_tau_lastInc - close (777) - if (worldrank == 0_pInt) then - call IO_write_jobRealFile(777,'F_aim',size(F_aim)) - write (777,rec=1) F_aim - close(777) - call IO_write_jobRealFile(777,'F_aim_lastInc',size(F_aim_lastInc)) - write (777,rec=1) F_aim_lastInc - close (777) - 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,FandF_tau,ierr); CHKERRQ(ierr) + F => FandF_tau( 0: 8,:,:,:) + F_tau => FandF_tau( 9:17,:,:,:) - call utilities_updateIPcoords(F) + 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) - if (cutBack) then - F_aim = F_aim_lastInc - F_tau= reshape(F_tau_lastInc,[9,grid(1),grid(2),grid3]) - 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 + 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,'F_aimDot',size(F_aimDot)) + write (777,rec=1) F_aimDot; close(777) + endif + + 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) + call IO_write_jobRealFile(777,'F_tau'//trim(rankStr),size(F_tau)) ! writing deformation gradient field to file + write (777,rec=1) F_tau; close (777) + call IO_write_jobRealFile(777,'F_tau_lastInc'//trim(rankStr),size(F_tau_lastInc)) ! writing F_tau_lastInc field to file + write (777,rec=1) F_tau_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_tauDot = Utilities_calculateRate(guess, & + F_tau_lastInc,reshape(F_tau,[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 + F_tau_lastInc = reshape(F_tau, [3,3,grid(1),grid(2),grid3]) ! winding F_tau 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_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]) - F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3]) - endif +! update average and local deformation gradients + F_aim = F_aim_lastInc + F_aimDot * timeinc - 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]) - F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), & ! does not have any average value as boundary condition - [9,grid(1),grid(2),grid3]) - if (.not. guess) then ! large strain forwarding + 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]) + if (guess) then + F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), & + [9,grid(1),grid(2),grid3]) ! does not have any average value as boundary condition + else do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) F_lambda33 = reshape(F_tau(1:9,i,j,k)-F(1:9,i,j,k),[3,3]) F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, & math_mul3333xx33(C_scale,& math_mul33x33(math_transpose33(F_lambda33),& - F_lambda33) -math_I3))*0.5_pReal)& + F_lambda33)-math_I3))*0.5_pReal)& + math_I3 F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(1:9,i,j,k) enddo; enddo; enddo endif - call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) + + nullify(F) + nullify(F_tau) + call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) end subroutine Polarisation_forward diff --git a/src/spectral_utilities.f90 b/src/spectral_utilities.f90 old mode 100644 new mode 100755 index 5a3361fcf..7b6f08a7d --- 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 @@ -145,8 +145,7 @@ module spectral_utilities FIELD_UNDEFINED_ID, & FIELD_MECH_ID, & FIELD_THERMAL_ID, & - FIELD_DAMAGE_ID, & - utilities_calcPlasticity + FIELD_DAMAGE_ID private :: & utilities_getFreqDerivative @@ -414,8 +413,7 @@ subroutine utilities_updateGamma(C,saveReference) write(6,'(/,a)') ' writing reference stiffness to file' flush(6) call IO_write_jobRealFile(777,'C_ref',size(C_ref)) - write (777,rec=1) C_ref - close(777) + write (777,rec=1) C_ref; close(777) endif endif @@ -800,7 +798,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 @@ -818,28 +816,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 @@ -925,10 +925,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: & @@ -941,31 +941,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 @@ -976,17 +967,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 @@ -1003,7 +986,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 @@ -1021,157 +1016,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 yield stress, plastic strain, total strain and their equivalent values -!-------------------------------------------------------------------------------------------------- -subroutine utilities_calcPlasticity(yieldStress, plasticStrain, eqStress, eqTotalStrain, & - eqPlasticStrain, plasticWork, rotation_BC) - use crystallite, only: & - crystallite_Fe, & - crystallite_P, & - crystallite_subF - use material, only: & - homogenization_maxNgrains - use mesh, only: & - mesh_maxNips,& - mesh_NcpElems - use math, only: & - math_det33, & - math_inv33, & - math_mul33x33, & - math_trace33, & - math_transpose33, & - math_equivStrain33, & - math_equivStress33, & - math_rotate_forward33, & - math_identity2nd, & - math_crossproduct, & - math_eigenvectorBasisSym, & - math_eigenvectorBasisSym33, & - math_eigenvectorBasisSym33_log, & - math_eigenValuesVectorsSym33 - - implicit none - - real(pReal), intent(inout) :: eqStress, eqPlasticStrain, plasticWork - real(pReal), intent(out) :: eqTotalStrain - real(pReal), dimension(3,3),intent(out) :: yieldStress, plasticStrain - real(pReal), intent(in), dimension(3,3) :: rotation_BC !< rotation of load frame - real(pReal), dimension(3,3) :: cauchy, P_av, F_av, Ve_av !< average - real(pReal), dimension(3) :: Values, S - real(pReal), dimension(3,3) :: Vectors, diag - real(pReal), dimension(3,3) :: & - Vp, F_temp, U, VT, R, V, V_total - real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & - Be, Ve, Fe - real(pReal), dimension(15) :: WORK !< previous deformation gradient - integer(pInt) :: INFO, i, j, k, l, ierr - real(pReal) :: wgtm - real(pReal) :: eqStressOld, eqPlasticStrainOld, plasticWorkOld - - external :: dgesvd - - eqStressOld = eqStress - eqPlasticStrainOld = eqPlasticStrain - plasticWorkOld = plasticWork - wgtm = 1.0_pReal/real(mesh_NcpElems*mesh_maxNips*homogenization_maxNgrains,pReal) - diag = 0.0_pReal - - P_av = sum(sum(sum(crystallite_P,dim=5),dim=4),dim=3) * wgtm - call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - P_av = math_rotate_forward33(P_av,rotation_BC) - - F_av = sum(sum(sum(crystallite_subF,dim=5),dim=4),dim=3) * wgtm - call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - F_av = math_rotate_forward33(F_av,rotation_BC) - - cauchy = 1.0_pReal/math_det33(F_av)*math_mul33x33(P_av,transpose(F_av)) - yieldStress = cauchy - eqStress = math_equivStress33(cauchy) - - F_temp = F_av - call dgesvd ('A', 'A', 3, 3, F_temp, 3, S, U, 3, VT, 3, WORK, 15, INFO) ! singular value decomposition - - R = math_mul33x33(U, VT) ! rotation of polar decomposition - V = math_mul33x33(F_av,math_inv33(R)) - - call math_eigenValuesVectorsSym33(V,Values,Vectors) - do l = 1_pInt, 3_pInt - if (Values(l) < 0.0_pReal) then - Values(l) = -Values(l) - Vectors(1:3, l) = -Vectors(1:3, l) - endif - Values(l) = log(Values(l)) - diag(l,l) = Values(l) - enddo - if (dot_product(Vectors(1:3,1),Vectors(1:3,2)) /= 0) then - Vectors(1:3,2) = math_crossproduct(Vectors(1:3,3), Vectors(1:3,1)) - Vectors(1:3,2) = Vectors(1:3,2)/sqrt(dot_product(Vectors(1:3,2),Vectors(1:3,2))) - endif - if (dot_product(Vectors(1:3,2),Vectors(1:3,3)) /= 0) then - Vectors(1:3,3) = math_crossproduct(Vectors(1:3,1), Vectors(1:3,2)) - Vectors(1:3,3) = Vectors(1:3,3)/sqrt(dot_product(Vectors(1:3,3),Vectors(1:3,3))) - endif - if (dot_product(Vectors(1:3,3),Vectors(1:3,1)) /= 0) then - Vectors(1:3,1) = math_crossproduct(Vectors(1:3,2), Vectors(1:3,3)) - Vectors(1:3,1) = Vectors(1:3,1)/sqrt(dot_product(Vectors(1:3,1),Vectors(1:3,1))) - endif - - V_total = REAL(math_mul33x33(Vectors, math_mul33x33(diag, transpose(Vectors)))) - eqTotalStrain = math_equivStrain33(V_total) - - do k = 1_pInt, mesh_NcpElems; do j = 1_pInt, mesh_maxNips; do i = 1_pInt,homogenization_maxNgrains - Fe(1:3,1:3,i,j,k) = crystallite_Fe(1:3,1:3,i,j,k) - Fe(1:3,1:3,i,j,k) = math_rotate_forward33(Fe(1:3,1:3,i,j,k),rotation_BC) - Be(1:3,1:3,i,j,k) = math_mul33x33(Fe(1:3,1:3,i,j,k),math_transpose33(Fe(1:3,1:3,i,j,k))) ! elastic part of left Cauchy–Green deformation tensor - Ve(1:3,1:3,i,j,k) = math_eigenvectorBasisSym33_log(Be(1:3,1:3,i,j,k)) - enddo; enddo; enddo - - Ve_av = sum(sum(sum(Ve,dim=5),dim=4),dim=3) * wgtm - call MPI_Allreduce(MPI_IN_PLACE,Ve_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - - Vp = V_total - Ve_av - - eqPlasticStrain = math_equivStrain33(Vp) - - plasticStrain = Vp - - plasticWork = plasticWorkOld + 0.5*(eqStressOld + eqStress) * (eqPlasticStrain - eqPlasticStrainOld) - -end subroutine utilities_calcPlasticity - !-------------------------------------------------------------------------------------------------- !> @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 @@ -1179,17 +1041,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