fixing spectral cutback hiccup and multiple cleanups
flush(6) at better places, added dedicated CPFEM_age subroutine, cleaned up cutback logic, fixed broken assignment of old timeinc, continueCalculation is now a logical, rearrnaged interfaces for utilities_constitutiveResponse and utilities_calculateRate, handling of stressBC more understandable, added more comments and explanations
This commit is contained in:
Philip Eisenlohr
parent
93073ed661
commit
b36151cc32
|
|
@ -162,6 +162,7 @@ subroutine CPFEM_init
|
|||
write(6,'(/,a)') ' <<<+- CPFEM init -+>>>'
|
||||
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
|
||||
#include "compilation_info.f90"
|
||||
flush(6)
|
||||
endif mainProcess
|
||||
|
||||
! initialize stress and jacobian to zero
|
||||
|
|
@ -242,8 +243,8 @@ subroutine CPFEM_init
|
|||
write(6,'(a32,1x,6(i8,1x))') 'CPFEM_dcsdE: ', shape(CPFEM_dcsdE)
|
||||
write(6,'(a32,1x,6(i8,1x),/)') 'CPFEM_dcsdE_knownGood: ', shape(CPFEM_dcsdE_knownGood)
|
||||
write(6,'(a32,l1)') 'symmetricSolver: ', symmetricSolver
|
||||
flush(6)
|
||||
endif
|
||||
flush(6)
|
||||
|
||||
end subroutine CPFEM_init
|
||||
|
||||
|
|
|
|||
158
src/CPFEM2.f90
158
src/CPFEM2.f90
|
|
@ -9,7 +9,7 @@ module CPFEM2
|
|||
private
|
||||
|
||||
public :: &
|
||||
CPFEM_general, &
|
||||
CPFEM_age, &
|
||||
CPFEM_initAll
|
||||
|
||||
contains
|
||||
|
|
@ -127,6 +127,7 @@ subroutine CPFEM_init
|
|||
write(6,'(/,a)') ' <<<+- CPFEM init -+>>>'
|
||||
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
|
||||
#include "compilation_info.f90"
|
||||
flush(6)
|
||||
endif mainProcess
|
||||
|
||||
! *** restore the last converged values of each essential variable from the binary file
|
||||
|
|
@ -194,7 +195,6 @@ subroutine CPFEM_init
|
|||
|
||||
restartRead = .false.
|
||||
endif
|
||||
flush(6)
|
||||
|
||||
end subroutine CPFEM_init
|
||||
|
||||
|
|
@ -202,7 +202,7 @@ end subroutine CPFEM_init
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief perform initialization at first call, update variables and call the actual material model
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine CPFEM_general(age, dt)
|
||||
subroutine CPFEM_age()
|
||||
use prec, only: &
|
||||
pReal, &
|
||||
pInt
|
||||
|
|
@ -215,7 +215,6 @@ subroutine CPFEM_general(age, dt)
|
|||
debug_levelExtensive, &
|
||||
debug_levelSelective
|
||||
use FEsolving, only: &
|
||||
terminallyIll, &
|
||||
restartWrite
|
||||
use math, only: &
|
||||
math_identity2nd, &
|
||||
|
|
@ -254,114 +253,99 @@ subroutine CPFEM_general(age, dt)
|
|||
crystallite_dPdF, &
|
||||
crystallite_Tstar0_v, &
|
||||
crystallite_Tstar_v
|
||||
use homogenization, only: &
|
||||
materialpoint_stressAndItsTangent, &
|
||||
materialpoint_postResults
|
||||
use IO, only: &
|
||||
IO_write_jobRealFile, &
|
||||
IO_warning
|
||||
use DAMASK_interface
|
||||
|
||||
implicit none
|
||||
real(pReal), intent(in) :: dt !< time increment
|
||||
logical, intent(in) :: age !< age results
|
||||
|
||||
integer(pInt) :: i, k, l, m, ph, homog, mySource
|
||||
character(len=1024) :: rankStr
|
||||
character(len=32) :: rankStr
|
||||
|
||||
!*** age results and write restart data if requested
|
||||
if (age) then
|
||||
crystallite_F0 = crystallite_partionedF ! crystallite deformation (_subF is perturbed...)
|
||||
crystallite_Fp0 = crystallite_Fp ! crystallite plastic deformation
|
||||
crystallite_Lp0 = crystallite_Lp ! crystallite plastic velocity
|
||||
crystallite_Fi0 = crystallite_Fi ! crystallite intermediate deformation
|
||||
crystallite_Li0 = crystallite_Li ! crystallite intermediate velocity
|
||||
crystallite_dPdF0 = crystallite_dPdF ! crystallite stiffness
|
||||
crystallite_Tstar0_v = crystallite_Tstar_v ! crystallite 2nd Piola Kirchhoff stress
|
||||
if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) &
|
||||
write(6,'(a)') '<< CPFEM >> aging states'
|
||||
|
||||
forall ( i = 1:size(plasticState )) plasticState(i)%state0 = plasticState(i)%state ! copy state in this lenghty way because: A component cannot be an array if the encompassing structure is an array
|
||||
do i = 1, size(sourceState)
|
||||
do mySource = 1,phase_Nsources(i)
|
||||
sourceState(i)%p(mySource)%state0 = sourceState(i)%p(mySource)%state ! copy state in this lenghty way because: A component cannot be an array if the encompassing structure is an array
|
||||
enddo; enddo
|
||||
if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) &
|
||||
write(6,'(a)') '<< CPFEM >> aging states'
|
||||
crystallite_F0 = crystallite_partionedF ! crystallite deformation (_subF is perturbed...)
|
||||
crystallite_Fp0 = crystallite_Fp ! crystallite plastic deformation
|
||||
crystallite_Lp0 = crystallite_Lp ! crystallite plastic velocity
|
||||
crystallite_Fi0 = crystallite_Fi ! crystallite intermediate deformation
|
||||
crystallite_Li0 = crystallite_Li ! crystallite intermediate velocity
|
||||
crystallite_dPdF0 = crystallite_dPdF ! crystallite stiffness
|
||||
crystallite_Tstar0_v = crystallite_Tstar_v ! crystallite 2nd Piola Kirchhoff stress
|
||||
|
||||
do homog = 1_pInt, material_Nhomogenization
|
||||
homogState (homog)%state0 = homogState (homog)%state
|
||||
thermalState (homog)%state0 = thermalState (homog)%state
|
||||
damageState (homog)%state0 = damageState (homog)%state
|
||||
vacancyfluxState (homog)%state0 = vacancyfluxState (homog)%state
|
||||
hydrogenfluxState(homog)%state0 = hydrogenfluxState(homog)%state
|
||||
enddo
|
||||
forall (i = 1:size(plasticState)) plasticState(i)%state0 = plasticState(i)%state ! copy state in this lengthy way because: A component cannot be an array if the encompassing structure is an array
|
||||
|
||||
do i = 1, size(sourceState)
|
||||
do mySource = 1,phase_Nsources(i)
|
||||
sourceState(i)%p(mySource)%state0 = sourceState(i)%p(mySource)%state ! copy state in this lengthy way because: A component cannot be an array if the encompassing structure is an array
|
||||
enddo; enddo
|
||||
|
||||
if (restartWrite) then
|
||||
if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) &
|
||||
write(6,'(a)') '<< CPFEM >> writing state variables of last converged step to binary files'
|
||||
|
||||
write(rankStr,'(a1,i0)')'_',worldrank
|
||||
do homog = 1_pInt, material_Nhomogenization
|
||||
homogState (homog)%state0 = homogState (homog)%state
|
||||
thermalState (homog)%state0 = thermalState (homog)%state
|
||||
damageState (homog)%state0 = damageState (homog)%state
|
||||
vacancyfluxState (homog)%state0 = vacancyfluxState (homog)%state
|
||||
hydrogenfluxState(homog)%state0 = hydrogenfluxState(homog)%state
|
||||
enddo
|
||||
|
||||
call IO_write_jobRealFile(777,'recordedPhase'//trim(rankStr),size(material_phase))
|
||||
write (777,rec=1) material_phase
|
||||
close (777)
|
||||
if (restartWrite) then
|
||||
if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) &
|
||||
write(6,'(a)') '<< CPFEM >> writing state variables of last converged step to binary files'
|
||||
|
||||
call IO_write_jobRealFile(777,'convergedF'//trim(rankStr),size(crystallite_F0))
|
||||
write (777,rec=1) crystallite_F0
|
||||
close (777)
|
||||
write(rankStr,'(a1,i0)')'_',worldrank
|
||||
|
||||
call IO_write_jobRealFile(777,'convergedFp'//trim(rankStr),size(crystallite_Fp0))
|
||||
write (777,rec=1) crystallite_Fp0
|
||||
close (777)
|
||||
call IO_write_jobRealFile(777,'recordedPhase'//trim(rankStr),size(material_phase))
|
||||
write (777,rec=1) material_phase; close (777)
|
||||
|
||||
call IO_write_jobRealFile(777,'convergedFi'//trim(rankStr),size(crystallite_Fi0))
|
||||
write (777,rec=1) crystallite_Fi0
|
||||
close (777)
|
||||
call IO_write_jobRealFile(777,'convergedF'//trim(rankStr),size(crystallite_F0))
|
||||
write (777,rec=1) crystallite_F0; close (777)
|
||||
|
||||
call IO_write_jobRealFile(777,'convergedLp'//trim(rankStr),size(crystallite_Lp0))
|
||||
write (777,rec=1) crystallite_Lp0
|
||||
close (777)
|
||||
call IO_write_jobRealFile(777,'convergedFp'//trim(rankStr),size(crystallite_Fp0))
|
||||
write (777,rec=1) crystallite_Fp0; close (777)
|
||||
|
||||
call IO_write_jobRealFile(777,'convergedLi'//trim(rankStr),size(crystallite_Li0))
|
||||
write (777,rec=1) crystallite_Li0
|
||||
close (777)
|
||||
call IO_write_jobRealFile(777,'convergedFi'//trim(rankStr),size(crystallite_Fi0))
|
||||
write (777,rec=1) crystallite_Fi0; close (777)
|
||||
|
||||
call IO_write_jobRealFile(777,'convergeddPdF'//trim(rankStr),size(crystallite_dPdF0))
|
||||
write (777,rec=1) crystallite_dPdF0
|
||||
close (777)
|
||||
call IO_write_jobRealFile(777,'convergedLp'//trim(rankStr),size(crystallite_Lp0))
|
||||
write (777,rec=1) crystallite_Lp0; close (777)
|
||||
|
||||
call IO_write_jobRealFile(777,'convergedTstar'//trim(rankStr),size(crystallite_Tstar0_v))
|
||||
write (777,rec=1) crystallite_Tstar0_v
|
||||
close (777)
|
||||
call IO_write_jobRealFile(777,'convergedLi'//trim(rankStr),size(crystallite_Li0))
|
||||
write (777,rec=1) crystallite_Li0; close (777)
|
||||
|
||||
call IO_write_jobRealFile(777,'convergedStateConst'//trim(rankStr))
|
||||
m = 0_pInt
|
||||
writePlasticityInstances: do ph = 1_pInt, size(phase_plasticity)
|
||||
do k = 1_pInt, plasticState(ph)%sizeState
|
||||
do l = 1, size(plasticState(ph)%state0(1,:))
|
||||
m = m+1_pInt
|
||||
write(777,rec=m) plasticState(ph)%state0(k,l)
|
||||
enddo; enddo
|
||||
enddo writePlasticityInstances
|
||||
close (777)
|
||||
call IO_write_jobRealFile(777,'convergeddPdF'//trim(rankStr),size(crystallite_dPdF0))
|
||||
write (777,rec=1) crystallite_dPdF0; close (777)
|
||||
|
||||
call IO_write_jobRealFile(777,'convergedStateHomog'//trim(rankStr))
|
||||
m = 0_pInt
|
||||
writeHomogInstances: do homog = 1_pInt, material_Nhomogenization
|
||||
do k = 1_pInt, homogState(homog)%sizeState
|
||||
do l = 1, size(homogState(homog)%state0(1,:))
|
||||
m = m+1_pInt
|
||||
write(777,rec=m) homogState(homog)%state0(k,l)
|
||||
enddo; enddo
|
||||
enddo writeHomogInstances
|
||||
close (777)
|
||||
call IO_write_jobRealFile(777,'convergedTstar'//trim(rankStr),size(crystallite_Tstar0_v))
|
||||
write (777,rec=1) crystallite_Tstar0_v; close (777)
|
||||
|
||||
endif
|
||||
endif
|
||||
call IO_write_jobRealFile(777,'convergedStateConst'//trim(rankStr))
|
||||
m = 0_pInt
|
||||
writePlasticityInstances: do ph = 1_pInt, size(phase_plasticity)
|
||||
do k = 1_pInt, plasticState(ph)%sizeState
|
||||
do l = 1, size(plasticState(ph)%state0(1,:))
|
||||
m = m+1_pInt
|
||||
write(777,rec=m) plasticState(ph)%state0(k,l)
|
||||
enddo; enddo
|
||||
enddo writePlasticityInstances
|
||||
close (777)
|
||||
|
||||
if (.not. terminallyIll) &
|
||||
call materialpoint_stressAndItsTangent(.True., dt)
|
||||
call IO_write_jobRealFile(777,'convergedStateHomog'//trim(rankStr))
|
||||
m = 0_pInt
|
||||
writeHomogInstances: do homog = 1_pInt, material_Nhomogenization
|
||||
do k = 1_pInt, homogState(homog)%sizeState
|
||||
do l = 1, size(homogState(homog)%state0(1,:))
|
||||
m = m+1_pInt
|
||||
write(777,rec=m) homogState(homog)%state0(k,l)
|
||||
enddo; enddo
|
||||
enddo writeHomogInstances
|
||||
close (777)
|
||||
|
||||
end subroutine CPFEM_general
|
||||
endif
|
||||
|
||||
if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) &
|
||||
write(6,'(a)') '<< CPFEM >> done aging states'
|
||||
|
||||
end subroutine CPFEM_age
|
||||
|
||||
end module CPFEM2
|
||||
|
|
|
|||
|
|
@ -442,8 +442,9 @@ program DAMASK_spectral
|
|||
if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_seek')
|
||||
|
||||
if (.not. appendToOutFile) then ! if not restarting, write 0th increment
|
||||
write(6,'(1/,a)') ' ... writing initial configuration to file ........................'
|
||||
do i = 1, size(materialpoint_results,3)/(maxByteOut/(materialpoint_sizeResults*pReal))+1 ! slice the output of my process in chunks not exceeding the limit for one output
|
||||
outputIndex = int([(i-1_pInt)*((maxRealOut)/materialpoint_sizeResults)+1_pInt, &
|
||||
outputIndex = int([(i-1_pInt)*((maxRealOut)/materialpoint_sizeResults)+1_pInt, & ! QUESTION: why not starting i at 0 instead of murky 1?
|
||||
min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt)
|
||||
call MPI_file_write(resUnit, &
|
||||
reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)), &
|
||||
|
|
@ -453,24 +454,23 @@ program DAMASK_spectral
|
|||
if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_write')
|
||||
enddo
|
||||
fileOffset = fileOffset + sum(outputSize) ! forward to current file position
|
||||
write(6,'(1/,a)') ' ... writing initial configuration to file ........................'
|
||||
endif
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! loopping over loadcases
|
||||
! looping over loadcases
|
||||
loadCaseLooping: do currentLoadCase = 1_pInt, size(loadCases)
|
||||
time0 = time ! currentLoadCase start time
|
||||
guess = loadCases(currentLoadCase)%followFormerTrajectory ! change of load case? homogeneous guess for the first inc
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! loop oper incs defined in input file for current currentLoadCase
|
||||
! loop over incs defined in input file for current currentLoadCase
|
||||
incLooping: do inc = 1_pInt, loadCases(currentLoadCase)%incs
|
||||
totalIncsCounter = totalIncsCounter + 1_pInt
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! forwarding time
|
||||
timeIncOld = timeinc
|
||||
timeIncOld = timeinc ! last timeinc that brought former inc to an end
|
||||
if (loadCases(currentLoadCase)%logscale == 0_pInt) then ! linear scale
|
||||
timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal) ! only valid for given linear time scale. will be overwritten later in case loglinear scale is used
|
||||
timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal)
|
||||
else
|
||||
if (currentLoadCase == 1_pInt) then ! 1st currentLoadCase of logarithmic scale
|
||||
if (inc == 1_pInt) then ! 1st inc of 1st currentLoadCase of logarithmic scale
|
||||
|
|
@ -486,20 +486,23 @@ program DAMASK_spectral
|
|||
real(loadCases(currentLoadCase)%incs ,pReal)))
|
||||
endif
|
||||
endif
|
||||
timeinc = timeinc / 2.0_pReal**real(cutBackLevel,pReal) ! depending on cut back level, decrease time step
|
||||
timeinc = timeinc / real(subStepFactor,pReal)**real(cutBackLevel,pReal) ! depending on cut back level, decrease time step
|
||||
|
||||
forwarding: if (totalIncsCounter >= restartInc) then
|
||||
stepFraction = 0_pInt
|
||||
skipping: if (totalIncsCounter < restartInc) then ! not yet at restart inc?
|
||||
time = time + timeinc ! just advance time, skip already performed calculation
|
||||
guess = .true. ! QUESTION:why forced guessing instead of inheriting loadcase preference
|
||||
else skipping
|
||||
stepFraction = 0_pInt ! fraction scaled by stepFactor**cutLevel
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! loop over sub incs
|
||||
subIncLooping: do while (stepFraction/subStepFactor**cutBackLevel <1_pInt)
|
||||
time = time + timeinc ! forward time
|
||||
stepFraction = stepFraction + 1_pInt
|
||||
remainingLoadCaseTime = time0 - time + loadCases(currentLoadCase)%time + timeInc
|
||||
! loop over sub step
|
||||
subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel)
|
||||
remainingLoadCaseTime = loadCases(currentLoadCase)%time+time0 - time
|
||||
time = time + timeinc ! forward target time
|
||||
stepFraction = stepFraction + 1_pInt ! count step
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! report begin of new increment
|
||||
! report begin of new step
|
||||
write(6,'(/,a)') ' ###########################################################################'
|
||||
write(6,'(1x,a,es12.5'//&
|
||||
',a,'//IO_intOut(inc)//',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//&
|
||||
|
|
@ -509,11 +512,11 @@ program DAMASK_spectral
|
|||
's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,&
|
||||
'-', stepFraction, '/', subStepFactor**cutBackLevel,&
|
||||
' of load case ', currentLoadCase,'/',size(loadCases)
|
||||
flush(6)
|
||||
write(incInfo,'(a,'//IO_intOut(totalIncsCounter)//',a,'//IO_intOut(sum(loadCases%incs))//&
|
||||
',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//')') &
|
||||
'Increment ',totalIncsCounter,'/',sum(loadCases%incs),&
|
||||
'-',stepFraction, '/', subStepFactor**cutBackLevel
|
||||
flush(6)
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! forward fields
|
||||
|
|
@ -542,7 +545,7 @@ program DAMASK_spectral
|
|||
end select
|
||||
|
||||
case(FIELD_THERMAL_ID); call spectral_thermal_forward()
|
||||
case(FIELD_DAMAGE_ID); call spectral_damage_forward()
|
||||
case(FIELD_DAMAGE_ID); call spectral_damage_forward()
|
||||
end select
|
||||
enddo
|
||||
|
||||
|
|
@ -582,65 +585,64 @@ program DAMASK_spectral
|
|||
solres(field) = spectral_damage_solution(timeinc,timeIncOld,remainingLoadCaseTime)
|
||||
|
||||
end select
|
||||
|
||||
if (.not. solres(field)%converged) exit ! no solution found
|
||||
|
||||
enddo
|
||||
stagIter = stagIter + 1_pInt
|
||||
stagIterate = stagIter < stagItMax .and. &
|
||||
all(solres(:)%converged) .and. &
|
||||
.not. all(solres(:)%stagConverged)
|
||||
stagIterate = stagIter < stagItMax &
|
||||
.and. all(solres(:)%converged) &
|
||||
.and. .not. all(solres(:)%stagConverged) ! stationary with respect to staggered iteration
|
||||
enddo
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! check solution
|
||||
cutBack = .False.
|
||||
if(solres(1)%termIll .or. .not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found
|
||||
if (cutBackLevel < maxCutBack) then ! do cut back
|
||||
write(6,'(/,a)') ' cut back detected'
|
||||
cutBack = .True.
|
||||
stepFraction = (stepFraction - 1_pInt) * subStepFactor ! adjust to new denominator
|
||||
cutBackLevel = cutBackLevel + 1_pInt
|
||||
time = time - timeinc ! rewind time
|
||||
timeinc = timeinc/2.0_pReal
|
||||
elseif (solres(1)%termIll) then ! material point model cannot find a solution, exit in any casy
|
||||
call IO_warning(850_pInt)
|
||||
call MPI_file_close(resUnit,ierr)
|
||||
close(statUnit)
|
||||
call quit(-1_pInt*(lastRestartWritten+1_pInt)) ! quit and provide information about last restart inc written
|
||||
elseif (continueCalculation == 1_pInt) then
|
||||
guess = .true. ! accept non converged BVP solution
|
||||
else ! default behavior, exit if spectral solver does not converge
|
||||
call IO_warning(850_pInt)
|
||||
call MPI_file_close(resUnit,ierr)
|
||||
close(statUnit)
|
||||
call quit(-1_pInt*(lastRestartWritten+1_pInt)) ! quit and provide information about last restart inc written
|
||||
endif
|
||||
else
|
||||
! check solution for either advance or retry
|
||||
|
||||
if ( (continueCalculation .or. all(solres(:)%converged .and. solres(:)%stagConverged)) & ! don't care or did converge
|
||||
.and. .not. solres(1)%termIll) then ! and acceptable solution found
|
||||
timeIncOld = timeinc
|
||||
cutBack = .false.
|
||||
guess = .true. ! start guessing after first converged (sub)inc
|
||||
endif
|
||||
if (.not. cutBack) then
|
||||
if (worldrank == 0) then
|
||||
write(statUnit,*) totalIncsCounter, time, cutBackLevel, &
|
||||
solres%converged, solres%iterationsNeeded ! write statistics about accepted solution
|
||||
solres%converged, solres%iterationsNeeded
|
||||
flush(statUnit)
|
||||
endif
|
||||
elseif (cutBackLevel < maxCutBack) then ! further cutbacking tolerated?
|
||||
cutBack = .true.
|
||||
stepFraction = (stepFraction - 1_pInt) * subStepFactor ! adjust to new denominator
|
||||
cutBackLevel = cutBackLevel + 1_pInt
|
||||
time = time - timeinc ! rewind time
|
||||
timeinc = timeinc/real(subStepFactor,pReal) ! cut timestep
|
||||
write(6,'(/,a)') ' cutting back '
|
||||
else ! no more options to continue
|
||||
call IO_warning(850_pInt)
|
||||
call MPI_file_close(resUnit,ierr)
|
||||
close(statUnit)
|
||||
call quit(-1_pInt*(lastRestartWritten+1_pInt)) ! quit and provide information about last restart inc written
|
||||
endif
|
||||
enddo subIncLooping
|
||||
|
||||
enddo subStepLooping
|
||||
|
||||
cutBackLevel = max(0_pInt, cutBackLevel - 1_pInt) ! try half number of subincs next inc
|
||||
if(all(solres(:)%converged)) then ! report converged inc
|
||||
|
||||
if (all(solres(:)%converged)) then
|
||||
convergedCounter = convergedCounter + 1_pInt
|
||||
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') &
|
||||
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report converged inc
|
||||
' increment ', totalIncsCounter, ' converged'
|
||||
else
|
||||
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report non-converged inc
|
||||
' increment ', totalIncsCounter, ' NOT converged'
|
||||
notConvergedCounter = notConvergedCounter + 1_pInt
|
||||
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report non-converged inc
|
||||
' increment ', totalIncsCounter, ' NOT converged'
|
||||
endif; flush(6)
|
||||
|
||||
if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency
|
||||
if (worldrank == 0) &
|
||||
write(6,'(1/,a)') ' ... writing results to file ......................................'
|
||||
flush(6)
|
||||
call materialpoint_postResults()
|
||||
call MPI_file_seek (resUnit,fileOffset,MPI_SEEK_SET,ierr)
|
||||
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_file_seek')
|
||||
if (ierr /= 0_pInt) call IO_error(894_pInt, ext_msg='MPI_file_seek')
|
||||
do i=1, size(materialpoint_results,3)/(maxByteOut/(materialpoint_sizeResults*pReal))+1 ! slice the output of my process in chunks not exceeding the limit for one output
|
||||
outputIndex=int([(i-1_pInt)*((maxRealOut)/materialpoint_sizeResults)+1_pInt, &
|
||||
min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt)
|
||||
|
|
@ -652,15 +654,12 @@ program DAMASK_spectral
|
|||
enddo
|
||||
fileOffset = fileOffset + sum(outputSize) ! forward to current file position
|
||||
endif
|
||||
if( loadCases(currentLoadCase)%restartFrequency > 0_pInt .and. & ! at frequency of writing restart information set restart parameter for FEsolving
|
||||
mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0_pInt) then ! first call to CPFEM_general will write?
|
||||
restartWrite = .true.
|
||||
lastRestartWritten = inc
|
||||
if ( loadCases(currentLoadCase)%restartFrequency > 0_pInt & ! writing of restart info requested ...
|
||||
.and. mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0_pInt) then ! ... and at frequency of writing restart information
|
||||
restartWrite = .true. ! set restart parameter for FEsolving
|
||||
lastRestartWritten = inc ! QUESTION: first call to CPFEM_general will write?
|
||||
endif
|
||||
else forwarding
|
||||
time = time + timeinc
|
||||
guess = .true.
|
||||
endif forwarding
|
||||
endif skipping
|
||||
|
||||
enddo incLooping
|
||||
enddo loadCaseLooping
|
||||
|
|
@ -673,6 +672,7 @@ program DAMASK_spectral
|
|||
real(convergedCounter, pReal)/&
|
||||
real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, &
|
||||
' %) increments converged!'
|
||||
flush(6)
|
||||
call MPI_file_close(resUnit,ierr)
|
||||
close(statUnit)
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -213,8 +213,9 @@ subroutine AL_init
|
|||
endif restart
|
||||
|
||||
call Utilities_updateIPcoords(reshape(F,shape(F_lastInc)))
|
||||
call Utilities_constitutiveResponse(F_lastInc, reshape(F,shape(F_lastInc)), &
|
||||
0.0_pReal,P,C_volAvg,C_minMaxAvg,temp33_Real,.false.,math_I3)
|
||||
call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, &
|
||||
reshape(F,shape(F_lastInc)), 0.0_pReal, math_I3)
|
||||
|
||||
nullify(F)
|
||||
nullify(F_lambda)
|
||||
call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! write data back to PETSc
|
||||
|
|
@ -364,12 +365,10 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
|
|||
DMDALocalInfo, dimension(&
|
||||
DMDA_LOCAL_INFO_SIZE) :: &
|
||||
in
|
||||
PetscScalar, target, dimension(3,3,2, &
|
||||
XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: &
|
||||
x_scal
|
||||
PetscScalar, target, dimension(3,3,2, &
|
||||
X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: &
|
||||
f_scal
|
||||
PetscScalar, &
|
||||
target, dimension(3,3,2, XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: x_scal
|
||||
PetscScalar, &
|
||||
target, dimension(3,3,2, X_RANGE, Y_RANGE, Z_RANGE), intent(out) :: f_scal
|
||||
PetscScalar, pointer, dimension(:,:,:,:,:) :: &
|
||||
F, &
|
||||
F_lambda, &
|
||||
|
|
@ -441,8 +440,9 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
! evaluate constitutive response
|
||||
P_avLastEval = P_av
|
||||
call Utilities_constitutiveResponse(F_lastInc,F - residual_F_lambda/polarBeta,params%timeinc, &
|
||||
residual_F,C_volAvg,C_minMaxAvg,P_av,ForwardData,params%rotation_BC)
|
||||
|
||||
call Utilities_constitutiveResponse(residual_F,P_av,C_volAvg,C_minMaxAvg, &
|
||||
F - residual_F_lambda/polarBeta,params%timeinc, params%rotation_BC)
|
||||
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
|
||||
ForwardData = .False.
|
||||
|
||||
|
|
@ -655,10 +655,12 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stre
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
! update coordinates and rate and forward last inc
|
||||
call utilities_updateIPcoords(F)
|
||||
Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
|
||||
timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]))
|
||||
F_lambdaDot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
|
||||
timeinc_old,guess,F_lambda_lastInc,reshape(F_lambda,[3,3,grid(1),grid(2),grid3]))
|
||||
Fdot = Utilities_calculateRate(guess, &
|
||||
F_lastInc, reshape(F, [3,3,grid(1),grid(2),grid3]), timeinc_old, &
|
||||
math_rotate_backward33(f_aimDot,rotation_BC))
|
||||
F_lambdaDot = Utilities_calculateRate(guess, &
|
||||
F_lambda_lastInc,reshape(F_lambda,[3,3,grid(1),grid(2),grid3]), timeinc_old, &
|
||||
math_rotate_backward33(f_aimDot,rotation_BC))
|
||||
F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3])
|
||||
F_lambda_lastInc = reshape(F_lambda,[3,3,grid(1),grid(2),grid3])
|
||||
endif
|
||||
|
|
|
|||
|
|
@ -39,16 +39,16 @@ module spectral_mech_basic
|
|||
! stress, stiffness and compliance average etc.
|
||||
real(pReal), private, dimension(3,3) :: &
|
||||
F_aim = math_I3, &
|
||||
F_aim_lastIter = math_I3, &
|
||||
F_aim_lastInc = math_I3, &
|
||||
P_av = 0.0_pReal, &
|
||||
F_aimDot=0.0_pReal
|
||||
F_aimDot = 0.0_pReal
|
||||
character(len=1024), private :: incInfo
|
||||
real(pReal), private, dimension(3,3,3,3) :: &
|
||||
C_volAvg = 0.0_pReal, & !< current volume average stiffness
|
||||
C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness
|
||||
C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 stiffness
|
||||
S = 0.0_pReal !< current compliance (filled up with zeros)
|
||||
C_minMaxAvgLastInc = 0.0_pReal, & !< previous (min+max)/2 stiffness
|
||||
S = 0.0_pReal !< current compliance (filled up with zeros)
|
||||
real(pReal), private :: err_stress, err_div
|
||||
logical, private :: ForwardData
|
||||
integer(pInt), private :: &
|
||||
|
|
@ -69,7 +69,7 @@ module spectral_mech_basic
|
|||
contains
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief allocates all neccessary fields and fills them with data, potentially from restart info
|
||||
!> @brief allocates all necessary fields and fills them with data, potentially from restart info
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine basicPETSc_init
|
||||
#ifdef __GFORTRAN__
|
||||
|
|
@ -90,6 +90,8 @@ subroutine basicPETSc_init
|
|||
use numerics, only: &
|
||||
worldrank, &
|
||||
worldsize
|
||||
use homogenization, only: &
|
||||
materialpoint_F0
|
||||
use DAMASK_interface, only: &
|
||||
getSolverJobName
|
||||
use spectral_utilities, only: &
|
||||
|
|
@ -172,14 +174,11 @@ subroutine basicPETSc_init
|
|||
flush(6)
|
||||
write(rankStr,'(a1,i0)')'_',worldrank
|
||||
call IO_read_realFile(777,'F'//trim(rankStr),trim(getSolverJobName()),size(F))
|
||||
read (777,rec=1) F
|
||||
close (777)
|
||||
read (777,rec=1) F; close (777)
|
||||
call IO_read_realFile(777,'F_lastInc'//trim(rankStr),trim(getSolverJobName()),size(F_lastInc))
|
||||
read (777,rec=1) F_lastInc
|
||||
close (777)
|
||||
read (777,rec=1) F_lastInc; close (777)
|
||||
call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot))
|
||||
read (777,rec=1) f_aimDot
|
||||
close (777)
|
||||
read (777,rec=1) f_aimDot; close (777)
|
||||
F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F
|
||||
F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc
|
||||
elseif (restartInc == 1_pInt) then restart
|
||||
|
|
@ -187,41 +186,36 @@ subroutine basicPETSc_init
|
|||
F = reshape(F_lastInc,[9,grid(1),grid(2),grid3])
|
||||
endif restart
|
||||
|
||||
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
|
||||
call Utilities_updateIPcoords(reshape(F,shape(F_lastInc)))
|
||||
call Utilities_constitutiveResponse(F_lastInc, reshape(F,shape(F_lastInc)), &
|
||||
0.0_pReal, &
|
||||
P, &
|
||||
C_volAvg,C_minMaxAvg, & ! global average of stiffness and (min+max)/2
|
||||
temp33_Real, &
|
||||
.false., &
|
||||
math_I3)
|
||||
call Utilities_constitutiveResponse(P, temp33_Real, C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
|
||||
reshape(F,shape(F_lastInc)), & ! target F
|
||||
0.0_pReal, & ! time increment
|
||||
math_I3) ! no rotation of boundary condition
|
||||
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! write data back to PETSc
|
||||
! QUESTION: why not writing back right after reading (l.189)?
|
||||
|
||||
restartRead: if (restartInc > 1_pInt) then
|
||||
restartRead: if (restartInc > 1_pInt) then ! QUESTION: are those values not calc'ed by constitutiveResponse? why reading from file?
|
||||
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) &
|
||||
write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
|
||||
'reading more values of increment', restartInc - 1_pInt, 'from file'
|
||||
'reading more values of increment', restartInc-1_pInt, 'from file'
|
||||
flush(6)
|
||||
call IO_read_realFile(777,'C_volAvg',trim(getSolverJobName()),size(C_volAvg))
|
||||
read (777,rec=1) C_volAvg
|
||||
close (777)
|
||||
read (777,rec=1) C_volAvg; close (777)
|
||||
call IO_read_realFile(777,'C_volAvgLastInc',trim(getSolverJobName()),size(C_volAvgLastInc))
|
||||
read (777,rec=1) C_volAvgLastInc
|
||||
close (777)
|
||||
read (777,rec=1) C_volAvgLastInc; close (777)
|
||||
call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg))
|
||||
read (777,rec=1) C_minMaxAvg
|
||||
close (777)
|
||||
read (777,rec=1) C_minMaxAvg; close (777)
|
||||
endif restartRead
|
||||
|
||||
call Utilities_updateGamma(C_minmaxAvg,.True.)
|
||||
call Utilities_updateGamma(C_minmaxAvg,.true.)
|
||||
|
||||
end subroutine basicPETSc_init
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief solution for the Basic PETSC scheme with internal iterations
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
type(tSolutionState) function &
|
||||
basicPETSc_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC)
|
||||
type(tSolutionState) function basicPETSc_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC)
|
||||
use IO, only: &
|
||||
IO_error
|
||||
use numerics, only: &
|
||||
|
|
@ -238,13 +232,13 @@ type(tSolutionState) function &
|
|||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! input data for solution
|
||||
real(pReal), intent(in) :: &
|
||||
timeinc, & !< increment in time for current solution
|
||||
timeinc_old !< increment in time of last increment
|
||||
type(tBoundaryCondition), intent(in) :: &
|
||||
stress_BC
|
||||
character(len=*), intent(in) :: &
|
||||
incInfoIn
|
||||
real(pReal), intent(in) :: &
|
||||
timeinc, & !< increment time for current solution
|
||||
timeinc_old !< increment time of last successful increment
|
||||
type(tBoundaryCondition), intent(in) :: &
|
||||
stress_BC
|
||||
real(pReal), dimension(3,3), intent(in) :: rotation_BC
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
@ -279,14 +273,13 @@ type(tSolutionState) function &
|
|||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! check convergence
|
||||
call SNESGetConvergedReason(snes,reason,ierr)
|
||||
CHKERRQ(ierr)
|
||||
call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr)
|
||||
|
||||
BasicPETSc_solution%converged = reason > 0
|
||||
basicPETSC_solution%iterationsNeeded = totalIter
|
||||
basicPETSc_solution%termIll = terminallyIll
|
||||
terminallyIll = .false.
|
||||
BasicPETSc_solution%converged =.true.
|
||||
if (reason == -4) call IO_error(893_pInt)
|
||||
if (reason < 1) basicPETSC_solution%converged = .false.
|
||||
basicPETSC_solution%iterationsNeeded = totalIter
|
||||
if (reason == -4) call IO_error(893_pInt) ! MPI error
|
||||
|
||||
end function BasicPETSc_solution
|
||||
|
||||
|
|
@ -322,19 +315,18 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
|
|||
terminallyIll
|
||||
|
||||
implicit none
|
||||
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &
|
||||
in
|
||||
PetscScalar, dimension(3,3, &
|
||||
XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: &
|
||||
x_scal
|
||||
PetscScalar, dimension(3,3, &
|
||||
X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: &
|
||||
f_scal
|
||||
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in
|
||||
PetscScalar, &
|
||||
dimension(3,3, XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: x_scal !< what is this?
|
||||
PetscScalar, &
|
||||
dimension(3,3, X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: f_scal !< what is this?
|
||||
PetscInt :: &
|
||||
PETScIter, &
|
||||
nfuncs
|
||||
PetscObject :: dummy
|
||||
PetscErrorCode :: ierr
|
||||
real(pReal), dimension(3,3) :: &
|
||||
deltaF_aim
|
||||
|
||||
external :: &
|
||||
SNESGetNumberFunctionEvals, &
|
||||
|
|
@ -343,46 +335,45 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
|
|||
call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr)
|
||||
call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr)
|
||||
|
||||
if(nfuncs== 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment
|
||||
newIteration: if(totalIter <= PETScIter) then
|
||||
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! report begin of new iteration
|
||||
! begin of new iteration
|
||||
newIteration: if (totalIter <= PETScIter) then
|
||||
totalIter = totalIter + 1_pInt
|
||||
write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), &
|
||||
' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
|
||||
write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') &
|
||||
trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
|
||||
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
|
||||
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', &
|
||||
math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC))
|
||||
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', &
|
||||
math_transpose33(F_aim)
|
||||
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
|
||||
' deformation gradient aim (lab) =', math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC))
|
||||
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
|
||||
' deformation gradient aim =', math_transpose33(F_aim)
|
||||
flush(6)
|
||||
endif newIteration
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! evaluate constitutive response
|
||||
call Utilities_constitutiveResponse(F_lastInc,x_scal,params%timeinc, &
|
||||
f_scal,C_volAvg,C_minmaxAvg,P_av,ForwardData,params%rotation_BC)
|
||||
call Utilities_constitutiveResponse(f_scal,P_av,C_volAvg,C_minmaxAvg, &
|
||||
x_scal,params%timeinc, params%rotation_BC)
|
||||
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
|
||||
ForwardData = .false.
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! stress BC handling
|
||||
F_aim_lastIter = F_aim
|
||||
F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc
|
||||
err_stress = maxval(abs(mask_stress * (P_av - params%stress_BC))) ! mask = 0.0 for no bc
|
||||
deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC)
|
||||
F_aim = F_aim - deltaF_aim
|
||||
err_stress = maxval(abs(mask_stress * (P_av - params%stress_BC))) ! mask = 0.0 when no stress bc
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! updated deformation gradient using fix point algorithm of basic scheme
|
||||
tensorField_real = 0.0_pReal
|
||||
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = f_scal
|
||||
call utilities_FFTtensorForward()
|
||||
err_div = Utilities_divergenceRMS()
|
||||
call utilities_fourierGammaConvolution(math_rotate_backward33(F_aim_lastIter-F_aim,params%rotation_BC))
|
||||
call utilities_FFTtensorBackward()
|
||||
call utilities_FFTtensorForward() ! FFT forward of global "tensorField_real"
|
||||
err_div = Utilities_divergenceRMS() ! divRMS of tensorField_fourier
|
||||
call utilities_fourierGammaConvolution(math_rotate_backward33(deltaF_aim,params%rotation_BC)) ! convolution of Gamma and tensorField_fourier, with arg
|
||||
call utilities_FFTtensorBackward() ! FFT backward of global tensorField_fourier
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! constructing residual
|
||||
f_scal = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3)
|
||||
f_scal = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too
|
||||
|
||||
end subroutine BasicPETSc_formResidual
|
||||
|
||||
|
|
@ -443,106 +434,120 @@ end subroutine BasicPETSc_converged
|
|||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief forwarding routine
|
||||
!> @details find new boundary conditions and best F estimate for end of current timestep
|
||||
!> possibly writing restart information, triggering of state increment in DAMASK, and updating of IPcoordinates
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC)
|
||||
use math, only: &
|
||||
math_mul33x33 ,&
|
||||
math_rotate_backward33
|
||||
use numerics, only: &
|
||||
worldrank
|
||||
use mesh, only: &
|
||||
grid, &
|
||||
grid3
|
||||
use spectral_utilities, only: &
|
||||
Utilities_calculateRate, &
|
||||
Utilities_forwardField, &
|
||||
Utilities_updateIPcoords, &
|
||||
tBoundaryCondition, &
|
||||
cutBack
|
||||
use IO, only: &
|
||||
IO_write_JobRealFile
|
||||
use FEsolving, only: &
|
||||
restartWrite
|
||||
use math, only: &
|
||||
math_mul33x33 ,&
|
||||
math_rotate_backward33
|
||||
use numerics, only: &
|
||||
worldrank
|
||||
use homogenization, only: &
|
||||
materialpoint_F0
|
||||
use mesh, only: &
|
||||
grid, &
|
||||
grid3
|
||||
use CPFEM2, only: &
|
||||
CPFEM_age
|
||||
use spectral_utilities, only: &
|
||||
Utilities_calculateRate, &
|
||||
Utilities_forwardField, &
|
||||
Utilities_updateIPcoords, &
|
||||
tBoundaryCondition, &
|
||||
cutBack
|
||||
use IO, only: &
|
||||
IO_write_JobRealFile
|
||||
use FEsolving, only: &
|
||||
restartWrite
|
||||
|
||||
implicit none
|
||||
real(pReal), intent(in) :: &
|
||||
timeinc_old, &
|
||||
timeinc, &
|
||||
loadCaseTime !< remaining time of current load case
|
||||
type(tBoundaryCondition), intent(in) :: &
|
||||
stress_BC, &
|
||||
deformation_BC
|
||||
real(pReal), dimension(3,3), intent(in) :: rotation_BC
|
||||
logical, intent(in) :: &
|
||||
guess
|
||||
PetscErrorCode :: ierr
|
||||
PetscScalar, pointer :: F(:,:,:,:)
|
||||
implicit none
|
||||
logical, intent(in) :: &
|
||||
guess
|
||||
real(pReal), intent(in) :: &
|
||||
timeinc_old, &
|
||||
timeinc, &
|
||||
loadCaseTime !< remaining time of current load case
|
||||
type(tBoundaryCondition), intent(in) :: &
|
||||
stress_BC, &
|
||||
deformation_BC
|
||||
real(pReal), dimension(3,3), intent(in) ::&
|
||||
rotation_BC
|
||||
PetscErrorCode :: ierr
|
||||
PetscScalar, pointer :: F(:,:,:,:)
|
||||
|
||||
character(len=1024) :: rankStr
|
||||
character(len=32) :: rankStr
|
||||
|
||||
call DMDAVecGetArrayF90(da,solution_vec,F,ierr)
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! restart information for spectral solver
|
||||
if (restartWrite) then
|
||||
write(6,'(/,a)') ' writing converged results for restart'
|
||||
flush(6)
|
||||
write(rankStr,'(a1,i0)')'_',worldrank
|
||||
call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file
|
||||
write (777,rec=1) F
|
||||
close (777)
|
||||
call IO_write_jobRealFile(777,'F_lastInc'//trim(rankStr),size(F_lastInc)) ! writing F_lastInc field to file
|
||||
write (777,rec=1) F_lastInc
|
||||
close (777)
|
||||
if (worldrank == 0_pInt) then
|
||||
call IO_write_jobRealFile(777,'F_aimDot',size(F_aimDot))
|
||||
write (777,rec=1) F_aimDot
|
||||
close(777)
|
||||
call IO_write_jobRealFile(777,'C_volAvg',size(C_volAvg))
|
||||
write (777,rec=1) C_volAvg
|
||||
close(777)
|
||||
call IO_write_jobRealFile(777,'C_volAvgLastInc',size(C_volAvgLastInc))
|
||||
write (777,rec=1) C_volAvgLastInc
|
||||
close(777)
|
||||
endif
|
||||
endif
|
||||
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
|
||||
|
||||
if (cutBack) then
|
||||
C_volAvg = C_volAvgLastInc ! QUESTION: where is this required?
|
||||
C_minMaxAvg = C_minMaxAvgLastInc ! QUESTION: where is this required?
|
||||
else
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! restart information for spectral solver
|
||||
if (restartWrite) then ! QUESTION: where is this logical properly set?
|
||||
write(6,'(/,a)') ' writing converged results for restart'
|
||||
flush(6)
|
||||
|
||||
call utilities_updateIPcoords(F)
|
||||
if (worldrank == 0_pInt) then
|
||||
call IO_write_jobRealFile(777,'C_volAvg',size(C_volAvg))
|
||||
write (777,rec=1) C_volAvg; close(777)
|
||||
call IO_write_jobRealFile(777,'C_volAvgLastInc',size(C_volAvgLastInc))
|
||||
write (777,rec=1) C_volAvgLastInc; close(777)
|
||||
call IO_write_jobRealFile(777,'C_minMaxAvg',size(C_volAvg))
|
||||
write (777,rec=1) C_minMaxAvg; close(777)
|
||||
call IO_write_jobRealFile(777,'C_minMaxAvgLastInc',size(C_volAvgLastInc))
|
||||
write (777,rec=1) C_minMaxAvgLastInc; close(777)
|
||||
endif
|
||||
|
||||
if (cutBack) then
|
||||
F_aim = F_aim_lastInc
|
||||
F = reshape(F_lastInc, [9,grid(1),grid(2),grid3])
|
||||
C_volAvg = C_volAvgLastInc
|
||||
else
|
||||
ForwardData = .True.
|
||||
C_volAvgLastInc = C_volAvg
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! calculate rate for aim
|
||||
if (deformation_BC%myType=='l') then ! calculate f_aimDot from given L and current F
|
||||
f_aimDot = deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim)
|
||||
elseif(deformation_BC%myType=='fdot') then ! f_aimDot is prescribed
|
||||
f_aimDot = deformation_BC%maskFloat * deformation_BC%values
|
||||
elseif(deformation_BC%myType=='f') then ! aim at end of load case is prescribed
|
||||
f_aimDot = deformation_BC%maskFloat * (deformation_BC%values -F_aim)/loadCaseTime
|
||||
endif
|
||||
if (guess) f_aimDot = f_aimDot + stress_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old
|
||||
F_aim_lastInc = F_aim
|
||||
write(rankStr,'(a1,i0)')'_',worldrank
|
||||
call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file
|
||||
write (777,rec=1) F; close (777)
|
||||
call IO_write_jobRealFile(777,'F_lastInc'//trim(rankStr),size(F_lastInc)) ! writing F_lastInc field to file
|
||||
write (777,rec=1) F_lastInc; close (777)
|
||||
endif
|
||||
|
||||
call CPFEM_age() ! age state and kinematics
|
||||
call utilities_updateIPcoords(F)
|
||||
|
||||
C_volAvgLastInc = C_volAvg
|
||||
C_minMaxAvgLastInc = C_minMaxAvg
|
||||
|
||||
if (guess) then ! QUESTION: better with a = L ? x:y
|
||||
F_aimDot = stress_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old ! initialize with correction based on last inc
|
||||
else
|
||||
F_aimDot = 0.0_pReal
|
||||
endif
|
||||
F_aim_lastInc = F_aim
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! calculate rate for aim
|
||||
if (deformation_BC%myType=='l') then ! calculate f_aimDot from given L and current F
|
||||
F_aimDot = &
|
||||
F_aimDot + deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim_lastInc)
|
||||
elseif(deformation_BC%myType=='fdot') then ! f_aimDot is prescribed
|
||||
F_aimDot = &
|
||||
F_aimDot + deformation_BC%maskFloat * deformation_BC%values
|
||||
elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed
|
||||
F_aimDot = &
|
||||
F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime
|
||||
endif
|
||||
|
||||
|
||||
Fdot = Utilities_calculateRate(guess, &
|
||||
F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, &
|
||||
math_rotate_backward33(f_aimDot,rotation_BC))
|
||||
F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) ! winding F forward
|
||||
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
|
||||
endif
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! update coordinates and rate and forward last inc
|
||||
call utilities_updateIPcoords(F)
|
||||
Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
|
||||
timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]))
|
||||
F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3])
|
||||
endif
|
||||
|
||||
F_aim = F_aim + f_aimDot * timeinc
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! update local deformation gradient
|
||||
F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! ensure that it matches rotated F_aim
|
||||
math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid3])
|
||||
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
|
||||
|
||||
! update average and local deformation gradients
|
||||
F_aim = F_aim_lastInc + f_aimDot * timeinc
|
||||
F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average
|
||||
math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid3])
|
||||
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
|
||||
|
||||
end subroutine BasicPETSc_forward
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
|
|||
|
|
@ -213,8 +213,8 @@ subroutine Polarisation_init
|
|||
endif restart
|
||||
|
||||
call Utilities_updateIPcoords(reshape(F,shape(F_lastInc)))
|
||||
call Utilities_constitutiveResponse(F_lastInc, reshape(F,shape(F_lastInc)), &
|
||||
0.0_pReal,P,C_volAvg,C_minMaxAvg,temp33_Real,.false.,math_I3)
|
||||
call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, &
|
||||
reshape(F,shape(F_lastInc)),0.0_pReal,math_I3)
|
||||
nullify(F)
|
||||
nullify(F_tau)
|
||||
call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! write data back to PETSc
|
||||
|
|
@ -364,12 +364,10 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
|
|||
DMDALocalInfo, dimension(&
|
||||
DMDA_LOCAL_INFO_SIZE) :: &
|
||||
in
|
||||
PetscScalar, target, dimension(3,3,2, &
|
||||
XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: &
|
||||
x_scal
|
||||
PetscScalar, target, dimension(3,3,2, &
|
||||
X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: &
|
||||
f_scal
|
||||
PetscScalar, &
|
||||
target, dimension(3,3,2, XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: x_scal
|
||||
PetscScalar, &
|
||||
target, dimension(3,3,2, X_RANGE, Y_RANGE, Z_RANGE), intent(out) :: f_scal
|
||||
PetscScalar, pointer, dimension(:,:,:,:,:) :: &
|
||||
F, &
|
||||
F_tau, &
|
||||
|
|
@ -440,8 +438,8 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
! evaluate constitutive response
|
||||
P_avLastEval = P_av
|
||||
call Utilities_constitutiveResponse(F_lastInc,F - residual_F_tau/polarBeta,params%timeinc, &
|
||||
residual_F,C_volAvg,C_minMaxAvg,P_av,ForwardData,params%rotation_BC)
|
||||
call Utilities_constitutiveResponse(residual_F,P_av,C_volAvg,C_minMaxAvg, &
|
||||
F - residual_F_tau/polarBeta,params%timeinc,params%rotation_BC)
|
||||
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
|
||||
ForwardData = .False.
|
||||
|
||||
|
|
@ -654,13 +652,13 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,deformati
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
! update coordinates and rate and forward last inc
|
||||
call utilities_updateIPcoords(F)
|
||||
Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
|
||||
timeinc_old,guess,F_lastInc, &
|
||||
reshape(F,[3,3,grid(1),grid(2),grid3]))
|
||||
F_tauDot = Utilities_calculateRate(math_rotate_backward33(2.0_pReal*f_aimDot,rotation_BC), &
|
||||
timeinc_old,guess,F_tau_lastInc, &
|
||||
reshape(F_tau,[3,3,grid(1),grid(2),grid3]))
|
||||
F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3])
|
||||
Fdot = Utilities_calculateRate(guess, &
|
||||
F_lastInc, reshape(F, [3,3,grid(1),grid(2),grid3]), timeinc_old, &
|
||||
math_rotate_backward33( f_aimDot,rotation_BC))
|
||||
F_tauDot = Utilities_calculateRate(guess, &
|
||||
F_tau_lastInc, reshape(F_tau,[3,3,grid(1),grid(2),grid3]), timeinc_old, &
|
||||
math_rotate_backward33(2.0_pReal*f_aimDot,rotation_BC))
|
||||
F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3])
|
||||
F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3])
|
||||
endif
|
||||
|
||||
|
|
|
|||
|
|
@ -16,7 +16,7 @@ module spectral_utilities
|
|||
#include <petsc/finclude/petscsys.h>
|
||||
include 'fftw3-mpi.f03'
|
||||
|
||||
logical, public :: cutBack =.false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill
|
||||
logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill
|
||||
integer(pInt), public, parameter :: maxPhaseFields = 2_pInt
|
||||
integer(pInt), public :: nActiveFields = 0_pInt
|
||||
|
||||
|
|
@ -799,7 +799,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
|
|||
|
||||
call math_invert(size_reduced, c_reduced, s_reduced, errmatinv) ! invert reduced stiffness
|
||||
if (any(IEEE_is_NaN(s_reduced))) errmatinv = .true.
|
||||
if(errmatinv) call IO_error(error_ID=400_pInt,ext_msg='utilities_maskedCompliance')
|
||||
if (errmatinv) call IO_error(error_ID=400_pInt,ext_msg='utilities_maskedCompliance')
|
||||
temp99_Real = 0.0_pReal ! fill up compliance with zeros
|
||||
k = 0_pInt
|
||||
do n = 1_pInt,9_pInt
|
||||
|
|
@ -817,28 +817,30 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
|
|||
sTimesC = matmul(c_reduced,s_reduced)
|
||||
do m=1_pInt, size_reduced
|
||||
do n=1_pInt, size_reduced
|
||||
if(m==n .and. abs(sTimesC(m,n)) > (1.0_pReal + 10.0e-12_pReal)) errmatinv = .true. ! diagonal elements of S*C should be 1
|
||||
if(m/=n .and. abs(sTimesC(m,n)) > (0.0_pReal + 10.0e-12_pReal)) errmatinv = .true. ! off diagonal elements of S*C should be 0
|
||||
errmatinv = errmatinv &
|
||||
.or. (m==n .and. abs(sTimesC(m,n)-1.0_pReal) > 1.0e-12_pReal) & ! diagonal elements of S*C should be 1
|
||||
.or. (m/=n .and. abs(sTimesC(m,n)) > 1.0e-12_pReal) ! off-diagonal elements of S*C should be 0
|
||||
enddo
|
||||
enddo
|
||||
if(debugGeneral .or. errmatinv) then
|
||||
write(formatString, '(I16.16)') size_reduced
|
||||
if (debugGeneral .or. errmatinv) then
|
||||
write(formatString, '(i2)') size_reduced
|
||||
formatString = '(/,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))'
|
||||
write(6,trim(formatString),advance='no') ' C * S (load) ', &
|
||||
transpose(matmul(c_reduced,s_reduced))
|
||||
write(6,trim(formatString),advance='no') ' S (load) ', transpose(s_reduced)
|
||||
if(errmatinv) call IO_error(error_ID=400_pInt,ext_msg='utilities_maskedCompliance')
|
||||
endif
|
||||
if(errmatinv) call IO_error(error_ID=400_pInt,ext_msg='utilities_maskedCompliance')
|
||||
deallocate(c_reduced)
|
||||
deallocate(s_reduced)
|
||||
deallocate(sTimesC)
|
||||
else
|
||||
temp99_real = 0.0_pReal
|
||||
endif
|
||||
if(debugGeneral) &
|
||||
write(6,'(/,a,/,9(9(2x,f12.7,1x)/),/)',advance='no') ' Masked Compliance (load) * GPa =', &
|
||||
transpose(temp99_Real*1.e9_pReal)
|
||||
flush(6)
|
||||
if(debugGeneral) then
|
||||
write(6,'(/,a,/,9(9(2x,f10.5,1x)/),/)',advance='no') &
|
||||
' Masked Compliance (load) / GPa =', transpose(temp99_Real*1.e-9_pReal)
|
||||
flush(6)
|
||||
endif
|
||||
utilities_maskedCompliance = math_Plain99to3333(temp99_Real)
|
||||
|
||||
end function utilities_maskedCompliance
|
||||
|
|
@ -924,10 +926,10 @@ end subroutine utilities_fourierTensorDivergence
|
|||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief calculates constitutive response
|
||||
!> @brief calculate constitutive response from materialpoint_F0 to F during timeinc
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
|
||||
P,C_volAvg,C_minmaxAvg,P_av,forwardData,rotation_BC)
|
||||
subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
|
||||
F,timeinc,rotation_BC)
|
||||
use IO, only: &
|
||||
IO_error
|
||||
use debug, only: &
|
||||
|
|
@ -940,31 +942,22 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
|
|||
use mesh, only: &
|
||||
grid,&
|
||||
grid3
|
||||
use FEsolving, only: &
|
||||
restartWrite
|
||||
use CPFEM2, only: &
|
||||
CPFEM_general
|
||||
use homogenization, only: &
|
||||
materialpoint_F0, &
|
||||
materialpoint_F, &
|
||||
materialpoint_P, &
|
||||
materialpoint_dPdF
|
||||
materialpoint_dPdF, &
|
||||
materialpoint_stressAndItsTangent
|
||||
|
||||
implicit none
|
||||
real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: &
|
||||
F_lastInc, & !< target deformation gradient
|
||||
F !< previous deformation gradient
|
||||
real(pReal), intent(in) :: timeinc !< loading time
|
||||
logical, intent(in) :: forwardData !< age results
|
||||
real(pReal), intent(in), dimension(3,3) :: rotation_BC !< rotation of load frame
|
||||
|
||||
real(pReal),intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness
|
||||
real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress
|
||||
real(pReal),intent(out), dimension(3,3,grid(1),grid(2),grid3) :: P !< PK stress
|
||||
|
||||
logical :: &
|
||||
age
|
||||
real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: F !< deformation gradient target !< previous deformation gradient
|
||||
real(pReal), intent(in) :: timeinc !< loading time
|
||||
real(pReal), intent(in), dimension(3,3) :: rotation_BC !< rotation of load frame
|
||||
|
||||
|
||||
integer(pInt) :: &
|
||||
j,k,ierr
|
||||
real(pReal), dimension(3,3,3,3) :: max_dPdF, min_dPdF
|
||||
|
|
@ -975,17 +968,9 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
|
|||
|
||||
write(6,'(/,a)') ' ... evaluating constitutive response ......................................'
|
||||
flush(6)
|
||||
age = .False.
|
||||
|
||||
if (forwardData) then ! aging results
|
||||
age = .True.
|
||||
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3])
|
||||
endif
|
||||
if (cutBack) age = .False. ! restore saved variables
|
||||
|
||||
materialpoint_F = reshape(F,[3,3,1,product(grid(1:2))*grid3])
|
||||
call debug_reset() ! this has no effect on rank >0
|
||||
|
||||
|
||||
materialpoint_F = reshape(F,[3,3,1,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! calculate bounds of det(F) and report
|
||||
if(debugGeneral) then
|
||||
|
|
@ -1002,7 +987,19 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
|
|||
flush(6)
|
||||
endif
|
||||
|
||||
call CPFEM_general(age,timeinc)
|
||||
call debug_reset() ! this has no effect on rank >0
|
||||
call materialpoint_stressAndItsTangent(.true.,timeinc) ! calculate P field
|
||||
|
||||
P = reshape(materialpoint_P, [3,3,grid(1),grid(2),grid3])
|
||||
P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P
|
||||
call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
|
||||
if (debugRotation) &
|
||||
write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress (lab) / MPa =',&
|
||||
math_transpose33(P_av)*1.e-6_pReal
|
||||
P_av = math_rotate_forward33(P_av,rotation_BC)
|
||||
write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',&
|
||||
math_transpose33(P_av)*1.e-6_pReal
|
||||
flush(6)
|
||||
|
||||
max_dPdF = 0.0_pReal
|
||||
max_dPdF_norm = 0.0_pReal
|
||||
|
|
@ -1020,38 +1017,24 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
|
|||
end do
|
||||
|
||||
call MPI_Allreduce(MPI_IN_PLACE,max_dPdF,81,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr)
|
||||
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce max')
|
||||
if (ierr /= 0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce max')
|
||||
call MPI_Allreduce(MPI_IN_PLACE,min_dPdF,81,MPI_DOUBLE,MPI_MIN,PETSC_COMM_WORLD,ierr)
|
||||
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce min')
|
||||
if (ierr /= 0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce min')
|
||||
|
||||
C_minmaxAvg = 0.5_pReal*(max_dPdF + min_dPdF)
|
||||
C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) * wgt
|
||||
|
||||
C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) * wgt
|
||||
call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
|
||||
|
||||
call debug_info() ! this has no effect on rank >0
|
||||
|
||||
restartWrite = .false. ! reset restartWrite status
|
||||
cutBack = .false. ! reset cutBack status
|
||||
|
||||
P = reshape(materialpoint_P, [3,3,grid(1),grid(2),grid3])
|
||||
P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P
|
||||
call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
|
||||
if (debugRotation) &
|
||||
write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress (lab) / MPa =',&
|
||||
math_transpose33(P_av)*1.e-6_pReal
|
||||
P_av = math_rotate_forward33(P_av,rotation_BC)
|
||||
write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',&
|
||||
math_transpose33(P_av)*1.e-6_pReal
|
||||
flush(6)
|
||||
|
||||
end subroutine utilities_constitutiveResponse
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief calculates forward rate, either guessing or just add delta/timeinc
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
pure function utilities_calculateRate(avRate,timeinc_old,guess,field_lastInc,field)
|
||||
pure function utilities_calculateRate(heterogeneous,field0,field,dt,avRate)
|
||||
use mesh, only: &
|
||||
grid3, &
|
||||
grid
|
||||
|
|
@ -1059,17 +1042,17 @@ pure function utilities_calculateRate(avRate,timeinc_old,guess,field_lastInc,fie
|
|||
implicit none
|
||||
real(pReal), intent(in), dimension(3,3) :: avRate !< homogeneous addon
|
||||
real(pReal), intent(in) :: &
|
||||
timeinc_old !< timeinc of last step
|
||||
dt !< timeinc between field0 and field
|
||||
logical, intent(in) :: &
|
||||
guess !< guess along former trajectory
|
||||
heterogeneous !< calculate field of rates
|
||||
real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: &
|
||||
field_lastInc, & !< data of previous step
|
||||
field0, & !< data of previous step
|
||||
field !< data of current step
|
||||
real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: &
|
||||
utilities_calculateRate
|
||||
|
||||
if (guess) then
|
||||
utilities_calculateRate = (field-field_lastInc) / timeinc_old
|
||||
if (heterogeneous) then
|
||||
utilities_calculateRate = (field-field0) / dt
|
||||
else
|
||||
utilities_calculateRate = spread(spread(spread(avRate,3,grid(1)),4,grid(2)),5,grid3)
|
||||
endif
|
||||
|
|
|
|||
Loading…
Reference in New Issue