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

included changes to correct cutback issue of spectral solver

This commit is contained in:
Philip Eisenlohr 2018-02-16 09:36:18 -05:00
parent 15b0fa90a1
commit 0a763ff116
19 changed files with 952 additions and 1179 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

0
examples/SpectralMethod/20grains.seeds → examples/SpectralMethod/Polycrystal/20grains.seeds
View File

0
examples/SpectralMethod/20grains16x16x16.geom → examples/SpectralMethod/Polycrystal/20grains16x16x16.geom
View File

0
examples/SpectralMethod/20grains32x32x32.geom → examples/SpectralMethod/Polycrystal/20grains32x32x32.geom
View File

0
examples/SpectralMethod/20grains64x64x64.geom → examples/SpectralMethod/Polycrystal/20grains64x64x64.geom
View File

0
examples/SpectralMethod/material.config → examples/SpectralMethod/Polycrystal/material.config
View File

0
examples/SpectralMethod/numerics.config → examples/SpectralMethod/Polycrystal/numerics.config
View File

0
examples/SpectralMethod/shearXY.load → examples/SpectralMethod/Polycrystal/shearXY.load
View File

0
examples/SpectralMethod/shearZX.load → examples/SpectralMethod/Polycrystal/shearZX.load
View File

0
examples/SpectralMethod/tensionX.load → examples/SpectralMethod/Polycrystal/tensionX.load
View File

182
src/CPFEM2.f90
View File

@ -9,7 +9,7 @@ module CPFEM2
private private
public :: & public :: &
CPFEM_general, & CPFEM_age, &
CPFEM_initAll CPFEM_initAll
contains contains
@ -127,6 +127,7 @@ subroutine CPFEM_init
write(6,'(/,a)') ' <<<+- CPFEM init -+>>>' write(6,'(/,a)') ' <<<+- CPFEM init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
flush(6)
endif mainProcess endif mainProcess
! *** restore the last converged values of each essential variable from the binary file ! *** 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 write(rankStr,'(a1,i0)')'_',worldrank
call IO_read_intFile(777,'recordedPhase'//trim(rankStr),modelName,size(material_phase)) call IO_read_intFile(777,'recordedPhase'//trim(rankStr),modelName,size(material_phase))
read (777,rec=1) material_phase read (777,rec=1) material_phase; close (777)
close (777)
call IO_read_realFile(777,'convergedF'//trim(rankStr),modelName,size(crystallite_F0)) call IO_read_realFile(777,'convergedF'//trim(rankStr),modelName,size(crystallite_F0))
read (777,rec=1) crystallite_F0 read (777,rec=1) crystallite_F0; close (777)
close (777)
call IO_read_realFile(777,'convergedFp'//trim(rankStr),modelName,size(crystallite_Fp0)) call IO_read_realFile(777,'convergedFp'//trim(rankStr),modelName,size(crystallite_Fp0))
read (777,rec=1) crystallite_Fp0 read (777,rec=1) crystallite_Fp0; close (777)
close (777)
call IO_read_realFile(777,'convergedFi'//trim(rankStr),modelName,size(crystallite_Fi0)) call IO_read_realFile(777,'convergedFi'//trim(rankStr),modelName,size(crystallite_Fi0))
read (777,rec=1) crystallite_Fi0 read (777,rec=1) crystallite_Fi0; close (777)
close (777)
call IO_read_realFile(777,'convergedLp'//trim(rankStr),modelName,size(crystallite_Lp0)) call IO_read_realFile(777,'convergedLp'//trim(rankStr),modelName,size(crystallite_Lp0))
read (777,rec=1) crystallite_Lp0 read (777,rec=1) crystallite_Lp0; close (777)
close (777)
call IO_read_realFile(777,'convergedLi'//trim(rankStr),modelName,size(crystallite_Li0)) call IO_read_realFile(777,'convergedLi'//trim(rankStr),modelName,size(crystallite_Li0))
read (777,rec=1) crystallite_Li0 read (777,rec=1) crystallite_Li0; close (777)
close (777)
call IO_read_realFile(777,'convergeddPdF'//trim(rankStr),modelName,size(crystallite_dPdF0)) call IO_read_realFile(777,'convergeddPdF'//trim(rankStr),modelName,size(crystallite_dPdF0))
read (777,rec=1) crystallite_dPdF0 read (777,rec=1) crystallite_dPdF0; close (777)
close (777)
call IO_read_realFile(777,'convergedTstar'//trim(rankStr),modelName,size(crystallite_Tstar0_v)) call IO_read_realFile(777,'convergedTstar'//trim(rankStr),modelName,size(crystallite_Tstar0_v))
read (777,rec=1) crystallite_Tstar0_v read (777,rec=1) crystallite_Tstar0_v; close (777)
close (777)
call IO_read_realFile(777,'convergedStateConst'//trim(rankStr),modelName) call IO_read_realFile(777,'convergedStateConst'//trim(rankStr),modelName)
m = 0_pInt m = 0_pInt
@ -194,7 +187,6 @@ subroutine CPFEM_init
restartRead = .false. restartRead = .false.
endif endif
flush(6)
end subroutine CPFEM_init 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 !> @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: & use prec, only: &
pReal, & pReal, &
pInt pInt
@ -215,7 +207,6 @@ subroutine CPFEM_general(age, dt)
debug_levelExtensive, & debug_levelExtensive, &
debug_levelSelective debug_levelSelective
use FEsolving, only: & use FEsolving, only: &
terminallyIll, &
restartWrite restartWrite
use math, only: & use math, only: &
math_identity2nd, & math_identity2nd, &
@ -254,114 +245,99 @@ subroutine CPFEM_general(age, dt)
crystallite_dPdF, & crystallite_dPdF, &
crystallite_Tstar0_v, & crystallite_Tstar0_v, &
crystallite_Tstar_v crystallite_Tstar_v
use homogenization, only: &
materialpoint_stressAndItsTangent, &
materialpoint_postResults
use IO, only: & use IO, only: &
IO_write_jobRealFile, & IO_write_jobRealFile, &
IO_warning IO_warning
use DAMASK_interface use DAMASK_interface
implicit none implicit none
real(pReal), intent(in) :: dt !< time increment
logical, intent(in) :: age !< age results
integer(pInt) :: i, k, l, m, ph, homog, mySource 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 (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) &
if (age) then 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
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 crystallite_F0 = crystallite_partionedF ! crystallite deformation (_subF is perturbed...)
do i = 1, size(sourceState) crystallite_Fp0 = crystallite_Fp ! crystallite plastic deformation
do mySource = 1,phase_Nsources(i) crystallite_Lp0 = crystallite_Lp ! crystallite plastic velocity
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 crystallite_Fi0 = crystallite_Fi ! crystallite intermediate deformation
enddo; enddo crystallite_Li0 = crystallite_Li ! crystallite intermediate velocity
if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) & crystallite_dPdF0 = crystallite_dPdF ! crystallite stiffness
write(6,'(a)') '<< CPFEM >> aging states' crystallite_Tstar0_v = crystallite_Tstar_v ! crystallite 2nd Piola Kirchhoff stress
do homog = 1_pInt, material_Nhomogenization 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
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
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 do homog = 1_pInt, material_Nhomogenization
if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) & homogState (homog)%state0 = homogState (homog)%state
write(6,'(a)') '<< CPFEM >> writing state variables of last converged step to binary files' thermalState (homog)%state0 = thermalState (homog)%state
damageState (homog)%state0 = damageState (homog)%state
write(rankStr,'(a1,i0)')'_',worldrank vacancyfluxState (homog)%state0 = vacancyfluxState (homog)%state
hydrogenfluxState(homog)%state0 = hydrogenfluxState(homog)%state
enddo
call IO_write_jobRealFile(777,'recordedPhase'//trim(rankStr),size(material_phase)) if (restartWrite) then
write (777,rec=1) material_phase if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) &
close (777) 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(rankStr,'(a1,i0)')'_',worldrank
write (777,rec=1) crystallite_F0
close (777)
call IO_write_jobRealFile(777,'convergedFp'//trim(rankStr),size(crystallite_Fp0)) call IO_write_jobRealFile(777,'recordedPhase'//trim(rankStr),size(material_phase))
write (777,rec=1) crystallite_Fp0 write (777,rec=1) material_phase; close (777)
close (777)
call IO_write_jobRealFile(777,'convergedFi'//trim(rankStr),size(crystallite_Fi0)) call IO_write_jobRealFile(777,'convergedF'//trim(rankStr),size(crystallite_F0))
write (777,rec=1) crystallite_Fi0 write (777,rec=1) crystallite_F0; close (777)
close (777)
call IO_write_jobRealFile(777,'convergedLp'//trim(rankStr),size(crystallite_Lp0)) call IO_write_jobRealFile(777,'convergedFp'//trim(rankStr),size(crystallite_Fp0))
write (777,rec=1) crystallite_Lp0 write (777,rec=1) crystallite_Fp0; close (777)
close (777)
call IO_write_jobRealFile(777,'convergedLi'//trim(rankStr),size(crystallite_Li0)) call IO_write_jobRealFile(777,'convergedFi'//trim(rankStr),size(crystallite_Fi0))
write (777,rec=1) crystallite_Li0 write (777,rec=1) crystallite_Fi0; close (777)
close (777)
call IO_write_jobRealFile(777,'convergeddPdF'//trim(rankStr),size(crystallite_dPdF0)) call IO_write_jobRealFile(777,'convergedLp'//trim(rankStr),size(crystallite_Lp0))
write (777,rec=1) crystallite_dPdF0 write (777,rec=1) crystallite_Lp0; close (777)
close (777)
call IO_write_jobRealFile(777,'convergedTstar'//trim(rankStr),size(crystallite_Tstar0_v)) call IO_write_jobRealFile(777,'convergedLi'//trim(rankStr),size(crystallite_Li0))
write (777,rec=1) crystallite_Tstar0_v write (777,rec=1) crystallite_Li0; close (777)
close (777)
call IO_write_jobRealFile(777,'convergedStateConst'//trim(rankStr)) call IO_write_jobRealFile(777,'convergeddPdF'//trim(rankStr),size(crystallite_dPdF0))
m = 0_pInt write (777,rec=1) crystallite_dPdF0; close (777)
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,'convergedStateHomog'//trim(rankStr)) call IO_write_jobRealFile(777,'convergedTstar'//trim(rankStr),size(crystallite_Tstar0_v))
m = 0_pInt write (777,rec=1) crystallite_Tstar0_v; close (777)
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)
endif call IO_write_jobRealFile(777,'convergedStateConst'//trim(rankStr))
endif 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 IO_write_jobRealFile(777,'convergedStateHomog'//trim(rankStr))
call materialpoint_stressAndItsTangent(.True., dt) 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 end module CPFEM2

258
src/DAMASK_spectral.f90 Normal file → Executable file
View File

@ -78,8 +78,7 @@ program DAMASK_spectral
FIELD_UNDEFINED_ID, & FIELD_UNDEFINED_ID, &
FIELD_MECH_ID, & FIELD_MECH_ID, &
FIELD_THERMAL_ID, & FIELD_THERMAL_ID, &
FIELD_DAMAGE_ID, & FIELD_DAMAGE_ID
utilities_calcPlasticity
use spectral_mech_Basic use spectral_mech_Basic
use spectral_mech_AL use spectral_mech_AL
use spectral_mech_Polarisation use spectral_mech_Polarisation
@ -157,19 +156,6 @@ program DAMASK_spectral
MPI_finalize, & MPI_finalize, &
MPI_allreduce, & MPI_allreduce, &
PETScFinalize 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) ! init DAMASK (all modules)
@ -227,8 +213,6 @@ program DAMASK_spectral
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! reading the load case and assign values to the allocated data structure ! reading the load case and assign values to the allocated data structure
yieldStop = .False.
yieldStopSatisfied = .False.
rewind(FILEUNIT) rewind(FILEUNIT)
do do
line = IO_read(FILEUNIT) line = IO_read(FILEUNIT)
@ -303,30 +287,10 @@ program DAMASK_spectral
temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j)
enddo enddo
loadCases(currentLoadCase)%rotation = math_plain9to33(temp_valueVector) 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 end select
enddo; enddo enddo; enddo
close(FILEUNIT) 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 ! consistency checks and output of load case
loadCases(1)%followFormerTrajectory = .false. ! cannot guess along trajectory for first inc of first currentLoadCase 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 if (loadCases(currentLoadCase)%outputfrequency < 1_pInt) errorID = 836_pInt ! non-positive result frequency
write(6,'(2x,a,i5)') 'output frequency: ', & write(6,'(2x,a,i5)') 'output frequency: ', &
loadCases(currentLoadCase)%outputfrequency loadCases(currentLoadCase)%outputfrequency
write(6,'(2x,a,i5,/)') 'restart frequency: ', & write(6,'(2x,a,i5,/)') 'restart frequency: ', &
loadCases(currentLoadCase)%restartfrequency loadCases(currentLoadCase)%restartfrequency
if (errorID > 0_pInt) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message if (errorID > 0_pInt) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message
enddo checkLoadcases 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 (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_seek')
if (.not. appendToOutFile) then ! if not restarting, write 0th increment if (.not. appendToOutFile) then ! if not restarting, write 0th increment
write(6,'(1/,a)') ' ... writing initial configuration to file ........................'
do i = 1, size(materialpoint_results,3)/(maxByteOut/(materialpoint_sizeResults*pReal))+1 ! slice the output of my process in chunks not exceeding the limit for one output do i = 1, size(materialpoint_results,3)/(maxByteOut/(materialpoint_sizeResults*pReal))+1 ! slice the output of my process in chunks not exceeding the limit for one output
outputIndex = int([(i-1_pInt)*((maxRealOut)/materialpoint_sizeResults)+1_pInt, & outputIndex = int([(i-1_pInt)*((maxRealOut)/materialpoint_sizeResults)+1_pInt, & ! QUESTION: why not starting i at 0 instead of murky 1?
min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt) min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt)
call MPI_file_write(resUnit, & call MPI_file_write(resUnit, &
reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)), & reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)), &
@ -489,24 +454,23 @@ program DAMASK_spectral
if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_write') if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_write')
enddo enddo
fileOffset = fileOffset + sum(outputSize) ! forward to current file position fileOffset = fileOffset + sum(outputSize) ! forward to current file position
write(6,'(1/,a)') ' ... writing initial configuration to file ........................'
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! loopping over loadcases ! looping over loadcases
loadCaseLooping: do currentLoadCase = 1_pInt, size(loadCases) loadCaseLooping: do currentLoadCase = 1_pInt, size(loadCases)
time0 = time ! currentLoadCase start time time0 = time ! currentLoadCase start time
guess = loadCases(currentLoadCase)%followFormerTrajectory ! change of load case? homogeneous guess for the first inc 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 incLooping: do inc = 1_pInt, loadCases(currentLoadCase)%incs
totalIncsCounter = totalIncsCounter + 1_pInt totalIncsCounter = totalIncsCounter + 1_pInt
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! forwarding time ! forwarding time
timeIncOld = timeinc timeIncOld = timeinc ! last timeinc that brought former inc to an end
if (loadCases(currentLoadCase)%logscale == 0_pInt) then ! linear scale 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 else
if (currentLoadCase == 1_pInt) then ! 1st currentLoadCase of logarithmic scale if (currentLoadCase == 1_pInt) then ! 1st currentLoadCase of logarithmic scale
if (inc == 1_pInt) then ! 1st inc of 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)/& ( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc,pReal)/&
real(loadCases(currentLoadCase)%incs ,pReal))& real(loadCases(currentLoadCase)%incs ,pReal))&
-(1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( (inc-1_pInt),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
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 skipping: if (totalIncsCounter < restartInc) then ! not yet at restart inc?
stepFraction = 0_pInt time = time + timeinc ! just advance time, skip already performed calculation
guess = .true. ! QUESTION:why forced guessing instead of inheriting loadcase preference
else skipping
stepFraction = 0_pInt ! fraction scaled by stepFactor**cutLevel
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! loop over sub incs ! loop over sub step
subIncLooping: do while (stepFraction/subStepFactor**cutBackLevel <1_pInt) subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel)
time = time + timeinc ! forward time remainingLoadCaseTime = loadCases(currentLoadCase)%time+time0 - time
stepFraction = stepFraction + 1_pInt time = time + timeinc ! forward target time
remainingLoadCaseTime = time0 - time + loadCases(currentLoadCase)%time + timeInc stepFraction = stepFraction + 1_pInt ! count step
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report begin of new increment ! report begin of new step
write(6,'(/,a)') ' ###########################################################################' write(6,'(/,a)') ' ###########################################################################'
write(6,'(1x,a,es12.5'//& write(6,'(1x,a,es12.5'//&
',a,'//IO_intOut(inc)//',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//& ',a,'//IO_intOut(inc) //',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//&
',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//& ',a,'//IO_intOut(stepFraction) //',a,'//IO_intOut(subStepFactor**cutBackLevel)//&
',a,'//IO_intOut(currentLoadCase)//',a,'//IO_intOut(size(loadCases))//')') & ',a,'//IO_intOut(currentLoadCase)//',a,'//IO_intOut(size(loadCases))//')') &
'Time', time, & 'Time', time, &
's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,& 's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,&
'-', stepFraction, '/', subStepFactor**cutBackLevel,& '-', stepFraction, '/', subStepFactor**cutBackLevel,&
' of load case ', currentLoadCase,'/',size(loadCases) ' 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) 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 ! forward fields
@ -578,7 +548,7 @@ program DAMASK_spectral
end select end select
case(FIELD_THERMAL_ID); call spectral_thermal_forward() 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 end select
enddo enddo
@ -618,65 +588,63 @@ program DAMASK_spectral
solres(field) = spectral_damage_solution(timeinc,timeIncOld,remainingLoadCaseTime) solres(field) = spectral_damage_solution(timeinc,timeIncOld,remainingLoadCaseTime)
end select end select
if (.not. solres(field)%converged) exit ! no solution found if (.not. solres(field)%converged) exit ! no solution found
enddo enddo
stagIter = stagIter + 1_pInt stagIter = stagIter + 1_pInt
stagIterate = stagIter < stagItMax .and. & stagIterate = stagIter < stagItMax &
all(solres(:)%converged) .and. & .and. all(solres(:)%converged) &
.not. all(solres(:)%stagConverged) .and. .not. all(solres(:)%stagConverged) ! stationary with respect to staggered iteration
enddo enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! check solution ! check solution for either advance or retry
cutBack = .False.
if(solres(1)%termIll .or. .not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found if ( (continueCalculation .or. all(solres(:)%converged .and. solres(:)%stagConverged)) & ! don't care or did converge
if (cutBackLevel < maxCutBack) then ! do cut back .and. .not. solres(1)%termIll) then ! and acceptable solution found
write(6,'(/,a)') ' cut back detected' timeIncOld = timeinc
cutBack = .True. cutBack = .false.
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
guess = .true. ! start guessing after first converged (sub)inc guess = .true. ! start guessing after first converged (sub)inc
endif
if (.not. cutBack) then
if (worldrank == 0) then if (worldrank == 0) then
write(statUnit,*) totalIncsCounter, time, cutBackLevel, & write(statUnit,*) totalIncsCounter, time, cutBackLevel, &
solres%converged, solres%iterationsNeeded ! write statistics about accepted solution solres%converged, solres%iterationsNeeded
flush(statUnit) flush(statUnit)
endif 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 endif
enddo subIncLooping
enddo subStepLooping
cutBackLevel = max(0_pInt, cutBackLevel - 1_pInt) ! try half number of subincs next inc cutBackLevel = max(0_pInt, cutBackLevel - 1_pInt) ! try half number of subincs next inc
if(all(solres(:)%converged)) then ! report converged inc
if (all(solres(:)%converged)) then
convergedCounter = convergedCounter + 1_pInt convergedCounter = convergedCounter + 1_pInt
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report converged inc
' increment ', totalIncsCounter, ' converged' ' increment ', totalIncsCounter, ' converged'
else else
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report non-converged inc
' increment ', totalIncsCounter, ' NOT converged'
notConvergedCounter = notConvergedCounter + 1_pInt notConvergedCounter = notConvergedCounter + 1_pInt
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report non-converged inc
' increment ', totalIncsCounter, ' NOT converged'
endif; flush(6) endif; flush(6)
if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency
if (worldrank == 0) & write(6,'(1/,a)') ' ... writing results to file ......................................'
write(6,'(1/,a)') ' ... writing results to file ......................................' flush(6)
call materialpoint_postResults() call materialpoint_postResults()
call MPI_file_seek (resUnit,fileOffset,MPI_SEEK_SET,ierr) 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 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, &
min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt) min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt)
@ -688,93 +656,29 @@ program DAMASK_spectral
enddo enddo
fileOffset = fileOffset + sum(outputSize) ! forward to current file position fileOffset = fileOffset + sum(outputSize) ! forward to current file position
endif endif
if( loadCases(currentLoadCase)%restartFrequency > 0_pInt .and. & ! at frequency of writing restart information set restart parameter for FEsolving if ( loadCases(currentLoadCase)%restartFrequency > 0_pInt & ! writing of restart info requested ...
mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0_pInt) then ! first call to CPFEM_general will write? .and. mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0_pInt) then ! ... and at frequency of writing restart information
restartWrite = .true. restartWrite = .true. ! set restart parameter for FEsolving
lastRestartWritten = inc lastRestartWritten = inc ! QUESTION: first call to CPFEM_general will write?
endif endif
else forwarding
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 endif skipping
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
enddo incLooping enddo incLooping
enddo loadCaseLooping enddo loadCaseLooping
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report summary of whole calculation ! report summary of whole calculation
write(6,'(/,a)') ' ###########################################################################' write(6,'(/,a)') ' ###########################################################################'
write(6,'(1x,i6.6,a,i6.6,a,f5.1,a)') convergedCounter, ' out of ', & write(6,'(1x,'//IO_intOut(convergedCounter)//',a,'//IO_intOut(notConvergedCounter + convergedCounter)//',a,f5.1,a)') &
notConvergedCounter + convergedCounter, ' (', & convergedCounter, ' out of ', &
real(convergedCounter, pReal)/& notConvergedCounter + convergedCounter, ' (', &
real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, & real(convergedCounter, pReal)/&
' %) increments converged!' real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, &
' %) increments converged!'
flush(6)
call MPI_file_close(resUnit,ierr) call MPI_file_close(resUnit,ierr)
close(statUnit) close(statUnit)
@ -831,8 +735,6 @@ subroutine quit(stop_id)
call PETScFinalize(ierr) call PETScFinalize(ierr)
if (ierr /= 0) write(6,'(a)') ' Error in PETScFinalize' if (ierr /= 0) write(6,'(a)') ' Error in PETScFinalize'
#ifdef _OPENMP #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) call MPI_finalize(error)
if (error /= 0) write(6,'(a)') ' Error in MPI_finalize' if (error /= 0) write(6,'(a)') ' Error in MPI_finalize'
#endif #endif

6
src/IO.f90
View File

@ -1268,11 +1268,11 @@ integer(pInt) function IO_countNumericalDataLines(fileUnit)
line = IO_read(fileUnit) line = IO_read(fileUnit)
chunkPos = IO_stringPos(line) chunkPos = IO_stringPos(line)
tmp = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) 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 line = IO_read(fileUnit, .true.) ! reset IO_read
exit exit
else
IO_countNumericalDataLines = IO_countNumericalDataLines + 1_pInt
endif endif
enddo enddo
backspace(fileUnit) backspace(fileUnit)

58
src/crystallite.f90 Normal file → Executable file
View File

@ -55,14 +55,14 @@ module crystallite
crystallite_Li0, & !< intermediate velocitiy grad at start of FE inc crystallite_Li0, & !< intermediate velocitiy grad at start of FE inc
crystallite_partionedLi0,& !< intermediate velocity grad at start of homog inc crystallite_partionedLi0,& !< intermediate velocity grad at start of homog inc
crystallite_Fe, & !< current "elastic" def grad (end of converged time step) crystallite_Fe, & !< current "elastic" def grad (end of converged time step)
crystallite_P, & !< 1st Piola-Kirchhoff stress per grain crystallite_P !< 1st Piola-Kirchhoff stress per grain
crystallite_subF !< def grad to be reached at end of crystallite inc
real(pReal), dimension(:,:,:,:,:), allocatable, private :: & real(pReal), dimension(:,:,:,:,:), allocatable, private :: &
crystallite_subFe0,& !< "elastic" def grad at start of crystallite inc crystallite_subFe0,& !< "elastic" def grad at start of crystallite inc
crystallite_invFp, & !< inverse of current plastic def grad (end of converged time step) crystallite_invFp, & !< inverse of current plastic def grad (end of converged time step)
crystallite_subFp0,& !< plastic def grad at start of crystallite inc crystallite_subFp0,& !< plastic def grad at start of crystallite inc
crystallite_invFi, & !< inverse of current intermediate def grad (end of converged time step) crystallite_invFi, & !< inverse of current intermediate def grad (end of converged time step)
crystallite_subFi0,& !< intermediate def grad at start of crystallite inc 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_subF0, & !< def grad at start of crystallite inc
crystallite_subLp0,& !< plastic velocity 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 crystallite_subLi0,& !< intermediate velocity grad at start of crystallite inc
@ -193,7 +193,7 @@ subroutine crystallite_init
c, & !< counter in integration point component loop c, & !< counter in integration point component loop
i, & !< counter in integration point loop i, & !< counter in integration point loop
e, & !< counter in element loop e, & !< counter in element loop
o, & !< counter in output loop o = 0_pInt, & !< counter in output loop
r, & !< counter in crystallite loop r, & !< counter in crystallite loop
cMax, & !< maximum number of integration point components cMax, & !< maximum number of integration point components
iMax, & !< maximum number of integration points iMax, & !< maximum number of integration points
@ -1239,14 +1239,16 @@ subroutine crystallite_integrateStateRK4()
use numerics, only: & use numerics, only: &
numerics_integrationMode numerics_integrationMode
use debug, only: & use debug, only: &
#ifdef DEBUG
debug_e, &
debug_i, &
debug_g, &
#endif
debug_level, & debug_level, &
debug_crystallite, & debug_crystallite, &
debug_levelBasic, & debug_levelBasic, &
debug_levelExtensive, & debug_levelExtensive, &
debug_levelSelective, & debug_levelSelective, &
debug_e, &
debug_i, &
debug_g, &
debug_StateLoopDistribution debug_StateLoopDistribution
use FEsolving, only: & use FEsolving, only: &
FEsolving_execElem, & FEsolving_execElem, &
@ -1533,14 +1535,16 @@ subroutine crystallite_integrateStateRKCK45()
use, intrinsic :: & use, intrinsic :: &
IEEE_arithmetic IEEE_arithmetic
use debug, only: & use debug, only: &
#ifdef DEBUG
debug_e, &
debug_i, &
debug_g, &
#endif
debug_level, & debug_level, &
debug_crystallite, & debug_crystallite, &
debug_levelBasic, & debug_levelBasic, &
debug_levelExtensive, & debug_levelExtensive, &
debug_levelSelective, & debug_levelSelective, &
debug_e, &
debug_i, &
debug_g, &
debug_StateLoopDistribution debug_StateLoopDistribution
use numerics, only: & use numerics, only: &
rTol_crystalliteState, & rTol_crystalliteState, &
@ -2036,14 +2040,16 @@ subroutine crystallite_integrateStateAdaptiveEuler()
use, intrinsic :: & use, intrinsic :: &
IEEE_arithmetic IEEE_arithmetic
use debug, only: & use debug, only: &
#ifdef DEBUG
debug_e, &
debug_i, &
debug_g, &
#endif
debug_level, & debug_level, &
debug_crystallite, & debug_crystallite, &
debug_levelBasic, & debug_levelBasic, &
debug_levelExtensive, & debug_levelExtensive, &
debug_levelSelective, & debug_levelSelective, &
debug_e, &
debug_i, &
debug_g, &
debug_StateLoopDistribution debug_StateLoopDistribution
use numerics, only: & use numerics, only: &
rTol_crystalliteState, & rTol_crystalliteState, &
@ -2396,14 +2402,16 @@ subroutine crystallite_integrateStateEuler()
use, intrinsic :: & use, intrinsic :: &
IEEE_arithmetic IEEE_arithmetic
use debug, only: & use debug, only: &
#ifdef DEBUG
debug_e, &
debug_i, &
debug_g, &
#endif
debug_level, & debug_level, &
debug_crystallite, & debug_crystallite, &
debug_levelBasic, & debug_levelBasic, &
debug_levelExtensive, & debug_levelExtensive, &
debug_levelSelective, & debug_levelSelective, &
debug_e, &
debug_i, &
debug_g, &
debug_StateLoopDistribution debug_StateLoopDistribution
use numerics, only: & use numerics, only: &
numerics_integrationMode, & numerics_integrationMode, &
@ -2619,9 +2627,11 @@ subroutine crystallite_integrateStateFPI()
use, intrinsic :: & use, intrinsic :: &
IEEE_arithmetic IEEE_arithmetic
use debug, only: & use debug, only: &
#ifdef DEBUG
debug_e, & debug_e, &
debug_i, & debug_i, &
debug_g, & debug_g, &
#endif
debug_level,& debug_level,&
debug_crystallite, & debug_crystallite, &
debug_levelBasic, & debug_levelBasic, &
@ -3068,14 +3078,16 @@ logical function crystallite_stateJump(ipc,ip,el)
IEEE_arithmetic IEEE_arithmetic
use prec, only: & use prec, only: &
dNeq0 dNeq0
#ifdef DEBUG
use debug, only: & use debug, only: &
debug_e, &
debug_i, &
debug_g, &
debug_level, & debug_level, &
debug_crystallite, & debug_crystallite, &
debug_levelExtensive, & debug_levelExtensive, &
debug_levelSelective, & debug_levelSelective
debug_e, & #endif
debug_i, &
debug_g
use material, only: & use material, only: &
plasticState, & plasticState, &
sourceState, & sourceState, &
@ -3205,9 +3217,11 @@ logical function crystallite_integrateStress(&
debug_levelBasic, & debug_levelBasic, &
debug_levelExtensive, & debug_levelExtensive, &
debug_levelSelective, & debug_levelSelective, &
#ifdef DEBUG
debug_e, & debug_e, &
debug_i, & debug_i, &
debug_g, & debug_g, &
#endif
debug_cumLpCalls, & debug_cumLpCalls, &
debug_cumLpTicks, & debug_cumLpTicks, &
debug_StressLoopLpDistribution, & debug_StressLoopLpDistribution, &
@ -3233,7 +3247,9 @@ logical function crystallite_integrateStress(&
math_Plain33to9, & math_Plain33to9, &
math_Plain9to33, & math_Plain9to33, &
math_Plain99to3333 math_Plain99to3333
#ifdef DEBUG
use mesh, only: mesh_element use mesh, only: mesh_element
#endif
implicit none implicit none
integer(pInt), intent(in):: el, & ! element index integer(pInt), intent(in):: el, & ! element index
@ -3297,8 +3313,8 @@ logical function crystallite_integrateStress(&
p, & p, &
jacoCounterLp, & jacoCounterLp, &
jacoCounterLi ! counters to check for Jacobian update jacoCounterLi ! counters to check for Jacobian update
integer(pLongInt) tick, & integer(pLongInt) :: tick = 0_pLongInt, &
tock, & tock = 0_pLongInt, &
tickrate, & tickrate, &
maxticks maxticks

40
src/homogenization.f90
View File

@ -16,7 +16,7 @@ module homogenization
! General variables for the homogenization at a material point ! General variables for the homogenization at a material point
implicit none implicit none
private private
real(pReal), dimension(:,:,:,:), allocatable, public :: & real(pReal), dimension(:,:,:,:), allocatable, public :: &
materialpoint_F0, & !< def grad of IP at start of FE increment 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_F, & !< def grad of IP to be reached at end of FE increment
materialpoint_P !< first P--K stress of IP materialpoint_P !< first P--K stress of IP
@ -128,7 +128,7 @@ subroutine homogenization_init
integer(pInt), dimension(:) , pointer :: thisNoutput integer(pInt), dimension(:) , pointer :: thisNoutput
character(len=64), dimension(:,:), pointer :: thisOutput character(len=64), dimension(:,:), pointer :: thisOutput
character(len=32) :: outputName !< name of output, intermediate fix until HDF5 output is ready 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 do p = 1,material_Nhomogenization
if (any(material_homog == p)) then if (any(material_homog == p)) then
i = homogenization_typeInstance(p) ! which instance of this homogenization type 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 select case(homogenization_type(p)) ! split per homogenization type
case (HOMOGENIZATION_NONE_ID) case (HOMOGENIZATION_NONE_ID)
outputName = HOMOGENIZATION_NONE_label outputName = HOMOGENIZATION_NONE_label
@ -217,10 +217,10 @@ subroutine homogenization_init
thisOutput => homogenization_RGC_output thisOutput => homogenization_RGC_output
thisSize => homogenization_RGC_sizePostResult thisSize => homogenization_RGC_sizePostResult
case default case default
knownHomogenization = .false. valid = .false.
end select end select
write(FILEUNIT,'(/,a,/)') '['//trim(homogenization_name(p))//']' write(FILEUNIT,'(/,a,/)') '['//trim(homogenization_name(p))//']'
if (knownHomogenization) then if (valid) then
write(FILEUNIT,'(a)') '(type)'//char(9)//trim(outputName) write(FILEUNIT,'(a)') '(type)'//char(9)//trim(outputName)
write(FILEUNIT,'(a,i4)') '(ngrains)'//char(9),homogenization_Ngrains(p) write(FILEUNIT,'(a,i4)') '(ngrains)'//char(9),homogenization_Ngrains(p)
if (homogenization_type(p) /= HOMOGENIZATION_NONE_ID) then if (homogenization_type(p) /= HOMOGENIZATION_NONE_ID) then
@ -230,7 +230,7 @@ subroutine homogenization_init
endif endif
endif endif
i = thermal_typeInstance(p) ! which instance of this thermal type 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 select case(thermal_type(p)) ! split per thermal type
case (THERMAL_isothermal_ID) case (THERMAL_isothermal_ID)
outputName = THERMAL_isothermal_label outputName = THERMAL_isothermal_label
@ -248,9 +248,9 @@ subroutine homogenization_init
thisOutput => thermal_conduction_output thisOutput => thermal_conduction_output
thisSize => thermal_conduction_sizePostResult thisSize => thermal_conduction_sizePostResult
case default case default
knownThermal = .false. valid = .false.
end select end select
if (knownThermal) then if (valid) then
write(FILEUNIT,'(a)') '(thermal)'//char(9)//trim(outputName) write(FILEUNIT,'(a)') '(thermal)'//char(9)//trim(outputName)
if (thermal_type(p) /= THERMAL_isothermal_ID) then if (thermal_type(p) /= THERMAL_isothermal_ID) then
do e = 1,thisNoutput(i) do e = 1,thisNoutput(i)
@ -259,7 +259,7 @@ subroutine homogenization_init
endif endif
endif endif
i = damage_typeInstance(p) ! which instance of this damage type 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 select case(damage_type(p)) ! split per damage type
case (DAMAGE_none_ID) case (DAMAGE_none_ID)
outputName = DAMAGE_none_label outputName = DAMAGE_none_label
@ -277,9 +277,9 @@ subroutine homogenization_init
thisOutput => damage_nonlocal_output thisOutput => damage_nonlocal_output
thisSize => damage_nonlocal_sizePostResult thisSize => damage_nonlocal_sizePostResult
case default case default
knownDamage = .false. valid = .false.
end select end select
if (knownDamage) then if (valid) then
write(FILEUNIT,'(a)') '(damage)'//char(9)//trim(outputName) write(FILEUNIT,'(a)') '(damage)'//char(9)//trim(outputName)
if (damage_type(p) /= DAMAGE_none_ID) then if (damage_type(p) /= DAMAGE_none_ID) then
do e = 1,thisNoutput(i) do e = 1,thisNoutput(i)
@ -288,7 +288,7 @@ subroutine homogenization_init
endif endif
endif endif
i = vacancyflux_typeInstance(p) ! which instance of this vacancy flux type 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 select case(vacancyflux_type(p)) ! split per vacancy flux type
case (VACANCYFLUX_isoconc_ID) case (VACANCYFLUX_isoconc_ID)
outputName = VACANCYFLUX_isoconc_label outputName = VACANCYFLUX_isoconc_label
@ -306,9 +306,9 @@ subroutine homogenization_init
thisOutput => vacancyflux_cahnhilliard_output thisOutput => vacancyflux_cahnhilliard_output
thisSize => vacancyflux_cahnhilliard_sizePostResult thisSize => vacancyflux_cahnhilliard_sizePostResult
case default case default
knownVacancyflux = .false. valid = .false.
end select end select
if (knownVacancyflux) then if (valid) then
write(FILEUNIT,'(a)') '(vacancyflux)'//char(9)//trim(outputName) write(FILEUNIT,'(a)') '(vacancyflux)'//char(9)//trim(outputName)
if (vacancyflux_type(p) /= VACANCYFLUX_isoconc_ID) then if (vacancyflux_type(p) /= VACANCYFLUX_isoconc_ID) then
do e = 1,thisNoutput(i) do e = 1,thisNoutput(i)
@ -317,7 +317,7 @@ subroutine homogenization_init
endif endif
endif endif
i = porosity_typeInstance(p) ! which instance of this porosity type 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 select case(porosity_type(p)) ! split per porosity type
case (POROSITY_none_ID) case (POROSITY_none_ID)
outputName = POROSITY_none_label outputName = POROSITY_none_label
@ -330,9 +330,9 @@ subroutine homogenization_init
thisOutput => porosity_phasefield_output thisOutput => porosity_phasefield_output
thisSize => porosity_phasefield_sizePostResult thisSize => porosity_phasefield_sizePostResult
case default case default
knownPorosity = .false. valid = .false.
end select end select
if (knownPorosity) then if (valid) then
write(FILEUNIT,'(a)') '(porosity)'//char(9)//trim(outputName) write(FILEUNIT,'(a)') '(porosity)'//char(9)//trim(outputName)
if (porosity_type(p) /= POROSITY_none_ID) then if (porosity_type(p) /= POROSITY_none_ID) then
do e = 1,thisNoutput(i) do e = 1,thisNoutput(i)
@ -341,7 +341,7 @@ subroutine homogenization_init
endif endif
endif endif
i = hydrogenflux_typeInstance(p) ! which instance of this hydrogen flux type 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 select case(hydrogenflux_type(p)) ! split per hydrogen flux type
case (HYDROGENFLUX_isoconc_ID) case (HYDROGENFLUX_isoconc_ID)
outputName = HYDROGENFLUX_isoconc_label outputName = HYDROGENFLUX_isoconc_label
@ -354,9 +354,9 @@ subroutine homogenization_init
thisOutput => hydrogenflux_cahnhilliard_output thisOutput => hydrogenflux_cahnhilliard_output
thisSize => hydrogenflux_cahnhilliard_sizePostResult thisSize => hydrogenflux_cahnhilliard_sizePostResult
case default case default
knownHydrogenflux = .false. valid = .false.
end select end select
if (knownHydrogenflux) then if (valid) then
write(FILEUNIT,'(a)') '(hydrogenflux)'//char(9)//trim(outputName) write(FILEUNIT,'(a)') '(hydrogenflux)'//char(9)//trim(outputName)
if (hydrogenflux_type(p) /= HYDROGENFLUX_isoconc_ID) then if (hydrogenflux_type(p) /= HYDROGENFLUX_isoconc_ID) then
do e = 1,thisNoutput(i) do e = 1,thisNoutput(i)

14
src/numerics.f90
View File

@ -120,9 +120,9 @@ module numerics
petsc_options = '' petsc_options = ''
integer(pInt), protected, public :: & 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 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 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 :: & 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 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 update_gamma = .false. !< update gamma operator with current stiffness, Default .false.: use initial stiffness
#endif #endif
@ -424,9 +424,9 @@ subroutine numerics_init
case ('err_stress_tolabs') case ('err_stress_tolabs')
err_stress_tolabs = IO_floatValue(line,chunkPos,2_pInt) err_stress_tolabs = IO_floatValue(line,chunkPos,2_pInt)
case ('continuecalculation') case ('continuecalculation')
continueCalculation = IO_intValue(line,chunkPos,2_pInt) continueCalculation = IO_intValue(line,chunkPos,2_pInt) > 0_pInt
case ('memory_efficient') 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') case ('fftw_timelimit')
fftw_timelimit = IO_floatValue(line,chunkPos,2_pInt) fftw_timelimit = IO_floatValue(line,chunkPos,2_pInt)
case ('fftw_plan_mode') case ('fftw_plan_mode')
@ -436,7 +436,7 @@ subroutine numerics_init
case ('divergence_correction') case ('divergence_correction')
divergence_correction = IO_intValue(line,chunkPos,2_pInt) divergence_correction = IO_intValue(line,chunkPos,2_pInt)
case ('update_gamma') 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') case ('petsc_options')
petsc_options = trim(line(chunkPos(4):)) petsc_options = trim(line(chunkPos(4):))
case ('spectralsolver','myspectralsolver') case ('spectralsolver','myspectralsolver')
@ -599,7 +599,7 @@ subroutine numerics_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! spectral parameters ! spectral parameters
#ifdef Spectral #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,L8)') ' memory_efficient: ',memory_efficient
write(6,'(a24,1x,i8)') ' divergence_correction: ',divergence_correction write(6,'(a24,1x,i8)') ' divergence_correction: ',divergence_correction
write(6,'(a24,1x,a)') ' spectral_derivative: ',trim(spectral_derivative) 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_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') if (err_hydrogenflux_tolrel <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_hydrogenflux_tolrel')
#ifdef Spectral #ifdef Spectral
if (continueCalculation /= 0_pInt .and. &
continueCalculation /= 1_pInt) call IO_error(301_pInt,ext_msg='continueCalculation')
if (divergence_correction < 0_pInt .or. & if (divergence_correction < 0_pInt .or. &
divergence_correction > 2_pInt) call IO_error(301_pInt,ext_msg='divergence_correction') divergence_correction > 2_pInt) call IO_error(301_pInt,ext_msg='divergence_correction')
if (update_gamma .and. & if (update_gamma .and. &
@ -713,7 +711,7 @@ subroutine numerics_init
if (polarAlpha <= 0.0_pReal .or. & if (polarAlpha <= 0.0_pReal .or. &
polarAlpha > 2.0_pReal) call IO_error(301_pInt,ext_msg='polarAlpha') polarAlpha > 2.0_pReal) call IO_error(301_pInt,ext_msg='polarAlpha')
if (polarBeta < 0.0_pReal .or. & 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 #endif
end subroutine numerics_init end subroutine numerics_init

496
src/spectral_mech_AL.f90
View File

@ -43,17 +43,20 @@ module spectral_mech_AL
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc. ! stress, stiffness and compliance average etc.
real(pReal), private, dimension(3,3) :: & 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 = math_I3, & !< current prescribed deformation gradient
F_aim_lastInc = math_I3, & !< previous average deformation gradient F_aim_lastInc = math_I3, & !< previous average deformation gradient
F_av = 0.0_pReal, & !< average incompatible def grad field F_av = 0.0_pReal, & !< average incompatible def grad field
P_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress 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 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) :: & real(pReal), private, dimension(3,3,3,3) :: &
C_volAvg = 0.0_pReal, & !< current volume average stiffness C_volAvg = 0.0_pReal, & !< current volume average stiffness
C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness
C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 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) S = 0.0_pReal, & !< current compliance (filled up with zeros)
C_scale = 0.0_pReal, & C_scale = 0.0_pReal, &
S_scale = 0.0_pReal S_scale = 0.0_pReal
@ -62,7 +65,7 @@ module spectral_mech_AL
err_BC, & !< deviation from stress BC err_BC, & !< deviation from stress BC
err_curl, & !< RMS of curl of F err_curl, & !< RMS of curl of F
err_div !< RMS of div of P err_div !< RMS of div of P
logical, private :: ForwardData
integer(pInt), private :: & integer(pInt), private :: &
totalIter = 0_pInt !< total iteration in current increment totalIter = 0_pInt !< total iteration in current increment
@ -80,7 +83,7 @@ module spectral_mech_AL
contains 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) !> @todo use sourced allocation, e.g. allocate(Fdot,source = F_lastInc)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine AL_init subroutine AL_init
@ -102,12 +105,15 @@ subroutine AL_init
use numerics, only: & use numerics, only: &
worldrank, & worldrank, &
worldsize worldsize
use homogenization, only: &
materialpoint_F0
use DAMASK_interface, only: & use DAMASK_interface, only: &
getSolverJobName getSolverJobName
use spectral_utilities, only: & use spectral_utilities, only: &
Utilities_constitutiveResponse, & Utilities_constitutiveResponse, &
Utilities_updateGamma, & Utilities_updateGamma, &
Utilities_updateIPcoords Utilities_updateIPcoords, &
wgt
use mesh, only: & use mesh, only: &
grid, & grid, &
grid3 grid3
@ -120,7 +126,11 @@ subroutine AL_init
temp33_Real = 0.0_pReal temp33_Real = 0.0_pReal
PetscErrorCode :: ierr 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), dimension(:), allocatable :: localK
integer(pInt) :: proc integer(pInt) :: proc
character(len=1024) :: rankStr character(len=1024) :: rankStr
@ -160,51 +170,43 @@ subroutine AL_init
DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point
grid(1),grid(2),grid(3), & ! global grid grid(1),grid(2),grid(3), & ! global grid
1 , 1, worldsize, & 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 grid(1),grid(2),localK, & ! local grid
da,ierr) ! handle, error da,ierr) ! handle, error
CHKERRQ(ierr) 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) 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) 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) CHKERRQ(ierr)
call SNESSetConvergenceTest(snes,AL_converged,PETSC_NULL_OBJECT,PETSC_NULL_FUNCTION,ierr) call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional cli arguments
CHKERRQ(ierr)
call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init fields ! init fields
call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! places pointer xx_psc on PETSc data call DMDAVecGetArrayF90(da,solution_vec,FandF_lambda,ierr); CHKERRQ(ierr) ! places pointer on PETSc data
F => xx_psc(0:8,:,:,:) F => FandF_lambda( 0: 8,:,:,:)
F_lambda => xx_psc(9:17,:,:,:) F_lambda => FandF_lambda( 9:17,:,:,:)
restart: if (restartInc > 1_pInt) then 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)') & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'reading values of increment ', restartInc - 1_pInt, ' from file' 'reading values of increment ', restartInc-1_pInt, ' from file'
flush(6) flush(6)
endif
write(rankStr,'(a1,i0)')'_',worldrank write(rankStr,'(a1,i0)')'_',worldrank
call IO_read_realFile(777,'F'//trim(rankStr),trim(getSolverJobName()),size(F)) call IO_read_realFile(777,'F'//trim(rankStr),trim(getSolverJobName()),size(F))
read (777,rec=1) F read (777,rec=1) F; close (777)
close (777)
call IO_read_realFile(777,'F_lastInc'//trim(rankStr),trim(getSolverJobName()),size(F_lastInc)) call IO_read_realFile(777,'F_lastInc'//trim(rankStr),trim(getSolverJobName()),size(F_lastInc))
read (777,rec=1) F_lastInc read (777,rec=1) F_lastInc; close (777)
close (777)
call IO_read_realFile(777,'F_lambda'//trim(rankStr),trim(getSolverJobName()),size(F_lambda)) call IO_read_realFile(777,'F_lambda'//trim(rankStr),trim(getSolverJobName()),size(F_lambda))
read (777,rec=1) F_lambda read (777,rec=1) F_lambda; close (777)
close (777) call IO_read_realFile(777,'F_lambda_lastInc'//trim(rankStr),trim(getSolverJobName()),size(F_lambda_lastInc))
call IO_read_realFile(777,'F_lambda_lastInc'//trim(rankStr),& read (777,rec=1) F_lambda_lastInc; close (777)
trim(getSolverJobName()),size(F_lambda_lastInc)) call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(F_aimDot))
read (777,rec=1) F_lambda_lastInc read (777,rec=1) F_aimDot; close (777)
close (777) F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F
call IO_read_realFile(777,'F_aim', trim(getSolverJobName()),size(F_aim)) F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc
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)
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_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 = reshape(F_lastInc,[9,grid(1),grid(2),grid3])
@ -212,30 +214,30 @@ subroutine AL_init
F_lambda_lastInc = F_lastInc F_lambda_lastInc = F_lastInc
endif restart 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_updateIPcoords(reshape(F,shape(F_lastInc)))
call Utilities_constitutiveResponse(F_lastInc, reshape(F,shape(F_lastInc)), & call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
0.0_pReal,P,C_volAvg,C_minMaxAvg,temp33_Real,.false.,math_I3) reshape(F,shape(F_lastInc)), & ! target F
0.0_pReal, & ! time increment
math_I3) ! no rotation of boundary condition
nullify(F) nullify(F)
nullify(F_lambda) 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 restartRead: 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 .and. worldrank == 0_pInt) &
write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'reading more values of increment', restartInc - 1_pInt, 'from file' 'reading more values of increment ', restartInc-1_pInt, ' from file'
flush(6) flush(6)
call IO_read_realFile(777,'C_volAvg',trim(getSolverJobName()),size(C_volAvg)) call IO_read_realFile(777,'C_volAvg',trim(getSolverJobName()),size(C_volAvg))
read (777,rec=1) C_volAvg read (777,rec=1) C_volAvg; close (777)
close (777)
call IO_read_realFile(777,'C_volAvgLastInc',trim(getSolverJobName()),size(C_volAvgLastInc)) call IO_read_realFile(777,'C_volAvgLastInc',trim(getSolverJobName()),size(C_volAvgLastInc))
read (777,rec=1) C_volAvgLastInc read (777,rec=1) C_volAvgLastInc; close (777)
close (777)
call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg)) call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg))
read (777,rec=1) C_minMaxAvg read (777,rec=1) C_minMaxAvg; close (777)
close (777)
endif restartRead endif restartRead
call Utilities_updateGamma(C_minMaxAvg,.True.) call Utilities_updateGamma(C_minMaxAvg,.true.)
C_scale = C_minMaxAvg C_scale = C_minMaxAvg
S_scale = math_invSym3333(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 !> @brief solution for the AL scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
type(tSolutionState) function & type(tSolutionState) function AL_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC)
AL_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC)
use IO, only: & use IO, only: &
IO_error IO_error
use numerics, only: & use numerics, only: &
@ -265,13 +266,13 @@ type(tSolutionState) function &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! input data for solution ! 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) :: & character(len=*), intent(in) :: &
incInfoIn 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 real(pReal), dimension(3,3), intent(in) :: rotation_BC
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -304,18 +305,17 @@ type(tSolutionState) function &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! solve BVP ! solve BVP
call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr) call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr); CHKERRQ(ierr)
CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! check convergence ! check convergence
call SNESGetConvergedReason(snes,reason,ierr) call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr)
CHKERRQ(ierr)
AL_solution%converged = reason > 0
AL_solution%iterationsNeeded = totalIter
AL_solution%termIll = terminallyIll AL_solution%termIll = terminallyIll
terminallyIll = .false. terminallyIll = .false.
if (reason == -4) call IO_error(893_pInt) if (reason == -4) call IO_error(893_pInt) ! MPI error
if (reason < 1) AL_solution%converged = .false.
AL_solution%iterationsNeeded = totalIter
end function AL_solution end function AL_solution
@ -330,8 +330,8 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
polarAlpha, & polarAlpha, &
polarBeta polarBeta
use mesh, only: & use mesh, only: &
grid3, & grid, &
grid grid3
use IO, only: & use IO, only: &
IO_intOut IO_intOut
use math, only: & use math, only: &
@ -340,6 +340,10 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
math_mul3333xx33, & math_mul3333xx33, &
math_invSym3333, & math_invSym3333, &
math_mul33x33 math_mul33x33
use debug, only: &
debug_level, &
debug_spectral, &
debug_spectralRotation
use spectral_utilities, only: & use spectral_utilities, only: &
wgt, & wgt, &
tensorField_real, & tensorField_real, &
@ -349,27 +353,17 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
Utilities_constitutiveResponse, & Utilities_constitutiveResponse, &
Utilities_divergenceRMS, & Utilities_divergenceRMS, &
Utilities_curlRMS Utilities_curlRMS
use debug, only: &
debug_level, &
debug_spectral, &
debug_spectralRotation
use homogenization, only: & use homogenization, only: &
materialpoint_dPdF materialpoint_dPdF
use FEsolving, only: & use FEsolving, only: &
terminallyIll terminallyIll
implicit none implicit none
!-------------------------------------------------------------------------------------------------- DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in
! strange syntax in the next line because otherwise macros expand beyond 132 character limit PetscScalar, &
DMDALocalInfo, dimension(& target, dimension(3,3,2, XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: x_scal !< what is this?
DMDA_LOCAL_INFO_SIZE) :: & PetscScalar, &
in target, dimension(3,3,2, X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: f_scal !< what is this?
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(:,:,:,:,:) :: & PetscScalar, pointer, dimension(:,:,:,:,:) :: &
F, & F, &
F_lambda, & F_lambda, &
@ -387,33 +381,33 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
SNESGetNumberFunctionEvals, & SNESGetNumberFunctionEvals, &
SNESGetIterationNumber SNESGetIterationNumber
F => x_scal(1:3,1:3,1,& F => x_scal(1:3,1:3,1,&
XG_RANGE,YG_RANGE,ZG_RANGE) XG_RANGE,YG_RANGE,ZG_RANGE)
F_lambda => x_scal(1:3,1:3,2,& F_lambda => x_scal(1:3,1:3,2,&
XG_RANGE,YG_RANGE,ZG_RANGE) XG_RANGE,YG_RANGE,ZG_RANGE)
residual_F => f_scal(1:3,1:3,1,& residual_F => f_scal(1:3,1:3,1,&
X_RANGE,Y_RANGE,Z_RANGE) X_RANGE, Y_RANGE, Z_RANGE)
residual_F_lambda => f_scal(1:3,1:3,2,& residual_F_lambda => f_scal(1:3,1:3,2,&
X_RANGE,Y_RANGE,Z_RANGE) X_RANGE, Y_RANGE, Z_RANGE)
call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr)
call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr)
F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt 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) 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 call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr)
newIteration: if(totalIter <= PETScIter) then 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 totalIter = totalIter + 1_pInt
write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') &
' @ Iteration ', itmin, '≤',totalIter, '≤', itmax trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) ' 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 =', & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
math_transpose33(F_aim) ' deformation gradient aim =', math_transpose33(F_aim)
flush(6) flush(6)
endif newIteration 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) -& 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), & 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)) math_mul3333xx33(C_scale,F_lambda(1:3,1:3,i,j,k) - math_I3))
enddo; enddo; enddo enddo; enddo; enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -435,24 +428,23 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
call utilities_FFTtensorBackward() call utilities_FFTtensorBackward()
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! constructing residual ! constructing F_lambda residual
residual_F_lambda = polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) 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 ! evaluate constitutive response
P_avLastEval = P_av 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) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
ForwardData = .False.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate divergence ! calculate divergence
tensorField_real = 0.0_pReal 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() call utilities_FFTtensorForward()
err_div = Utilities_divergenceRMS() err_div = Utilities_divergenceRMS() !< root mean squared error in divergence of stress
call utilities_FFTtensorBackward()
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! constructing residual ! 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) - & residual_F(1:3,1:3,i,j,k) - &
math_mul33x33(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))) & 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 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 tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F
call utilities_FFTtensorForward() call utilities_FFTtensorForward()
err_curl = Utilities_curlRMS() err_curl = Utilities_curlRMS()
call utilities_FFTtensorBackward()
nullify(F)
nullify(F_lambda)
nullify(residual_F)
nullify(residual_F_lambda)
end subroutine AL_formResidual end subroutine AL_formResidual
@ -488,8 +483,8 @@ subroutine AL_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr
err_div_tolAbs, & err_div_tolAbs, &
err_curl_tolRel, & err_curl_tolRel, &
err_curl_tolAbs, & err_curl_tolAbs, &
err_stress_tolAbs, & err_stress_tolRel, &
err_stress_tolRel err_stress_tolAbs
use math, only: & use math, only: &
math_mul3333xx33 math_mul3333xx33
use FEsolving, only: & use FEsolving, only: &
@ -508,24 +503,24 @@ subroutine AL_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr
real(pReal) :: & real(pReal) :: &
curlTol, & curlTol, &
divTol, & divTol, &
BC_tol BCTol
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress BC handling ! stress BC handling
F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc 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) + & 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 mask_stress * (P_av-params%stress_BC))) ! mask = 0.0 for no bc
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! error calculation ! error calculation
curlTol = max(maxval(abs(F_aim-math_I3))*err_curl_tolRel,err_curl_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) 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) BCTol = max(maxval(abs(P_av)) *err_stress_tolRel,err_stress_tolAbs)
converged: if ((totalIter >= itmin .and. & converged: if ((totalIter >= itmin .and. &
all([ err_div/divTol, & all([ err_div /divTol, &
err_curl/curlTol, & err_curl/curlTol, &
err_BC/BC_tol ] < 1.0_pReal)) & err_BC /BCTol ] < 1.0_pReal)) &
.or. terminallyIll) then .or. terminallyIll) then
reason = 1 reason = 1
elseif (totalIter >= itmax) then converged elseif (totalIter >= itmax) then converged
@ -537,12 +532,12 @@ subroutine AL_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report ! report
write(6,'(1/,a)') ' ... reporting .............................................................' write(6,'(1/,a)') ' ... reporting .............................................................'
write(6,'(/,a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', & write(6,'(/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', &
err_curl/curlTol,' (',err_curl,' -, tol =',curlTol,')' err_div/divTol, ' (',err_div, ' / m, tol = ',divTol,')'
write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', &
err_div/divTol, ' (',err_div, ' / m, tol =',divTol,')' err_curl/curlTol,' (',err_curl,' -, tol = ',curlTol,')'
write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', & write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', &
err_BC/BC_tol, ' (',err_BC, ' Pa, tol =',BC_tol,')' err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')'
write(6,'(/,a)') ' ===========================================================================' write(6,'(/,a)') ' ==========================================================================='
flush(6) flush(6)
@ -550,130 +545,142 @@ end subroutine AL_converged
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief forwarding routine !> @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) subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC)
use math, only: & use math, only: &
math_mul33x33, & math_mul33x33, &
math_mul3333xx33, & math_mul3333xx33, &
math_transpose33, & math_transpose33, &
math_rotate_backward33 math_rotate_backward33
use numerics, only: & use numerics, only: &
worldrank worldrank
use mesh, only: & use homogenization, only: &
grid3, & materialpoint_F0
grid use mesh, only: &
use spectral_utilities, only: & grid, &
Utilities_calculateRate, & grid3
Utilities_forwardField, & use CPFEM2, only: &
Utilities_updateIPcoords, & CPFEM_age
tBoundaryCondition, & use spectral_utilities, only: &
cutBack Utilities_calculateRate, &
use IO, only: & Utilities_forwardField, &
IO_write_JobRealFile Utilities_updateIPcoords, &
use FEsolving, only: & tBoundaryCondition, &
restartWrite cutBack
use IO, only: &
IO_write_JobRealFile
use FEsolving, only: &
restartWrite
implicit none implicit none
real(pReal), intent(in) :: & logical, intent(in) :: &
timeinc_old, & guess
timeinc, & real(pReal), intent(in) :: &
loadCaseTime !< remaining time of current load case timeinc_old, &
type(tBoundaryCondition), intent(in) :: & timeinc, &
stress_BC, & loadCaseTime !< remaining time of current load case
deformation_BC type(tBoundaryCondition), intent(in) :: &
real(pReal), dimension(3,3), intent(in) :: rotation_BC stress_BC, &
logical, intent(in) :: & deformation_BC
guess real(pReal), dimension(3,3), intent(in) ::&
PetscErrorCode :: ierr rotation_BC
PetscScalar, dimension(:,:,:,:), pointer :: xx_psc, F, F_lambda PetscErrorCode :: ierr
integer(pInt) :: i, j, k PetscScalar, dimension(:,:,:,:), pointer :: FandF_lambda, F, F_lambda
real(pReal), dimension(3,3) :: F_lambda33 integer(pInt) :: i, j, k
character(len=1024) :: rankStr real(pReal), dimension(3,3) :: F_lambda33
character(len=32) :: rankStr
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update coordinates and rate and forward last inc ! update coordinates and rate and forward last inc
call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr) call DMDAVecGetArrayF90(da,solution_vec,FandF_lambda,ierr); CHKERRQ(ierr)
F => xx_psc(0:8,:,:,:) F => FandF_lambda( 0: 8,:,:,:)
F_lambda => xx_psc(9:17,:,:,:) F_lambda => FandF_lambda( 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 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 if (worldrank == 0_pInt) then
F_aim = F_aim_lastInc call IO_write_jobRealFile(777,'C_volAvg',size(C_volAvg))
F_lambda = reshape(F_lambda_lastInc,[9,grid(1),grid(2),grid3]) write (777,rec=1) C_volAvg; close(777)
F = reshape(F_lastInc, [9,grid(1),grid(2),grid3]) call IO_write_jobRealFile(777,'C_volAvgLastInc',size(C_volAvgLastInc))
C_volAvg = C_volAvgLastInc write (777,rec=1) C_volAvgLastInc; close(777)
else ! call IO_write_jobRealFile(777,'C_minMaxAvg',size(C_volAvg))
ForwardData = .True. ! write (777,rec=1) C_minMaxAvg; close(777)
C_volAvgLastInc = C_volAvg ! call IO_write_jobRealFile(777,'C_minMaxAvgLastInc',size(C_volAvgLastInc))
!-------------------------------------------------------------------------------------------------- ! write (777,rec=1) C_minMaxAvgLastInc; close(777)
! calculate rate for aim call IO_write_jobRealFile(777,'F_aimDot',size(F_aimDot))
if (deformation_BC%myType=='l') then ! calculate f_aimDot from given L and current F write (777,rec=1) F_aimDot; close(777)
f_aimDot = deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim) endif
elseif(deformation_BC%myType=='fdot') then ! f_aimDot is prescribed
f_aimDot = deformation_BC%maskFloat * deformation_BC%values write(rankStr,'(a1,i0)')'_',worldrank
elseif(deformation_BC%myType=='f') then ! aim at end of load case is prescribed call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file
f_aimDot = deformation_BC%maskFloat * (deformation_BC%values -F_aim)/loadCaseTime write (777,rec=1) F; close (777)
endif call IO_write_jobRealFile(777,'F_lastInc'//trim(rankStr),size(F_lastInc)) ! writing F_lastInc field to file
if (guess) f_aimDot = f_aimDot + stress_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old write (777,rec=1) F_lastInc; close (777)
F_aim_lastInc = F_aim 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 ! update average and local deformation gradients
call utilities_updateIPcoords(F) F_aim = F_aim_lastInc + F_aimDot * timeinc
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
F_aim = F_aim + 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])
! update local deformation gradient if (guess) then
F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! ensure that it matches rotated F_aim F_lambda = reshape(Utilities_forwardField(timeinc,F_lambda_lastInc,F_lambdadot), &
math_rotate_backward33(F_aim,rotation_BC)), & [9,grid(1),grid(2),grid3]) ! does not have any average value as boundary condition
[9,grid(1),grid(2),grid3]) else
F_lambda = reshape(Utilities_forwardField(timeinc,F_lambda_lastInc,F_lambdadot), & do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1)
[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_lambda33 = reshape(F_lambda(1:9,i,j,k),[3,3]) F_lambda33 = reshape(F_lambda(1:9,i,j,k),[3,3])
F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, & F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, &
math_mul3333xx33(C_scale,& 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)& F_lambda33) -math_I3))*0.5_pReal)&
+ math_I3 + math_I3
F_lambda(1:9,i,j,k) = reshape(F_lambda33,[9]) F_lambda(1:9,i,j,k) = reshape(F_lambda33,[9])
enddo; enddo; enddo enddo; enddo; enddo
endif endif
call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr)
nullify(F)
nullify(F_lambda)
call DMDAVecRestoreArrayF90(da,solution_vec,FandF_lambda,ierr); CHKERRQ(ierr)
end subroutine AL_forward end subroutine AL_forward

378
src/spectral_mech_Basic.f90
View File

@ -38,21 +38,27 @@ module spectral_mech_basic
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc. ! stress, stiffness and compliance average etc.
real(pReal), private, dimension(3,3) :: & real(pReal), private, dimension(3,3) :: &
F_aim = math_I3, & F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient
F_aim_lastIter = math_I3, & F_aim = math_I3, & !< current prescribed deformation gradient
F_aim_lastInc = math_I3, & F_aim_lastInc = math_I3, & !< previous average deformation gradient
P_av = 0.0_pReal, & P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress
F_aimDot=0.0_pReal
character(len=1024), private :: incInfo character(len=1024), private :: incInfo !< time and increment information
real(pReal), private, dimension(3,3,3,3) :: & real(pReal), private, dimension(3,3,3,3) :: &
C_volAvg = 0.0_pReal, & !< current volume average stiffness C_volAvg = 0.0_pReal, & !< current volume average stiffness
C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness
C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 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
real(pReal), private :: err_stress, err_div S = 0.0_pReal !< current compliance (filled up with zeros)
logical, private :: ForwardData
real(pReal), private :: &
err_BC, & !< deviation from stress BC
err_div !< RMS of div of P
integer(pInt), private :: & integer(pInt), private :: &
totalIter = 0_pInt !< total iteration in current increment totalIter = 0_pInt !< total iteration in current increment
real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal
public :: & public :: &
@ -69,7 +75,7 @@ module spectral_mech_basic
contains 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 subroutine basicPETSc_init
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 #if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
@ -90,6 +96,8 @@ subroutine basicPETSc_init
use numerics, only: & use numerics, only: &
worldrank, & worldrank, &
worldsize worldsize
use homogenization, only: &
materialpoint_F0
use DAMASK_interface, only: & use DAMASK_interface, only: &
getSolverJobName getSolverJobName
use spectral_utilities, only: & use spectral_utilities, only: &
@ -132,8 +140,8 @@ subroutine basicPETSc_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate global fields ! allocate global fields
allocate (F_lastInc (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) allocate (Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc ! initialize solver specific parts of PETSc
@ -152,11 +160,10 @@ subroutine basicPETSc_init
grid(1),grid(2),localK, & ! local grid grid(1),grid(2),localK, & ! local grid
da,ierr) ! handle, error da,ierr) ! handle, error
CHKERRQ(ierr) 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 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 call DMDASNESSetFunctionLocal(da,INSERT_VALUES,BasicPETSC_formResidual,PETSC_NULL_OBJECT,ierr) ! residual vector of same shape as solution vector
CHKERRQ(ierr) 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" call SNESSetConvergenceTest(snes,BasicPETSC_converged,PETSC_NULL_OBJECT,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged"
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional cli arguments 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 call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with
restart: if (restartInc > 1_pInt) then 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)') & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'reading values of increment ', restartInc - 1_pInt, ' from file' 'reading values of increment ', restartInc-1_pInt, ' from file'
flush(6) flush(6)
endif
write(rankStr,'(a1,i0)')'_',worldrank write(rankStr,'(a1,i0)')'_',worldrank
call IO_read_realFile(777,'F'//trim(rankStr),trim(getSolverJobName()),size(F)) call IO_read_realFile(777,'F'//trim(rankStr),trim(getSolverJobName()),size(F))
read (777,rec=1) F read (777,rec=1) F; close (777)
close (777)
call IO_read_realFile(777,'F_lastInc'//trim(rankStr),trim(getSolverJobName()),size(F_lastInc)) call IO_read_realFile(777,'F_lastInc'//trim(rankStr),trim(getSolverJobName()),size(F_lastInc))
read (777,rec=1) F_lastInc read (777,rec=1) F_lastInc; close (777)
close (777) call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(F_aimDot))
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 = 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 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_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 = reshape(F_lastInc,[9,grid(1),grid(2),grid3])
endif restart 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_updateIPcoords(reshape(F,shape(F_lastInc)))
call Utilities_constitutiveResponse(F_lastInc, reshape(F,shape(F_lastInc)), & call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
0.0_pReal, & reshape(F,shape(F_lastInc)), & ! target F
P, & 0.0_pReal, & ! time increment
C_volAvg,C_minMaxAvg, & ! global average of stiffness and (min+max)/2 math_I3) ! no rotation of boundary condition
temp33_Real, &
.false., &
math_I3)
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! write data back to PETSc call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! write data back to PETSc
! QUESTION: why not writing back right after reading (l.189)?
restartRead: if (restartInc > 1_pInt) then restartRead: if (restartInc > 1_pInt) then ! QUESTION: are those values not calc'ed by constitutiveResponse? why reading from file?
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) & if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) &
write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'reading more values of increment', restartInc - 1_pInt, 'from file' 'reading more values of increment ', restartInc-1_pInt, ' from file'
flush(6) flush(6)
call IO_read_realFile(777,'C_volAvg',trim(getSolverJobName()),size(C_volAvg)) call IO_read_realFile(777,'C_volAvg',trim(getSolverJobName()),size(C_volAvg))
read (777,rec=1) C_volAvg read (777,rec=1) C_volAvg; close (777)
close (777)
call IO_read_realFile(777,'C_volAvgLastInc',trim(getSolverJobName()),size(C_volAvgLastInc)) call IO_read_realFile(777,'C_volAvgLastInc',trim(getSolverJobName()),size(C_volAvgLastInc))
read (777,rec=1) C_volAvgLastInc read (777,rec=1) C_volAvgLastInc; close (777)
close (777)
call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg)) call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg))
read (777,rec=1) C_minMaxAvg read (777,rec=1) C_minMaxAvg; close (777)
close (777)
endif restartRead endif restartRead
call Utilities_updateGamma(C_minmaxAvg,.True.) call Utilities_updateGamma(C_minMaxAvg,.true.)
end subroutine basicPETSc_init end subroutine basicPETSc_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief solution for the Basic PETSC scheme with internal iterations !> @brief solution for the Basic PETSC scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
type(tSolutionState) function & type(tSolutionState) function basicPETSc_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC)
basicPETSc_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC)
use IO, only: & use IO, only: &
IO_error IO_error
use numerics, only: & use numerics, only: &
@ -238,13 +238,13 @@ type(tSolutionState) function &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! input data for solution ! 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) :: & character(len=*), intent(in) :: &
incInfoIn 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 real(pReal), dimension(3,3), intent(in) :: rotation_BC
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -261,7 +261,7 @@ type(tSolutionState) function &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator) ! update stiffness (and gamma operator)
S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) 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 ! solve BVP
call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr) call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr); CHKERRQ(ierr)
CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! check convergence ! check convergence
call SNESGetConvergedReason(snes,reason,ierr) call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr)
CHKERRQ(ierr)
BasicPETSc_solution%converged = reason > 0
basicPETSC_solution%iterationsNeeded = totalIter
basicPETSc_solution%termIll = terminallyIll basicPETSc_solution%termIll = terminallyIll
terminallyIll = .false. terminallyIll = .false.
BasicPETSc_solution%converged =.true. if (reason == -4) call IO_error(893_pInt) ! MPI error
if (reason == -4) call IO_error(893_pInt)
if (reason < 1) basicPETSC_solution%converged = .false.
basicPETSC_solution%iterationsNeeded = totalIter
end function BasicPETSc_solution end function BasicPETSc_solution
@ -322,19 +320,18 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
terminallyIll terminallyIll
implicit none implicit none
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in
in PetscScalar, &
PetscScalar, dimension(3,3, & dimension(3,3, XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: x_scal !< what is this?
XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: & PetscScalar, &
x_scal dimension(3,3, X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: f_scal !< what is this?
PetscScalar, dimension(3,3, &
X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: &
f_scal
PetscInt :: & PetscInt :: &
PETScIter, & PETScIter, &
nfuncs nfuncs
PetscObject :: dummy PetscObject :: dummy
PetscErrorCode :: ierr PetscErrorCode :: ierr
real(pReal), dimension(3,3) :: &
deltaF_aim
external :: & external :: &
SNESGetNumberFunctionEvals, & SNESGetNumberFunctionEvals, &
@ -343,46 +340,45 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr)
call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr)
if(nfuncs== 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment
newIteration: if(totalIter <= PETScIter) then
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report begin of new iteration ! begin of new iteration
newIteration: if (totalIter <= PETScIter) then
totalIter = totalIter + 1_pInt totalIter = totalIter + 1_pInt
write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') &
' @ Iteration ', itmin, '≤',totalIter, '≤', itmax trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) ' 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 =', & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
math_transpose33(F_aim) ' deformation gradient aim =', math_transpose33(F_aim)
flush(6) flush(6)
endif newIteration endif newIteration
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate constitutive response ! evaluate constitutive response
call Utilities_constitutiveResponse(F_lastInc,x_scal,params%timeinc, & call Utilities_constitutiveResponse(f_scal,P_av,C_volAvg,C_minMaxAvg, &
f_scal,C_volAvg,C_minmaxAvg,P_av,ForwardData,params%rotation_BC) x_scal,params%timeinc, params%rotation_BC)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
ForwardData = .false.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress BC handling ! stress BC handling
F_aim_lastIter = F_aim deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC)
F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc F_aim = F_aim - deltaF_aim
err_stress = maxval(abs(mask_stress * (P_av - params%stress_BC))) ! mask = 0.0 for no bc 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 ! updated deformation gradient using fix point algorithm of basic scheme
tensorField_real = 0.0_pReal tensorField_real = 0.0_pReal
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) = f_scal
call utilities_FFTtensorForward() call utilities_FFTtensorForward() ! FFT forward of global "tensorField_real"
err_div = Utilities_divergenceRMS() err_div = Utilities_divergenceRMS() ! divRMS of tensorField_fourier
call utilities_fourierGammaConvolution(math_rotate_backward33(F_aim_lastIter-F_aim,params%rotation_BC)) call utilities_fourierGammaConvolution(math_rotate_backward33(deltaF_aim,params%rotation_BC)) ! convolution of Gamma and tensorField_fourier, with arg
call utilities_FFTtensorBackward() call utilities_FFTtensorBackward() ! FFT backward of global tensorField_fourier
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! constructing residual ! 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 end subroutine BasicPETSc_formResidual
@ -413,14 +409,14 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du
PetscErrorCode :: ierr PetscErrorCode :: ierr
real(pReal) :: & real(pReal) :: &
divTol, & divTol, &
stressTol BCTol
divTol = max(maxval(abs(P_av))*err_div_tolRel,err_div_tolAbs) divTol = max(maxval(abs(P_av))*err_div_tolRel ,err_div_tolAbs)
stressTol = max(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs) BCTol = max(maxval(abs(P_av))*err_stress_tolRel,err_stress_tolAbs)
converged: if ((totalIter >= itmin .and. & converged: if ((totalIter >= itmin .and. &
all([ err_div/divTol, & all([ err_div/divTol, &
err_stress/stressTol ] < 1.0_pReal)) & err_BC /BCTol ] < 1.0_pReal)) &
.or. terminallyIll) then .or. terminallyIll) then
reason = 1 reason = 1
elseif (totalIter >= itmax) then converged elseif (totalIter >= itmax) then converged
@ -433,9 +429,9 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du
! report ! report
write(6,'(1/,a)') ' ... reporting .............................................................' write(6,'(1/,a)') ' ... reporting .............................................................'
write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', &
err_div/divTol, ' (',err_div,' / m, tol =',divTol,')' err_div/divTol, ' (',err_div,' / m, tol = ',divTol,')'
write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', & write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', &
err_stress/stressTol, ' (',err_stress, ' Pa, tol =',stressTol,')' err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')'
write(6,'(/,a)') ' ===========================================================================' write(6,'(/,a)') ' ==========================================================================='
flush(6) flush(6)
@ -443,106 +439,118 @@ end subroutine BasicPETSc_converged
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief forwarding routine !> @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) subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC)
use math, only: & use math, only: &
math_mul33x33 ,& math_mul33x33 ,&
math_rotate_backward33 math_rotate_backward33
use numerics, only: & use numerics, only: &
worldrank worldrank
use mesh, only: & use homogenization, only: &
grid, & materialpoint_F0
grid3 use mesh, only: &
use spectral_utilities, only: & grid, &
Utilities_calculateRate, & grid3
Utilities_forwardField, & use CPFEM2, only: &
Utilities_updateIPcoords, & CPFEM_age
tBoundaryCondition, & use spectral_utilities, only: &
cutBack Utilities_calculateRate, &
use IO, only: & Utilities_forwardField, &
IO_write_JobRealFile Utilities_updateIPcoords, &
use FEsolving, only: & tBoundaryCondition, &
restartWrite cutBack
use IO, only: &
IO_write_JobRealFile
use FEsolving, only: &
restartWrite
implicit none implicit none
real(pReal), intent(in) :: & logical, intent(in) :: &
timeinc_old, & guess
timeinc, & real(pReal), intent(in) :: &
loadCaseTime !< remaining time of current load case timeinc_old, &
type(tBoundaryCondition), intent(in) :: & timeinc, &
stress_BC, & loadCaseTime !< remaining time of current load case
deformation_BC type(tBoundaryCondition), intent(in) :: &
real(pReal), dimension(3,3), intent(in) :: rotation_BC stress_BC, &
logical, intent(in) :: & deformation_BC
guess real(pReal), dimension(3,3), intent(in) ::&
PetscErrorCode :: ierr rotation_BC
PetscScalar, pointer :: F(:,:,:,:) PetscErrorCode :: ierr
PetscScalar, dimension(:,:,:,:), pointer :: F
character(len=1024) :: rankStr character(len=32) :: rankStr
call DMDAVecGetArrayF90(da,solution_vec,F,ierr) call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! restart information for spectral solver if (cutBack) then
if (restartWrite) then C_volAvg = C_volAvgLastInc ! QUESTION: where is this required?
write(6,'(/,a)') ' writing converged results for restart' C_minMaxAvg = C_minMaxAvgLastInc ! QUESTION: where is this required?
flush(6) else
write(rankStr,'(a1,i0)')'_',worldrank !--------------------------------------------------------------------------------------------------
call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file ! restart information for spectral solver
write (777,rec=1) F if (restartWrite) then ! QUESTION: where is this logical properly set?
close (777) write(6,'(/,a)') ' writing converged results for restart'
call IO_write_jobRealFile(777,'F_lastInc'//trim(rankStr),size(F_lastInc)) ! writing F_lastInc field to file flush(6)
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 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 write(rankStr,'(a1,i0)')'_',worldrank
F_aim = F_aim_lastInc call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file
F = reshape(F_lastInc, [9,grid(1),grid(2),grid3]) write (777,rec=1) F; close (777)
C_volAvg = C_volAvgLastInc call IO_write_jobRealFile(777,'F_lastInc'//trim(rankStr),size(F_lastInc)) ! writing F_lastInc field to file
else write (777,rec=1) F_lastInc; close (777)
ForwardData = .True. endif
C_volAvgLastInc = C_volAvg
!-------------------------------------------------------------------------------------------------- call CPFEM_age() ! age state and kinematics
! calculate rate for aim call utilities_updateIPcoords(F)
if (deformation_BC%myType=='l') then ! calculate f_aimDot from given L and current F
f_aimDot = deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim) C_volAvgLastInc = C_volAvg
elseif(deformation_BC%myType=='fdot') then ! f_aimDot is prescribed C_minMaxAvgLastInc = C_minMaxAvg
f_aimDot = deformation_BC%maskFloat * deformation_BC%values
elseif(deformation_BC%myType=='f') then ! aim at end of load case is prescribed if (guess) then ! QUESTION: better with a = L ? x:y
f_aimDot = deformation_BC%maskFloat * (deformation_BC%values -F_aim)/loadCaseTime F_aimDot = stress_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old ! initialize with correction based on last inc
endif else
if (guess) f_aimDot = f_aimDot + stress_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old F_aimDot = 0.0_pReal
F_aim_lastInc = F_aim 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 ! update average and local deformation gradients
call utilities_updateIPcoords(F) F_aim = F_aim_lastInc + F_aimDot * timeinc
Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average
timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3])) math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid3])
F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
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)
end subroutine BasicPETSc_forward end subroutine BasicPETSc_forward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

463
src/spectral_mech_Polarisation.f90
View File

@ -43,7 +43,7 @@ module spectral_mech_Polarisation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc. ! stress, stiffness and compliance average etc.
real(pReal), private, dimension(3,3) :: & 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 = math_I3, & !< current prescribed deformation gradient
F_aim_lastInc = math_I3, & !< previous average deformation gradient F_aim_lastInc = math_I3, & !< previous average deformation gradient
F_av = 0.0_pReal, & !< average incompatible def grad field 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_volAvg = 0.0_pReal, & !< current volume average stiffness
C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness
C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 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) S = 0.0_pReal, & !< current compliance (filled up with zeros)
C_scale = 0.0_pReal, & C_scale = 0.0_pReal, &
S_scale = 0.0_pReal S_scale = 0.0_pReal
@ -62,7 +63,7 @@ module spectral_mech_Polarisation
err_BC, & !< deviation from stress BC err_BC, & !< deviation from stress BC
err_curl, & !< RMS of curl of F err_curl, & !< RMS of curl of F
err_div !< RMS of div of P err_div !< RMS of div of P
logical, private :: ForwardData
integer(pInt), private :: & integer(pInt), private :: &
totalIter = 0_pInt !< total iteration in current increment totalIter = 0_pInt !< total iteration in current increment
@ -80,7 +81,7 @@ module spectral_mech_Polarisation
contains 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) !> @todo use sourced allocation, e.g. allocate(Fdot,source = F_lastInc)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine Polarisation_init subroutine Polarisation_init
@ -102,12 +103,15 @@ subroutine Polarisation_init
use numerics, only: & use numerics, only: &
worldrank, & worldrank, &
worldsize worldsize
use homogenization, only: &
materialpoint_F0
use DAMASK_interface, only: & use DAMASK_interface, only: &
getSolverJobName getSolverJobName
use spectral_utilities, only: & use spectral_utilities, only: &
Utilities_constitutiveResponse, & Utilities_constitutiveResponse, &
Utilities_updateGamma, & Utilities_updateGamma, &
Utilities_updateIPcoords Utilities_updateIPcoords, &
wgt
use mesh, only: & use mesh, only: &
grid, & grid, &
grid3 grid3
@ -120,7 +124,11 @@ subroutine Polarisation_init
temp33_Real = 0.0_pReal temp33_Real = 0.0_pReal
PetscErrorCode :: ierr 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), dimension(:), allocatable :: localK
integer(pInt) :: proc integer(pInt) :: proc
character(len=1024) :: rankStr character(len=1024) :: rankStr
@ -164,78 +172,69 @@ subroutine Polarisation_init
grid(1),grid(2),localK, & ! local grid grid(1),grid(2),localK, & ! local grid
da,ierr) ! handle, error da,ierr) ! handle, error
CHKERRQ(ierr) 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) 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) 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) CHKERRQ(ierr)
call SNESSetConvergenceTest(snes,Polarisation_converged,PETSC_NULL_OBJECT,PETSC_NULL_FUNCTION,ierr) call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional cli arguments
CHKERRQ(ierr)
call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init fields ! init fields
call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! places pointer xx_psc on PETSc data call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! places pointer on PETSc data
F => xx_psc(0:8,:,:,:) F => FandF_tau( 0: 8,:,:,:)
F_tau => xx_psc(9:17,:,:,:) F_tau => FandF_tau( 9:17,:,:,:)
restart: if (restartInc > 1_pInt) then 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)') & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'reading values of increment ', restartInc - 1_pInt, ' from file' 'reading values of increment ', restartInc - 1_pInt, ' from file'
flush(6) flush(6)
endif
write(rankStr,'(a1,i0)')'_',worldrank write(rankStr,'(a1,i0)')'_',worldrank
call IO_read_realFile(777,'F'//trim(rankStr),trim(getSolverJobName()),size(F)) call IO_read_realFile(777,'F'//trim(rankStr),trim(getSolverJobName()),size(F))
read (777,rec=1) F read (777,rec=1) F; close (777)
close (777)
call IO_read_realFile(777,'F_lastInc'//trim(rankStr),trim(getSolverJobName()),size(F_lastInc)) call IO_read_realFile(777,'F_lastInc'//trim(rankStr),trim(getSolverJobName()),size(F_lastInc))
read (777,rec=1) F_lastInc read (777,rec=1) F_lastInc; close (777)
close (777)
call IO_read_realFile(777,'F_tau'//trim(rankStr),trim(getSolverJobName()),size(F_tau)) call IO_read_realFile(777,'F_tau'//trim(rankStr),trim(getSolverJobName()),size(F_tau))
read (777,rec=1) F_tau read (777,rec=1) F_tau; close (777)
close (777) call IO_read_realFile(777,'F_tau_lastInc'//trim(rankStr),trim(getSolverJobName()),size(F_tau_lastInc))
call IO_read_realFile(777,'F_tau_lastInc'//trim(rankStr),& read (777,rec=1) F_tau_lastInc; close (777)
trim(getSolverJobName()),size(F_tau_lastInc)) call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(F_aimDot))
read (777,rec=1) F_tau_lastInc read (777,rec=1) F_aimDot; close (777)
close (777) F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F
call IO_read_realFile(777,'F_aim', trim(getSolverJobName()),size(F_aim)) F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc
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)
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_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 = 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 F_tau_lastInc = 2.0_pReal*F_lastInc
endif restart 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_updateIPcoords(reshape(F,shape(F_lastInc)))
call Utilities_constitutiveResponse(F_lastInc, reshape(F,shape(F_lastInc)), & call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
0.0_pReal,P,C_volAvg,C_minMaxAvg,temp33_Real,.false.,math_I3) reshape(F,shape(F_lastInc)), & ! target F
0.0_pReal, & ! time increment
math_I3) ! no rotation of boundary condition
nullify(F) nullify(F)
nullify(F_tau) 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 restartRead: 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 .and. worldrank == 0_pInt) &
write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'reading more values of increment', restartInc - 1_pInt, 'from file' 'reading more values of increment ', restartInc-1_pInt, ' from file'
flush(6) flush(6)
call IO_read_realFile(777,'C_volAvg',trim(getSolverJobName()),size(C_volAvg)) call IO_read_realFile(777,'C_volAvg',trim(getSolverJobName()),size(C_volAvg))
read (777,rec=1) C_volAvg read (777,rec=1) C_volAvg; close (777)
close (777)
call IO_read_realFile(777,'C_volAvgLastInc',trim(getSolverJobName()),size(C_volAvgLastInc)) call IO_read_realFile(777,'C_volAvgLastInc',trim(getSolverJobName()),size(C_volAvgLastInc))
read (777,rec=1) C_volAvgLastInc read (777,rec=1) C_volAvgLastInc; close (777)
close (777)
call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg)) call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg))
read (777,rec=1) C_minMaxAvg read (777,rec=1) C_minMaxAvg; close (777)
close (777)
endif restartRead endif restartRead
call Utilities_updateGamma(C_minMaxAvg,.True.) call Utilities_updateGamma(C_minMaxAvg,.true.)
C_scale = C_minMaxAvg C_scale = C_minMaxAvg
S_scale = math_invSym3333(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 !> @brief solution for the Polarisation scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
type(tSolutionState) function & type(tSolutionState) function Polarisation_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC)
Polarisation_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC)
use IO, only: & use IO, only: &
IO_error IO_error
use numerics, only: & use numerics, only: &
@ -265,13 +263,13 @@ type(tSolutionState) function &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! input data for solution ! 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) :: & character(len=*), intent(in) :: &
incInfoIn 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 real(pReal), dimension(3,3), intent(in) :: rotation_BC
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -304,18 +302,17 @@ type(tSolutionState) function &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! solve BVP ! solve BVP
call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr) call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr); CHKERRQ(ierr)
CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! check convergence ! check convergence
call SNESGetConvergedReason(snes,reason,ierr) call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr)
CHKERRQ(ierr)
Polarisation_solution%converged = reason > 0
Polarisation_solution%iterationsNeeded = totalIter
Polarisation_solution%termIll = terminallyIll Polarisation_solution%termIll = terminallyIll
terminallyIll = .false. terminallyIll = .false.
if (reason == -4) call IO_error(893_pInt) if (reason == -4) call IO_error(893_pInt) ! MPI error
if (reason < 1) Polarisation_solution%converged = .false.
Polarisation_solution%iterationsNeeded = totalIter
end function Polarisation_solution end function Polarisation_solution
@ -330,8 +327,8 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
polarAlpha, & polarAlpha, &
polarBeta polarBeta
use mesh, only: & use mesh, only: &
grid3, & grid, &
grid grid3
use IO, only: & use IO, only: &
IO_intOut IO_intOut
use math, only: & use math, only: &
@ -340,6 +337,10 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
math_mul3333xx33, & math_mul3333xx33, &
math_invSym3333, & math_invSym3333, &
math_mul33x33 math_mul33x33
use debug, only: &
debug_level, &
debug_spectral, &
debug_spectralRotation
use spectral_utilities, only: & use spectral_utilities, only: &
wgt, & wgt, &
tensorField_real, & tensorField_real, &
@ -349,27 +350,17 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
Utilities_constitutiveResponse, & Utilities_constitutiveResponse, &
Utilities_divergenceRMS, & Utilities_divergenceRMS, &
Utilities_curlRMS Utilities_curlRMS
use debug, only: &
debug_level, &
debug_spectral, &
debug_spectralRotation
use homogenization, only: & use homogenization, only: &
materialpoint_dPdF materialpoint_dPdF
use FEsolving, only: & use FEsolving, only: &
terminallyIll terminallyIll
implicit none implicit none
!-------------------------------------------------------------------------------------------------- DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in
! strange syntax in the next line because otherwise macros expand beyond 132 character limit PetscScalar, &
DMDALocalInfo, dimension(& target, dimension(3,3,2, XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: x_scal
DMDA_LOCAL_INFO_SIZE) :: & PetscScalar, &
in 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(:,:,:,:,:) :: & PetscScalar, pointer, dimension(:,:,:,:,:) :: &
F, & F, &
F_tau, & F_tau, &
@ -387,33 +378,33 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
SNESGetNumberFunctionEvals, & SNESGetNumberFunctionEvals, &
SNESGetIterationNumber SNESGetIterationNumber
F => x_scal(1:3,1:3,1,& F => x_scal(1:3,1:3,1,&
XG_RANGE,YG_RANGE,ZG_RANGE) XG_RANGE,YG_RANGE,ZG_RANGE)
F_tau => x_scal(1:3,1:3,2,& F_tau => x_scal(1:3,1:3,2,&
XG_RANGE,YG_RANGE,ZG_RANGE) XG_RANGE,YG_RANGE,ZG_RANGE)
residual_F => f_scal(1:3,1:3,1,& residual_F => f_scal(1:3,1:3,1,&
X_RANGE,Y_RANGE,Z_RANGE) X_RANGE, Y_RANGE, Z_RANGE)
residual_F_tau => f_scal(1:3,1:3,2,& residual_F_tau => f_scal(1:3,1:3,2,&
X_RANGE,Y_RANGE,Z_RANGE) X_RANGE, Y_RANGE, Z_RANGE)
call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr)
call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr)
F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt 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) 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 call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr)
newIteration: if(totalIter <= PETScIter) then 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 totalIter = totalIter + 1_pInt
write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') &
' @ Iteration ', itmin, '≤',totalIter, '≤', itmax trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) ' 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 =', & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
math_transpose33(F_aim) ' deformation gradient aim =', math_transpose33(F_aim)
flush(6) flush(6)
endif newIteration endif newIteration
@ -440,18 +431,16 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate constitutive response ! evaluate constitutive response
P_avLastEval = P_av P_avLastEval = P_av
call Utilities_constitutiveResponse(F_lastInc,F - residual_F_tau/polarBeta,params%timeinc, & call Utilities_constitutiveResponse(residual_F,P_av,C_volAvg,C_minMaxAvg, &
residual_F,C_volAvg,C_minMaxAvg,P_av,ForwardData,params%rotation_BC) 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) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
ForwardData = .False.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate divergence ! calculate divergence
tensorField_real = 0.0_pReal 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() call utilities_FFTtensorForward()
err_div = Utilities_divergenceRMS() err_div = Utilities_divergenceRMS() !< root mean squared error in divergence of stress
call utilities_FFTtensorBackward()
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! constructing residual ! 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 tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F
call utilities_FFTtensorForward() call utilities_FFTtensorForward()
err_curl = Utilities_curlRMS() err_curl = Utilities_curlRMS()
call utilities_FFTtensorBackward()
nullify(F)
nullify(F_tau)
nullify(residual_F)
nullify(residual_F_tau)
end subroutine Polarisation_formResidual end subroutine Polarisation_formResidual
@ -487,8 +479,8 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,
err_div_tolAbs, & err_div_tolAbs, &
err_curl_tolRel, & err_curl_tolRel, &
err_curl_tolAbs, & err_curl_tolAbs, &
err_stress_tolAbs, & err_stress_tolRel, &
err_stress_tolRel err_stress_tolAbs
use math, only: & use math, only: &
math_mul3333xx33 math_mul3333xx33
use FEsolving, only: & use FEsolving, only: &
@ -507,24 +499,24 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,
real(pReal) :: & real(pReal) :: &
curlTol, & curlTol, &
divTol, & divTol, &
BC_tol BCTol
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress BC handling ! stress BC handling
F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc 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) + & 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 mask_stress * (P_av-params%stress_BC))) ! mask = 0.0 for no bc
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! error calculation ! error calculation
curlTol = max(maxval(abs(F_aim-math_I3))*err_curl_tolRel,err_curl_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) 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) BCTol = max(maxval(abs(P_av)) *err_stress_tolRel,err_stress_tolAbs)
converged: if ((totalIter >= itmin .and. & converged: if ((totalIter >= itmin .and. &
all([ err_div/divTol, & all([ err_div /divTol, &
err_curl/curlTol, & err_curl/curlTol, &
err_BC/BC_tol ] < 1.0_pReal)) & err_BC /BCTol ] < 1.0_pReal)) &
.or. terminallyIll) then .or. terminallyIll) then
reason = 1 reason = 1
elseif (totalIter >= itmax) then converged elseif (totalIter >= itmax) then converged
@ -536,12 +528,12 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report ! report
write(6,'(1/,a)') ' ... reporting .............................................................' write(6,'(1/,a)') ' ... reporting .............................................................'
write(6,'(/,a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', & write(6,'(/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', &
err_curl/curlTol,' (',err_curl,' -, tol =',curlTol,')' err_div/divTol, ' (',err_div, ' / m, tol = ',divTol,')'
write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', &
err_div/divTol, ' (',err_div, ' / m, tol =',divTol,')' err_curl/curlTol,' (',err_curl,' -, tol = ',curlTol,')'
write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', & write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', &
err_BC/BC_tol, ' (',err_BC, ' Pa, tol =',BC_tol,')' err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')'
write(6,'(/,a)') ' ===========================================================================' write(6,'(/,a)') ' ==========================================================================='
flush(6) flush(6)
@ -549,6 +541,8 @@ end subroutine Polarisation_converged
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief forwarding routine !> @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) subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC)
use math, only: & use math, only: &
@ -558,133 +552,140 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,deformati
math_rotate_backward33 math_rotate_backward33
use numerics, only: & use numerics, only: &
worldrank worldrank
use homogenization, only: &
materialpoint_F0
use mesh, only: & use mesh, only: &
grid3, & grid, &
grid grid3
use spectral_utilities, only: & use CPFEM2, only: &
Utilities_calculateRate, & CPFEM_age
Utilities_forwardField, & use spectral_utilities, only: &
Utilities_updateIPcoords, & Utilities_calculateRate, &
tBoundaryCondition, & Utilities_forwardField, &
cutBack Utilities_updateIPcoords, &
use IO, only: & tBoundaryCondition, &
IO_write_JobRealFile cutBack
use FEsolving, only: & use IO, only: &
restartWrite IO_write_JobRealFile
use FEsolving, only: &
restartWrite
implicit none implicit none
real(pReal), intent(in) :: & logical, intent(in) :: &
timeinc_old, & guess
timeinc, & real(pReal), intent(in) :: &
loadCaseTime !< remaining time of current load case timeinc_old, &
type(tBoundaryCondition), intent(in) :: & timeinc, &
stress_BC, & loadCaseTime !< remaining time of current load case
deformation_BC type(tBoundaryCondition), intent(in) :: &
real(pReal), dimension(3,3), intent(in) :: rotation_BC stress_BC, &
logical, intent(in) :: & deformation_BC
guess real(pReal), dimension(3,3), intent(in) ::&
PetscErrorCode :: ierr rotation_BC
PetscScalar, dimension(:,:,:,:), pointer :: xx_psc, F, F_tau PetscErrorCode :: ierr
integer(pInt) :: i, j, k PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau, F, F_tau
real(pReal), dimension(3,3) :: F_lambda33 integer(pInt) :: i, j, k
character(len=1024) :: rankStr real(pReal), dimension(3,3) :: F_lambda33
character(len=32) :: rankStr
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update coordinates and rate and forward last inc ! update coordinates and rate and forward last inc
call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr) call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
F => xx_psc(0:8,:,:,:) F => FandF_tau( 0: 8,:,:,:)
F_tau => xx_psc(9:17,:,:,:) F_tau => FandF_tau( 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 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 if (worldrank == 0_pInt) then
F_aim = F_aim_lastInc call IO_write_jobRealFile(777,'C_volAvg',size(C_volAvg))
F_tau= reshape(F_tau_lastInc,[9,grid(1),grid(2),grid3]) write (777,rec=1) C_volAvg; close(777)
F = reshape(F_lastInc, [9,grid(1),grid(2),grid3]) call IO_write_jobRealFile(777,'C_volAvgLastInc',size(C_volAvgLastInc))
C_volAvg = C_volAvgLastInc write (777,rec=1) C_volAvgLastInc; close(777)
else call IO_write_jobRealFile(777,'F_aimDot',size(F_aimDot))
ForwardData = .True. write (777,rec=1) F_aimDot; close(777)
C_volAvgLastInc = C_volAvg endif
!--------------------------------------------------------------------------------------------------
! calculate rate for aim write(rankStr,'(a1,i0)')'_',worldrank
if (deformation_BC%myType=='l') then ! calculate f_aimDot from given L and current F call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file
f_aimDot = deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim) write (777,rec=1) F; close (777)
elseif(deformation_BC%myType=='fdot') then ! f_aimDot is prescribed call IO_write_jobRealFile(777,'F_lastInc'//trim(rankStr),size(F_lastInc)) ! writing F_lastInc field to file
f_aimDot = deformation_BC%maskFloat * deformation_BC%values write (777,rec=1) F_lastInc; close (777)
elseif(deformation_BC%myType=='f') then ! aim at end of load case is prescribed call IO_write_jobRealFile(777,'F_tau'//trim(rankStr),size(F_tau)) ! writing deformation gradient field to file
f_aimDot = deformation_BC%maskFloat * (deformation_BC%values -F_aim)/loadCaseTime write (777,rec=1) F_tau; close (777)
endif call IO_write_jobRealFile(777,'F_tau_lastInc'//trim(rankStr),size(F_tau_lastInc)) ! writing F_tau_lastInc field to file
if (guess) f_aimDot = f_aimDot + stress_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old write (777,rec=1) F_tau_lastInc; close (777)
F_aim_lastInc = F_aim 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 ! update average and local deformation gradients
call utilities_updateIPcoords(F) F_aim = F_aim_lastInc + F_aimDot * timeinc
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
F_aim = F_aim + 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])
! update local deformation gradient if (guess) then
F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! ensure that it matches rotated F_aim F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), &
math_rotate_backward33(F_aim,rotation_BC)), & [9,grid(1),grid(2),grid3]) ! does not have any average value as boundary condition
[9,grid(1),grid(2),grid3]) else
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
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) 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 = 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, & F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, &
math_mul3333xx33(C_scale,& math_mul3333xx33(C_scale,&
math_mul33x33(math_transpose33(F_lambda33),& math_mul33x33(math_transpose33(F_lambda33),&
F_lambda33) -math_I3))*0.5_pReal)& F_lambda33)-math_I3))*0.5_pReal)&
+ math_I3 + math_I3
F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(1:9,i,j,k) F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(1:9,i,j,k)
enddo; enddo; enddo enddo; enddo; enddo
endif 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 end subroutine Polarisation_forward

236
src/spectral_utilities.f90 Normal file → Executable file
View File

@ -16,7 +16,7 @@ module spectral_utilities
#include <petsc/finclude/petscsys.h> #include <petsc/finclude/petscsys.h>
include 'fftw3-mpi.f03' 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, parameter :: maxPhaseFields = 2_pInt
integer(pInt), public :: nActiveFields = 0_pInt integer(pInt), public :: nActiveFields = 0_pInt
@ -145,8 +145,7 @@ module spectral_utilities
FIELD_UNDEFINED_ID, & FIELD_UNDEFINED_ID, &
FIELD_MECH_ID, & FIELD_MECH_ID, &
FIELD_THERMAL_ID, & FIELD_THERMAL_ID, &
FIELD_DAMAGE_ID, & FIELD_DAMAGE_ID
utilities_calcPlasticity
private :: & private :: &
utilities_getFreqDerivative utilities_getFreqDerivative
@ -414,8 +413,7 @@ subroutine utilities_updateGamma(C,saveReference)
write(6,'(/,a)') ' writing reference stiffness to file' write(6,'(/,a)') ' writing reference stiffness to file'
flush(6) flush(6)
call IO_write_jobRealFile(777,'C_ref',size(C_ref)) call IO_write_jobRealFile(777,'C_ref',size(C_ref))
write (777,rec=1) C_ref write (777,rec=1) C_ref; close(777)
close(777)
endif endif
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 call math_invert(size_reduced, c_reduced, s_reduced, errmatinv) ! invert reduced stiffness
if (any(IEEE_is_NaN(s_reduced))) errmatinv = .true. 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 temp99_Real = 0.0_pReal ! fill up compliance with zeros
k = 0_pInt k = 0_pInt
do n = 1_pInt,9_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) sTimesC = matmul(c_reduced,s_reduced)
do m=1_pInt, size_reduced do m=1_pInt, size_reduced
do n=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 errmatinv = errmatinv &
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 .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
enddo enddo
if(debugGeneral .or. errmatinv) then if (debugGeneral .or. errmatinv) then
write(formatString, '(I16.16)') size_reduced write(formatString, '(i2)') size_reduced
formatString = '(/,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))' formatString = '(/,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))'
write(6,trim(formatString),advance='no') ' C * S (load) ', & write(6,trim(formatString),advance='no') ' C * S (load) ', &
transpose(matmul(c_reduced,s_reduced)) transpose(matmul(c_reduced,s_reduced))
write(6,trim(formatString),advance='no') ' S (load) ', transpose(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 endif
if(errmatinv) call IO_error(error_ID=400_pInt,ext_msg='utilities_maskedCompliance')
deallocate(c_reduced) deallocate(c_reduced)
deallocate(s_reduced) deallocate(s_reduced)
deallocate(sTimesC) deallocate(sTimesC)
else else
temp99_real = 0.0_pReal temp99_real = 0.0_pReal
endif endif
if(debugGeneral) & if(debugGeneral) then
write(6,'(/,a,/,9(9(2x,f12.7,1x)/),/)',advance='no') ' Masked Compliance (load) * GPa =', & write(6,'(/,a,/,9(9(2x,f10.5,1x)/),/)',advance='no') &
transpose(temp99_Real*1.e9_pReal) ' Masked Compliance (load) / GPa =', transpose(temp99_Real*1.e-9_pReal)
flush(6) flush(6)
endif
utilities_maskedCompliance = math_Plain99to3333(temp99_Real) utilities_maskedCompliance = math_Plain99to3333(temp99_Real)
end function utilities_maskedCompliance 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, & subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
P,C_volAvg,C_minmaxAvg,P_av,forwardData,rotation_BC) F,timeinc,rotation_BC)
use IO, only: & use IO, only: &
IO_error IO_error
use debug, only: & use debug, only: &
@ -941,31 +941,22 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
use mesh, only: & use mesh, only: &
grid,& grid,&
grid3 grid3
use FEsolving, only: &
restartWrite
use CPFEM2, only: &
CPFEM_general
use homogenization, only: & use homogenization, only: &
materialpoint_F0, &
materialpoint_F, & materialpoint_F, &
materialpoint_P, & materialpoint_P, &
materialpoint_dPdF materialpoint_dPdF, &
materialpoint_stressAndItsTangent
implicit none 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,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) :: P_av !< average PK stress
real(pReal),intent(out), dimension(3,3,grid(1),grid(2),grid3) :: P !< PK stress real(pReal),intent(out), dimension(3,3,grid(1),grid(2),grid3) :: P !< PK stress
logical :: & real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: F !< deformation gradient target !< previous deformation gradient
age real(pReal), intent(in) :: timeinc !< loading time
real(pReal), intent(in), dimension(3,3) :: rotation_BC !< rotation of load frame
integer(pInt) :: & integer(pInt) :: &
j,k,ierr j,k,ierr
real(pReal), dimension(3,3,3,3) :: max_dPdF, min_dPdF 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 ......................................' write(6,'(/,a)') ' ... evaluating constitutive response ......................................'
flush(6) flush(6)
age = .False.
materialpoint_F = reshape(F,[3,3,1,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field
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
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate bounds of det(F) and report ! calculate bounds of det(F) and report
if(debugGeneral) then if(debugGeneral) then
@ -1003,7 +986,19 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
flush(6) flush(6)
endif 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 = 0.0_pReal
max_dPdF_norm = 0.0_pReal max_dPdF_norm = 0.0_pReal
@ -1021,157 +1016,24 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
end do end do
call MPI_Allreduce(MPI_IN_PLACE,max_dPdF,81,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) 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) 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_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 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 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 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 !> @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: & use mesh, only: &
grid3, & grid3, &
grid grid
@ -1179,17 +1041,17 @@ pure function utilities_calculateRate(avRate,timeinc_old,guess,field_lastInc,fie
implicit none implicit none
real(pReal), intent(in), dimension(3,3) :: avRate !< homogeneous addon real(pReal), intent(in), dimension(3,3) :: avRate !< homogeneous addon
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
timeinc_old !< timeinc of last step dt !< timeinc between field0 and field
logical, intent(in) :: & 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) :: & 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 field !< data of current step
real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: & real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: &
utilities_calculateRate utilities_calculateRate
if (guess) then if (heterogeneous) then
utilities_calculateRate = (field-field_lastInc) / timeinc_old utilities_calculateRate = (field-field0) / dt
else else
utilities_calculateRate = spread(spread(spread(avRate,3,grid(1)),4,grid(2)),5,grid3) utilities_calculateRate = spread(spread(spread(avRate,3,grid(1)),4,grid(2)),5,grid3)
endif endif