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

consistent indentation and line-spacings in reporting

This commit is contained in:
Philip Eisenlohr 2021-11-15 12:35:44 -05:00
parent 292449aca9
commit da9fdf53d2
60 changed files with 725 additions and 725 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/CPFEM.f90
View File

@ -103,7 +103,7 @@ subroutine CPFEM_init
class(tNode), pointer :: & class(tNode), pointer :: &
debug_CPFEM debug_CPFEM
print'(/,a)', ' <<<+- CPFEM init -+>>>'; flush(IO_STDOUT) print'(/,1x,a)', '<<<+- CPFEM init -+>>>'; flush(IO_STDOUT)
allocate(CPFEM_cs( 6,discretization_nIPs,discretization_Nelems), source= 0.0_pReal) allocate(CPFEM_cs( 6,discretization_nIPs,discretization_Nelems), source= 0.0_pReal)
allocate(CPFEM_dcsdE( 6,6,discretization_nIPs,discretization_Nelems), source= 0.0_pReal) allocate(CPFEM_dcsdE( 6,6,discretization_nIPs,discretization_Nelems), source= 0.0_pReal)

2
src/CPFEM2.f90
View File

@ -81,7 +81,7 @@ subroutine CPFEM_init
integer(HID_T) :: fileHandle integer(HID_T) :: fileHandle
print'(/,a)', ' <<<+- CPFEM init -+>>>'; flush(IO_STDOUT) print'(/,1x,a)', '<<<+- CPFEM init -+>>>'; flush(IO_STDOUT)
if (interface_restartInc > 0) then if (interface_restartInc > 0) then

2
src/DAMASK_Marc.f90
View File

@ -46,7 +46,7 @@ subroutine DAMASK_interface_init
integer :: ierr integer :: ierr
character(len=pPathLen) :: wd character(len=pPathLen) :: wd
print'(/,a)', ' <<<+- DAMASK_marc init -+>>>' print'(/,1x,a)', '<<<+- DAMASK_marc init -+>>>'
print*, 'Roters et al., Computational Materials Science 158:420–478, 2019' print*, 'Roters et al., Computational Materials Science 158:420–478, 2019'
print*, 'https://doi.org/10.1016/j.commatsci.2018.04.030' print*, 'https://doi.org/10.1016/j.commatsci.2018.04.030'

2
src/DAMASK_interface.f90
View File

@ -70,7 +70,7 @@ subroutine DAMASK_interface_init
external :: & external :: &
quit quit
print'(/,a)', ' <<<+- DAMASK_interface init -+>>>' print'(/,1x,a)', '<<<+- DAMASK_interface init -+>>>'
if(worldrank == 0) open(OUTPUT_UNIT, encoding='UTF-8') ! for special characters in output if(worldrank == 0) open(OUTPUT_UNIT, encoding='UTF-8') ! for special characters in output

2
src/HDF5_utilities.f90
View File

@ -109,7 +109,7 @@ subroutine HDF5_utilities_init
integer(SIZE_T) :: typeSize integer(SIZE_T) :: typeSize
print'(/,a)', ' <<<+- HDF5_Utilities init -+>>>' print'(/,1x,a)', '<<<+- HDF5_Utilities init -+>>>'
call h5open_f(hdferr) call h5open_f(hdferr)

2
src/IO.f90
View File

@ -56,7 +56,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine IO_init subroutine IO_init
print'(/,a)', ' <<<+- IO init -+>>>'; flush(IO_STDOUT) print'(/,1x,a)', '<<<+- IO init -+>>>'; flush(IO_STDOUT)
call selfTest call selfTest

2
src/YAML_types.f90
View File

@ -187,7 +187,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine YAML_types_init subroutine YAML_types_init
print'(/,a)', ' <<<+- YAML_types init -+>>>' print'(/,1x,a)', '<<<+- YAML_types init -+>>>'
call selfTest call selfTest

2
src/base64.f90
View File

@ -27,7 +27,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine base64_init subroutine base64_init
print'(/,a)', ' <<<+- base64 init -+>>>'; flush(IO_STDOUT) print'(/,1x,a)', '<<<+- base64 init -+>>>'; flush(IO_STDOUT)
call selfTest call selfTest

25
src/config.f90
View File

@ -30,8 +30,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine config_init subroutine config_init
print'(/,a)', ' <<<+- config init -+>>>'; flush(IO_STDOUT) print'(/,1x,a)', '<<<+- config init -+>>>'; flush(IO_STDOUT)
call parse_material call parse_material
call parse_numerics call parse_numerics
@ -50,15 +49,15 @@ subroutine parse_material()
inquire(file='material.yaml',exist=fileExists) inquire(file='material.yaml',exist=fileExists)
if(.not. fileExists) call IO_error(100,ext_msg='material.yaml') if (.not. fileExists) call IO_error(100,ext_msg='material.yaml')
if (worldrank == 0) then if (worldrank == 0) then
print*, 'reading material.yaml'; flush(IO_STDOUT) print'(/,1x,a)', 'reading material.yaml'; flush(IO_STDOUT)
fileContent = IO_read('material.yaml') fileContent = IO_read('material.yaml')
call results_openJobFile(parallel=.false.) call results_openJobFile(parallel=.false.)
call results_writeDataset_str(fileContent,'setup','material.yaml','main configuration') call results_writeDataset_str(fileContent,'setup','material.yaml','main configuration')
call results_closeJobFile call results_closeJobFile
endif end if
call parallelization_bcast_str(fileContent) call parallelization_bcast_str(fileContent)
config_material => YAML_parse_str(fileContent) config_material => YAML_parse_str(fileContent)
@ -81,19 +80,19 @@ subroutine parse_numerics()
if (fileExists) then if (fileExists) then
if (worldrank == 0) then if (worldrank == 0) then
print*, 'reading numerics.yaml'; flush(IO_STDOUT) print'(1x,a)', 'reading numerics.yaml'; flush(IO_STDOUT)
fileContent = IO_read('numerics.yaml') fileContent = IO_read('numerics.yaml')
if (len(fileContent) > 0) then if (len(fileContent) > 0) then
call results_openJobFile(parallel=.false.) call results_openJobFile(parallel=.false.)
call results_writeDataset_str(fileContent,'setup','numerics.yaml','numerics configuration') call results_writeDataset_str(fileContent,'setup','numerics.yaml','numerics configuration')
call results_closeJobFile call results_closeJobFile
endif end if
endif end if
call parallelization_bcast_str(fileContent) call parallelization_bcast_str(fileContent)
config_numerics => YAML_parse_str(fileContent) config_numerics => YAML_parse_str(fileContent)
endif end if
end subroutine parse_numerics end subroutine parse_numerics
@ -113,19 +112,19 @@ subroutine parse_debug()
if (fileExists) then if (fileExists) then
if (worldrank == 0) then if (worldrank == 0) then
print*, 'reading debug.yaml'; flush(IO_STDOUT) print'(1x,a)', 'reading debug.yaml'; flush(IO_STDOUT)
fileContent = IO_read('debug.yaml') fileContent = IO_read('debug.yaml')
if (len(fileContent) > 0) then if (len(fileContent) > 0) then
call results_openJobFile(parallel=.false.) call results_openJobFile(parallel=.false.)
call results_writeDataset_str(fileContent,'setup','debug.yaml','debug configuration') call results_writeDataset_str(fileContent,'setup','debug.yaml','debug configuration')
call results_closeJobFile call results_closeJobFile
endif end if
endif end if
call parallelization_bcast_str(fileContent) call parallelization_bcast_str(fileContent)
config_debug => YAML_parse_str(fileContent) config_debug => YAML_parse_str(fileContent)
endif end if
end subroutine parse_debug end subroutine parse_debug

2
src/discretization.f90
View File

@ -49,7 +49,7 @@ subroutine discretization_init(materialAt,&
integer, optional, intent(in) :: & integer, optional, intent(in) :: &
sharedNodesBegin !< index of first node shared among different processes (MPI) sharedNodesBegin !< index of first node shared among different processes (MPI)
print'(/,a)', ' <<<+- discretization init -+>>>'; flush(6) print'(/,1x,a)', '<<<+- discretization init -+>>>'; flush(6)
discretization_Nelems = size(materialAt,1) discretization_Nelems = size(materialAt,1)
discretization_nIPs = size(IPcoords0,2)/discretization_Nelems discretization_nIPs = size(IPcoords0,2)/discretization_Nelems

2
src/element.f90
View File

@ -923,7 +923,7 @@ subroutine tElement_init(self,elemType)
self%nIPneighbors = size(self%IPneighbor,1) self%nIPneighbors = size(self%IPneighbor,1)
print'(/,a)', ' <<<+- element_init -+>>>'; flush(IO_STDOUT) print'(/,1x,a)', '<<<+- element_init -+>>>'; flush(IO_STDOUT)
print*, 'element type: ',self%elemType print*, 'element type: ',self%elemType
print*, ' geom type: ',self%geomType print*, ' geom type: ',self%geomType

52
src/grid/DAMASK_grid.f90
View File

@ -113,10 +113,10 @@ program DAMASK_grid
! init DAMASK (all modules) ! init DAMASK (all modules)
call CPFEM_initAll call CPFEM_initAll
print'(/,a)', ' <<<+- DAMASK_grid init -+>>>'; flush(IO_STDOUT) print'(/,1x,a)', '<<<+- DAMASK_grid init -+>>>'; flush(IO_STDOUT)
print*, 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019' print'(/,1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019'
print*, 'https://doi.org/10.1007/978-981-10-6855-3_80' print'( 1x,a)', 'https://doi.org/10.1007/978-981-10-6855-3_80'
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
@ -227,12 +227,12 @@ program DAMASK_grid
loadCases(l)%estimate_rate = (load_step%get_asBool('estimate_rate',defaultVal=.true.) .and. l>1) loadCases(l)%estimate_rate = (load_step%get_asBool('estimate_rate',defaultVal=.true.) .and. l>1)
reportAndCheck: if (worldrank == 0) then reportAndCheck: if (worldrank == 0) then
print'(/,a,i0)', ' load case: ', l print'(/,1x,a,1x,i0)', 'load case:', l
print*, ' estimate_rate:', loadCases(l)%estimate_rate print'(2x,a,1x,l1)', 'estimate_rate:', loadCases(l)%estimate_rate
if (loadCases(l)%deformation%myType == 'F') then if (loadCases(l)%deformation%myType == 'F') then
print*, ' F:' print'(2x,a)', 'F:'
else else
print*, ' '//loadCases(l)%deformation%myType//' / 1/s:' print'(2x,a)', loadCases(l)%deformation%myType//' / 1/s:'
endif endif
do i = 1, 3; do j = 1, 3 do i = 1, 3; do j = 1, 3
if (loadCases(l)%deformation%mask(i,j)) then if (loadCases(l)%deformation%mask(i,j)) then
@ -246,8 +246,8 @@ program DAMASK_grid
if (any(.not.(loadCases(l)%stress%mask .or. transpose(loadCases(l)%stress%mask)) .and. (math_I3<1))) & if (any(.not.(loadCases(l)%stress%mask .or. transpose(loadCases(l)%stress%mask)) .and. (math_I3<1))) &
errorID = 838 ! no rotation is allowed by stress BC errorID = 838 ! no rotation is allowed by stress BC
if (loadCases(l)%stress%myType == 'P') print*, ' P / MPa:' if (loadCases(l)%stress%myType == 'P') print'(2x,a)', 'P / MPa:'
if (loadCases(l)%stress%myType == 'dot_P') print*, ' dot_P / MPa/s:' if (loadCases(l)%stress%myType == 'dot_P') print'(2x,a)', 'dot_P / MPa/s:'
if (loadCases(l)%stress%myType /= '') then if (loadCases(l)%stress%myType /= '') then
do i = 1, 3; do j = 1, 3 do i = 1, 3; do j = 1, 3
@ -270,15 +270,15 @@ program DAMASK_grid
if (loadCases(l)%f_restart < 1) errorID = 839 if (loadCases(l)%f_restart < 1) errorID = 839
if (dEq(loadCases(l)%r,1.0_pReal,1.e-9_pReal)) then if (dEq(loadCases(l)%r,1.0_pReal,1.e-9_pReal)) then
print'(a)', ' r: 1 (constant step width)' print'(2x,a)', 'r: 1 (constant step width)'
else else
print'(a,f0.3)', ' r: ', loadCases(l)%r print'(2x,a,1x,f0.3)', 'r:', loadCases(l)%r
endif endif
print'(a,f0.3)', ' t: ', loadCases(l)%t print'(2x,a,1x,f0.3)', 't:', loadCases(l)%t
print'(a,i0)', ' N: ', loadCases(l)%N print'(2x,a,1x,i0)', 'N:', loadCases(l)%N
print'(a,i0)', ' f_out: ', loadCases(l)%f_out print'(2x,a,1x,i0)', 'f_out:', loadCases(l)%f_out
if (loadCases(l)%f_restart < huge(0)) & if (loadCases(l)%f_restart < huge(0)) &
print'(a,i0)', ' f_restart: ', loadCases(l)%f_restart print'(2x,a,1x,i0)', 'f_restart:', loadCases(l)%f_restart
if (errorID > 0) call IO_error(error_ID = errorID, el = l) if (errorID > 0) call IO_error(error_ID = errorID, el = l)
@ -317,7 +317,7 @@ program DAMASK_grid
endif endif
writeUndeformed: if (interface_restartInc < 1) then writeUndeformed: if (interface_restartInc < 1) then
print'(/,a)', ' ... writing initial configuration to file ........................' print'(/,1x,a)', '... writing initial configuration to file ........................'
flush(IO_STDOUT) flush(IO_STDOUT)
call CPFEM_results(0,0.0_pReal) call CPFEM_results(0,0.0_pReal)
endif writeUndeformed endif writeUndeformed
@ -353,8 +353,8 @@ program DAMASK_grid
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report begin of new step ! report begin of new step
print'(/,a)', ' ###########################################################################' print'(/,1x,a)', '###########################################################################'
print'(1x,a,es12.5,6(a,i0))', & print'(1x,a,1x,es12.5,6(a,i0))', &
'Time', t, & 'Time', t, &
's: Increment ', inc,'/',loadCases(l)%N,& 's: Increment ', inc,'/',loadCases(l)%N,&
'-', stepFraction,'/',subStepFactor**cutBackLevel,& '-', stepFraction,'/',subStepFactor**cutBackLevel,&
@ -379,7 +379,7 @@ program DAMASK_grid
case(FIELD_DAMAGE_ID); call grid_damage_spectral_forward(cutBack) case(FIELD_DAMAGE_ID); call grid_damage_spectral_forward(cutBack)
end select end select
enddo enddo
if(.not. cutBack) call CPFEM_forward if (.not. cutBack) call CPFEM_forward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! solve fields ! solve fields
@ -425,7 +425,7 @@ program DAMASK_grid
cutBackLevel = cutBackLevel + 1 cutBackLevel = cutBackLevel + 1
t = t - Delta_t t = t - Delta_t
Delta_t = Delta_t/real(subStepFactor,pReal) ! cut timestep Delta_t = Delta_t/real(subStepFactor,pReal) ! cut timestep
print'(/,a)', ' cutting back ' print'(/,1x,a)', 'cutting back '
else ! no more options to continue else ! no more options to continue
if (worldrank == 0) close(statUnit) if (worldrank == 0) close(statUnit)
call IO_error(950) call IO_error(950)
@ -436,26 +436,26 @@ program DAMASK_grid
cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc
if (all(solres(:)%converged)) then if (all(solres(:)%converged)) then
print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' converged' print'(/,1x,a,i0,a)', 'increment ', totalIncsCounter, ' converged'
else else
print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' NOT converged' print'(/,1x,a,i0,a)', 'increment ', totalIncsCounter, ' NOT converged'
endif; flush(IO_STDOUT) endif; flush(IO_STDOUT)
call MPI_Allreduce(interface_SIGUSR1,signal,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr) call MPI_Allreduce(interface_SIGUSR1,signal,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr)
if (ierr /= 0) error stop 'MPI error' if (ierr /= 0) error stop 'MPI error'
if (mod(inc,loadCases(l)%f_out) == 0 .or. signal) then if (mod(inc,loadCases(l)%f_out) == 0 .or. signal) then
print'(1/,a)', ' ... writing results to file ......................................' print'(/,1x,a)', '... writing results to file ......................................'
flush(IO_STDOUT) flush(IO_STDOUT)
call CPFEM_results(totalIncsCounter,t) call CPFEM_results(totalIncsCounter,t)
endif endif
if(signal) call interface_setSIGUSR1(.false.) if (signal) call interface_setSIGUSR1(.false.)
call MPI_Allreduce(interface_SIGUSR2,signal,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr) call MPI_Allreduce(interface_SIGUSR2,signal,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr)
if (ierr /= 0) error stop 'MPI error' if (ierr /= 0) error stop 'MPI error'
if (mod(inc,loadCases(l)%f_restart) == 0 .or. signal) then if (mod(inc,loadCases(l)%f_restart) == 0 .or. signal) then
call mechanical_restartWrite call mechanical_restartWrite
call CPFEM_restartWrite call CPFEM_restartWrite
endif endif
if(signal) call interface_setSIGUSR2(.false.) if (signal) call interface_setSIGUSR2(.false.)
call MPI_Allreduce(interface_SIGTERM,signal,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr) call MPI_Allreduce(interface_SIGTERM,signal,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr)
if (ierr /= 0) error stop 'MPI error' if (ierr /= 0) error stop 'MPI error'
if (signal) exit loadCaseLooping if (signal) exit loadCaseLooping
@ -468,7 +468,7 @@ program DAMASK_grid
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report summary of whole calculation ! report summary of whole calculation
print'(/,a)', ' ###########################################################################' print'(/,1x,a)', '###########################################################################'
if (worldrank == 0) close(statUnit) if (worldrank == 0) close(statUnit)
call quit(0) ! no complains ;) call quit(0) ! no complains ;)

80
src/grid/discretization_grid.f90
View File

@ -72,10 +72,10 @@ subroutine discretization_grid_init(restart)
fileContent, fname fileContent, fname
print'(/,a)', ' <<<+- discretization_grid init -+>>>'; flush(IO_STDOUT) print'(/,1x,a)', '<<<+- discretization_grid init -+>>>'; flush(IO_STDOUT)
if(worldrank == 0) then if (worldrank == 0) then
fileContent = IO_read(interface_geomFile) fileContent = IO_read(interface_geomFile)
call readVTI(grid,geomSize,origin,materialAt_global,fileContent) call readVTI(grid,geomSize,origin,materialAt_global,fileContent)
fname = interface_geomFile fname = interface_geomFile
@ -85,7 +85,7 @@ subroutine discretization_grid_init(restart)
call results_closeJobFile call results_closeJobFile
else else
allocate(materialAt_global(0)) ! needed for IntelMPI allocate(materialAt_global(0)) ! needed for IntelMPI
endif end if
call MPI_Bcast(grid,3,MPI_INTEGER,0,MPI_COMM_WORLD, ierr) call MPI_Bcast(grid,3,MPI_INTEGER,0,MPI_COMM_WORLD, ierr)
@ -96,11 +96,11 @@ subroutine discretization_grid_init(restart)
call MPI_Bcast(origin,3,MPI_DOUBLE,0,MPI_COMM_WORLD, ierr) call MPI_Bcast(origin,3,MPI_DOUBLE,0,MPI_COMM_WORLD, ierr)
if (ierr /= 0) error stop 'MPI error' if (ierr /= 0) error stop 'MPI error'
print'(/,a,3(i12 ))', ' cells a b c: ', grid print'(/,1x,a,3(i12,1x))', 'cells a b c: ', grid
print'(a,3(es12.5))', ' size x y z: ', geomSize print '(1x,a,3(es12.5,1x))', 'size x y z: ', geomSize
print'(a,3(es12.5))', ' origin x y z: ', origin print '(1x,a,3(es12.5,1x))', 'origin x y z: ', origin
if(worldsize>grid(3)) call IO_error(894, ext_msg='number of processes exceeds grid(3)') if (worldsize>grid(3)) call IO_error(894, ext_msg='number of processes exceeds grid(3)')
call fftw_mpi_init call fftw_mpi_init
devNull = fftw_mpi_local_size_3d(int(grid(3),C_INTPTR_T), & devNull = fftw_mpi_local_size_3d(int(grid(3),C_INTPTR_T), &
@ -109,7 +109,7 @@ subroutine discretization_grid_init(restart)
PETSC_COMM_WORLD, & PETSC_COMM_WORLD, &
z, & ! domain grid size along z z, & ! domain grid size along z
z_offset) ! domain grid offset along z z_offset) ! domain grid offset along z
if(z==0_C_INTPTR_T) call IO_error(894, ext_msg='Cannot distribute MPI processes') if (z==0_C_INTPTR_T) call IO_error(894, ext_msg='Cannot distribute MPI processes')
grid3 = int(z) grid3 = int(z)
grid3Offset = int(z_offset) grid3Offset = int(z_offset)
@ -136,14 +136,14 @@ subroutine discretization_grid_init(restart)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! store geometry information for post processing ! store geometry information for post processing
if(.not. restart) then if (.not. restart) then
call results_openJobFile call results_openJobFile
call results_closeGroup(results_addGroup('geometry')) call results_closeGroup(results_addGroup('geometry'))
call results_addAttribute('cells', grid, '/geometry') call results_addAttribute('cells', grid, '/geometry')
call results_addAttribute('size', geomSize,'/geometry') call results_addAttribute('size', geomSize,'/geometry')
call results_addAttribute('origin',origin, '/geometry') call results_addAttribute('origin',origin, '/geometry')
call results_closeJobFile call results_closeJobFile
endif end if
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! geometry information required by the nonlocal CP model ! geometry information required by the nonlocal CP model
@ -202,18 +202,18 @@ subroutine readVTI(grid,geomSize,origin,material, &
if (endPos < startPos) endPos = len(fileContent,kind=pI64) ! end of file without new line if (endPos < startPos) endPos = len(fileContent,kind=pI64) ! end of file without new line
if (.not. inFile) then if (.not. inFile) then
if(index(fileContent(startPos:endPos),'<VTKFile',kind=pI64) /= 0_pI64) then if (index(fileContent(startPos:endPos),'<VTKFile',kind=pI64) /= 0_pI64) then
inFile = .true. inFile = .true.
if (.not. fileFormatOk(fileContent(startPos:endPos))) call IO_error(error_ID = 844, ext_msg='file format') if (.not. fileFormatOk(fileContent(startPos:endPos))) call IO_error(error_ID = 844, ext_msg='file format')
headerType = merge('UInt64','UInt32',getXMLValue(fileContent(startPos:endPos),'header_type')=='UInt64') headerType = merge('UInt64','UInt32',getXMLValue(fileContent(startPos:endPos),'header_type')=='UInt64')
compressed = getXMLValue(fileContent(startPos:endPos),'compressor') == 'vtkZLibDataCompressor' compressed = getXMLValue(fileContent(startPos:endPos),'compressor') == 'vtkZLibDataCompressor'
endif end if
else else
if (.not. inImage) then if (.not. inImage) then
if (index(fileContent(startPos:endPos),'<ImageData',kind=pI64) /= 0_pI64) then if (index(fileContent(startPos:endPos),'<ImageData',kind=pI64) /= 0_pI64) then
inImage = .true. inImage = .true.
call cellsSizeOrigin(grid,geomSize,origin,fileContent(startPos:endPos)) call cellsSizeOrigin(grid,geomSize,origin,fileContent(startPos:endPos))
endif end if
else else
if (index(fileContent(startPos:endPos),'<CellData>',kind=pI64) /= 0_pI64) then if (index(fileContent(startPos:endPos),'<CellData>',kind=pI64) /= 0_pI64) then
gotCellData = .true. gotCellData = .true.
@ -230,25 +230,25 @@ subroutine readVTI(grid,geomSize,origin,material, &
s = startPos + verify(fileContent(startPos:endPos),IO_WHITESPACE,kind=pI64) -1_pI64 ! start (no leading whitespace) s = startPos + verify(fileContent(startPos:endPos),IO_WHITESPACE,kind=pI64) -1_pI64 ! start (no leading whitespace)
material = as_Int(fileContent(s:endPos),headerType,compressed,dataType) material = as_Int(fileContent(s:endPos),headerType,compressed,dataType)
exit exit
endif end if
startPos = endPos + 2_pI64 startPos = endPos + 2_pI64
endPos = startPos + index(fileContent(startPos:),IO_EOL,kind=pI64) - 2_pI64 endPos = startPos + index(fileContent(startPos:),IO_EOL,kind=pI64) - 2_pI64
enddo end do
endif end if
endif end if
endif end if
if (gotCellData) exit if (gotCellData) exit
startPos = endPos + 2_pI64 startPos = endPos + 2_pI64
enddo end do
if(.not. allocated(material)) call IO_error(error_ID = 844, ext_msg='material data not found') if (.not. allocated(material)) call IO_error(error_ID = 844, ext_msg='material data not found')
if(size(material) /= product(grid)) call IO_error(error_ID = 844, ext_msg='size(material)') if (size(material) /= product(grid)) call IO_error(error_ID = 844, ext_msg='size(material)')
if(any(geomSize<=0)) call IO_error(error_ID = 844, ext_msg='size') if (any(geomSize<=0)) call IO_error(error_ID = 844, ext_msg='size')
if(any(grid<1)) call IO_error(error_ID = 844, ext_msg='grid') if (any(grid<1)) call IO_error(error_ID = 844, ext_msg='grid')
material = material + 1 material = material + 1
if(any(material<1)) call IO_error(error_ID = 844, ext_msg='material ID < 0') if (any(material<1)) call IO_error(error_ID = 844, ext_msg='material ID < 0')
contains contains
@ -352,11 +352,11 @@ subroutine readVTI(grid,geomSize,origin,material, &
integer(C_SIGNED_CHAR), dimension(:), allocatable :: bytes integer(C_SIGNED_CHAR), dimension(:), allocatable :: bytes
if(compressed) then if (compressed) then
bytes = asBytes_compressed(base64_str,headerType) bytes = asBytes_compressed(base64_str,headerType)
else else
bytes = asBytes_uncompressed(base64_str,headerType) bytes = asBytes_uncompressed(base64_str,headerType)
endif end if
end function asBytes end function asBytes
@ -384,12 +384,12 @@ subroutine readVTI(grid,geomSize,origin,material, &
nBlock = int(temp(1),pI64) nBlock = int(temp(1),pI64)
headerLen = 4_pI64 * (3_pI64 + nBlock) headerLen = 4_pI64 * (3_pI64 + nBlock)
temp = int(prec_bytesToC_INT32_T(base64_to_bytes(base64_str(:base64_nChar(headerLen)))),pI64) temp = int(prec_bytesToC_INT32_T(base64_to_bytes(base64_str(:base64_nChar(headerLen)))),pI64)
elseif(headerType == 'UInt64') then elseif (headerType == 'UInt64') then
temp = int(prec_bytesToC_INT64_T(base64_to_bytes(base64_str(:base64_nChar(8_pI64)))),pI64) temp = int(prec_bytesToC_INT64_T(base64_to_bytes(base64_str(:base64_nChar(8_pI64)))),pI64)
nBlock = int(temp(1),pI64) nBlock = int(temp(1),pI64)
headerLen = 8_pI64 * (3_pI64 + nBlock) headerLen = 8_pI64 * (3_pI64 + nBlock)
temp = int(prec_bytesToC_INT64_T(base64_to_bytes(base64_str(:base64_nChar(headerLen)))),pI64) temp = int(prec_bytesToC_INT64_T(base64_to_bytes(base64_str(:base64_nChar(headerLen)))),pI64)
endif end if
allocate(size_inflated(nBlock),source=temp(2)) allocate(size_inflated(nBlock),source=temp(2))
size_inflated(nBlock) = merge(temp(3),temp(2),temp(3)/=0_pI64) size_inflated(nBlock) = merge(temp(3),temp(2),temp(3)/=0_pI64)
@ -402,7 +402,7 @@ subroutine readVTI(grid,geomSize,origin,material, &
s = e + 1_pI64 s = e + 1_pI64
e = s + size_deflated(b) - 1_pI64 e = s + size_deflated(b) - 1_pI64
bytes(sum(size_inflated(:b-1))+1_pI64:sum(size_inflated(:b))) = zlib_inflate(bytes_inflated(s:e),size_inflated(b)) bytes(sum(size_inflated(:b-1))+1_pI64:sum(size_inflated(:b))) = zlib_inflate(bytes_inflated(s:e),size_inflated(b))
enddo end do
end function asBytes_compressed end function asBytes_compressed
@ -429,14 +429,14 @@ subroutine readVTI(grid,geomSize,origin,material, &
nByte = int(prec_bytesToC_INT32_T(base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(4_pI64)))),pI64) nByte = int(prec_bytesToC_INT32_T(base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(4_pI64)))),pI64)
bytes = [bytes,base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(4_pI64+nByte(1))),5_pI64)] bytes = [bytes,base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(4_pI64+nByte(1))),5_pI64)]
s = s + base64_nChar(4_pI64+nByte(1)) s = s + base64_nChar(4_pI64+nByte(1))
enddo end do
elseif(headerType == 'UInt64') then elseif (headerType == 'UInt64') then
do while(s+base64_nChar(8_pI64)<(len(base64_str,pI64))) do while(s+base64_nChar(8_pI64)<(len(base64_str,pI64)))
nByte = int(prec_bytesToC_INT64_T(base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(8_pI64)))),pI64) nByte = int(prec_bytesToC_INT64_T(base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(8_pI64)))),pI64)
bytes = [bytes,base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(8_pI64+nByte(1))),9_pI64)] bytes = [bytes,base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(8_pI64+nByte(1))),9_pI64)]
s = s + base64_nChar(8_pI64+nByte(1)) s = s + base64_nChar(8_pI64+nByte(1))
enddo end do
endif end if
end function asBytes_uncompressed end function asBytes_uncompressed
@ -455,11 +455,11 @@ subroutine readVTI(grid,geomSize,origin,material, &
#endif #endif
s = index(line," "//key,back=.true.) s = index(line," "//key,back=.true.)
if(s==0) then if (s==0) then
getXMLValue = '' getXMLValue = ''
else else
e = s + 1 + scan(line(s+1:),"'"//'"') e = s + 1 + scan(line(s+1:),"'"//'"')
if(scan(line(s:e-2),'=') == 0) then if (scan(line(s:e-2),'=') == 0) then
getXMLValue = '' getXMLValue = ''
else else
s = e s = e
@ -471,8 +471,8 @@ subroutine readVTI(grid,geomSize,origin,material, &
e = s + index(line(s:),merge("'",'"',line(s-1:s-1)=="'")) - 1 e = s + index(line(s:),merge("'",'"',line(s-1:s-1)=="'")) - 1
#endif #endif
getXMLValue = line(s:e-1) getXMLValue = line(s:e-1)
endif end if
endif end if
end function end function
@ -514,7 +514,7 @@ function IPcoordinates0(grid,geomSize,grid3Offset)
do c = 1, grid(3); do b = 1, grid(2); do a = 1, grid(1) do c = 1, grid(3); do b = 1, grid(2); do a = 1, grid(1)
i = i + 1 i = i + 1
IPcoordinates0(1:3,i) = geomSize/real(grid,pReal) * (real([a,b,grid3Offset+c],pReal) -0.5_pReal) IPcoordinates0(1:3,i) = geomSize/real(grid,pReal) * (real([a,b,grid3Offset+c],pReal) -0.5_pReal)
enddo; enddo; enddo end do; end do; end do
end function IPcoordinates0 end function IPcoordinates0
@ -538,7 +538,7 @@ pure function nodes0(grid,geomSize,grid3Offset)
do c = 0, grid3; do b = 0, grid(2); do a = 0, grid(1) do c = 0, grid3; do b = 0, grid(2); do a = 0, grid(1)
n = n + 1 n = n + 1
nodes0(1:3,n) = geomSize/real(grid,pReal) * real([a,b,grid3Offset+c],pReal) nodes0(1:3,n) = geomSize/real(grid,pReal) * real([a,b,grid3Offset+c],pReal)
enddo; enddo; enddo end do; end do; end do
end function nodes0 end function nodes0
@ -631,7 +631,7 @@ pure function IPneighborhood(grid)
IPneighborhood(3,5,1,e) = 6 IPneighborhood(3,5,1,e) = 6
IPneighborhood(3,6,1,e) = 5 IPneighborhood(3,6,1,e) = 5
enddo; enddo; enddo end do; end do; end do
end function IPneighborhood end function IPneighborhood

28
src/grid/grid_damage_spectral.f90
View File

@ -76,10 +76,10 @@ subroutine grid_damage_spectral_init()
character(len=pStringLen) :: & character(len=pStringLen) :: &
snes_type snes_type
print'(/,a)', ' <<<+- grid_spectral_damage init -+>>>' print'(/,1x,a)', '<<<+- grid_spectral_damage init -+>>>'
print*, 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019' print'(/,1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019'
print*, 'https://doi.org/10.1007/978-981-10-6855-3_80' print'( 1x,a)', 'https://doi.org/10.1007/978-981-10-6855-3_80'
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! read numerical parameters and do sanity checks ! read numerical parameters and do sanity checks
@ -137,7 +137,7 @@ subroutine grid_damage_spectral_init()
call SNESVISetVariableBounds(damage_snes,lBound,uBound,ierr) ! variable bounds for variational inequalities like contact mechanics, damage etc. call SNESVISetVariableBounds(damage_snes,lBound,uBound,ierr) ! variable bounds for variational inequalities like contact mechanics, damage etc.
call DMRestoreGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr) call DMRestoreGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr)
call DMRestoreGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr) call DMRestoreGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr)
endif end if
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init fields ! init fields
@ -187,7 +187,7 @@ function grid_damage_spectral_solution(Delta_t) result(solution)
else else
solution%converged = .true. solution%converged = .true.
solution%iterationsNeeded = totalIter solution%iterationsNeeded = totalIter
endif end if
stagNorm = maxval(abs(phi_current - phi_stagInc)) stagNorm = maxval(abs(phi_current - phi_stagInc))
solnNorm = maxval(abs(phi_current)) solnNorm = maxval(abs(phi_current))
call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD,ierr)
@ -201,14 +201,14 @@ function grid_damage_spectral_solution(Delta_t) result(solution)
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1 ce = ce + 1
call homogenization_set_phi(phi_current(i,j,k),ce) call homogenization_set_phi(phi_current(i,j,k),ce)
enddo; enddo; enddo end do; end do; end do
call VecMin(solution_vec,devNull,phi_min,ierr); CHKERRQ(ierr) call VecMin(solution_vec,devNull,phi_min,ierr); CHKERRQ(ierr)
call VecMax(solution_vec,devNull,phi_max,ierr); CHKERRQ(ierr) call VecMax(solution_vec,devNull,phi_max,ierr); CHKERRQ(ierr)
if (solution%converged) & if (solution%converged) &
print'(/,a)', ' ... nonlocal damage converged .....................................' print'(/,1x,a)', '... nonlocal damage converged .....................................'
print'(/,a,f8.6,2x,f8.6,2x,e11.4)', ' Minimum|Maximum|Delta Damage = ', phi_min, phi_max, stagNorm print'(/,1x,a,f8.6,2x,f8.6,2x,e11.4)', 'Minimum|Maximum|Delta Damage = ', phi_min, phi_max, stagNorm
print'(/,a)', ' ===========================================================================' print'(/,1x,a)', '==========================================================================='
flush(IO_STDOUT) flush(IO_STDOUT)
end function grid_damage_spectral_solution end function grid_damage_spectral_solution
@ -238,11 +238,11 @@ subroutine grid_damage_spectral_forward(cutBack)
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1 ce = ce + 1
call homogenization_set_phi(phi_current(i,j,k),ce) call homogenization_set_phi(phi_current(i,j,k),ce)
enddo; enddo; enddo end do; end do; end do
else else
phi_lastInc = phi_current phi_lastInc = phi_current
call updateReference call updateReference
endif end if
end subroutine grid_damage_spectral_forward end subroutine grid_damage_spectral_forward
@ -277,7 +277,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1 ce = ce + 1
vectorField_real(1:3,i,j,k) = matmul(homogenization_K_phi(ce) - K_ref, vectorField_real(1:3,i,j,k)) vectorField_real(1:3,i,j,k) = matmul(homogenization_K_phi(ce) - K_ref, vectorField_real(1:3,i,j,k))
enddo; enddo; enddo end do; end do; end do
call utilities_FFTvectorForward call utilities_FFTvectorForward
call utilities_fourierVectorDivergence !< calculate damage divergence in fourier field call utilities_fourierVectorDivergence !< calculate damage divergence in fourier field
call utilities_FFTscalarBackward call utilities_FFTscalarBackward
@ -287,7 +287,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
scalarField_real(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_phi(phi_current(i,j,k),ce)) & scalarField_real(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_phi(phi_current(i,j,k),ce)) &
+ homogenization_mu_phi(ce)*(phi_lastInc(i,j,k) - phi_current(i,j,k)) & + homogenization_mu_phi(ce)*(phi_lastInc(i,j,k) - phi_current(i,j,k)) &
+ mu_ref*phi_current(i,j,k) + mu_ref*phi_current(i,j,k)
enddo; enddo; enddo end do; end do; end do
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! convolution of damage field with green operator ! convolution of damage field with green operator
@ -320,7 +320,7 @@ subroutine updateReference()
do ce = 1, product(grid(1:2))*grid3 do ce = 1, product(grid(1:2))*grid3
K_ref = K_ref + homogenization_K_phi(ce) K_ref = K_ref + homogenization_K_phi(ce)
mu_ref = mu_ref + homogenization_mu_phi(ce) mu_ref = mu_ref + homogenization_mu_phi(ce)
enddo end do
K_ref = K_ref*wgt K_ref = K_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr)

24
src/grid/grid_mech_FEM.f90
View File

@ -116,7 +116,7 @@ subroutine grid_mechanical_FEM_init
num_grid, & num_grid, &
debug_grid debug_grid
print'(/,a)', ' <<<+- grid_mechanical_FEM init -+>>>'; flush(IO_STDOUT) print'(/,1x,a)', '<<<+- grid_mechanical_FEM init -+>>>'; flush(IO_STDOUT)
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! debugging options ! debugging options
@ -234,7 +234,7 @@ subroutine grid_mechanical_FEM_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init fields ! init fields
restartRead: if (interface_restartInc > 0) then restartRead: if (interface_restartInc > 0) then
print'(/,a,i0,a)', ' reading restart data of increment ', interface_restartInc, ' from file' print'(/,1x,a,i0,a)', 'reading restart data of increment ', interface_restartInc, ' from file'
fileHandle = HDF5_openFile(getSolverJobName()//'_restart.hdf5','r') fileHandle = HDF5_openFile(getSolverJobName()//'_restart.hdf5','r')
groupHandle = HDF5_openGroup(fileHandle,'solver') groupHandle = HDF5_openGroup(fileHandle,'solver')
@ -272,7 +272,7 @@ subroutine grid_mechanical_FEM_init
CHKERRQ(ierr) CHKERRQ(ierr)
restartRead2: if (interface_restartInc > 0) then restartRead2: if (interface_restartInc > 0) then
print'(a,i0,a)', ' reading more restart data of increment ', interface_restartInc, ' from file' print'(1x,a,i0,a)', 'reading more restart data of increment ', interface_restartInc, ' from file'
call HDF5_read(C_volAvg,groupHandle,'C_volAvg',.false.) call HDF5_read(C_volAvg,groupHandle,'C_volAvg',.false.)
call MPI_Bcast(C_volAvg,81,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr) call MPI_Bcast(C_volAvg,81,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
if(ierr /=0) error stop 'MPI error' if(ierr /=0) error stop 'MPI error'
@ -442,7 +442,7 @@ subroutine grid_mechanical_FEM_restartWrite
call DMDAVecGetArrayF90(mechanical_grid,solution_lastInc,u_lastInc,ierr) call DMDAVecGetArrayF90(mechanical_grid,solution_lastInc,u_lastInc,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
print*, 'writing solver data required for restart to file'; flush(IO_STDOUT) print'(1x,a)', 'writing solver data required for restart to file'; flush(IO_STDOUT)
fileHandle = HDF5_openFile(getSolverJobName()//'_restart.hdf5','w') fileHandle = HDF5_openFile(getSolverJobName()//'_restart.hdf5','w')
groupHandle = HDF5_addGroup(fileHandle,'solver') groupHandle = HDF5_addGroup(fileHandle,'solver')
@ -506,12 +506,12 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,i
reason = 0 reason = 0
endif endif
print'(1/,a)', ' ... reporting .............................................................' print'(/,1x,a)', '... reporting .............................................................'
print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', & print'(/,1x,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,')'
print'(a,f12.2,a,es8.2,a,es9.2,a)', ' error stress BC = ', & print'(1x,a,f12.2,a,es8.2,a,es9.2,a)', 'error stress BC = ', &
err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')'
print'(/,a)', ' ===========================================================================' print'(/,1x,a)', '==========================================================================='
flush(IO_STDOUT) flush(IO_STDOUT)
end subroutine converged end subroutine converged
@ -547,10 +547,10 @@ subroutine formResidual(da_local,x_local, &
newIteration: if (totalIter <= PETScIter) then newIteration: if (totalIter <= PETScIter) then
totalIter = totalIter + 1 totalIter = totalIter + 1
print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter+1, '≤', num%itmax print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter+1, '≤', num%itmax
if (debugRotation) print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', & if (debugRotation) print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) 'deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', & print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
' deformation gradient aim =', transpose(F_aim) 'deformation gradient aim =', transpose(F_aim)
flush(IO_STDOUT) flush(IO_STDOUT)
endif newIteration endif newIteration

60
src/grid/grid_mech_spectral_basic.f90
View File

@ -111,13 +111,13 @@ subroutine grid_mechanical_spectral_basic_init
num_grid, & num_grid, &
debug_grid debug_grid
print'(/,a)', ' <<<+- grid_mechanical_spectral_basic init -+>>>'; flush(IO_STDOUT) print'(/,1x,a)', '<<<+- grid_mechanical_spectral_basic init -+>>>'; flush(IO_STDOUT)
print*, 'P. Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013' print'(/,1x,a)', 'P. Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013'
print*, 'https://doi.org/10.1016/j.ijplas.2012.09.012'//IO_EOL print'( 1x,a)', 'https://doi.org/10.1016/j.ijplas.2012.09.012'//IO_EOL
print*, 'P. Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' print'( 1x,a)', 'P. Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006' print'( 1x,a)', 'https://doi.org/10.1016/j.ijplas.2014.02.006'
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! debugging options ! debugging options
@ -186,30 +186,30 @@ subroutine grid_mechanical_spectral_basic_init
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! places pointer on PETSc data call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! places pointer on PETSc data
restartRead: if (interface_restartInc > 0) then restartRead: if (interface_restartInc > 0) then
print'(/,a,i0,a)', ' reading restart data of increment ', interface_restartInc, ' from file' print'(/,1x,a,i0,a)', 'reading restart data of increment ', interface_restartInc, ' from file'
fileHandle = HDF5_openFile(getSolverJobName()//'_restart.hdf5','r') fileHandle = HDF5_openFile(getSolverJobName()//'_restart.hdf5','r')
groupHandle = HDF5_openGroup(fileHandle,'solver') groupHandle = HDF5_openGroup(fileHandle,'solver')
call HDF5_read(P_aim,groupHandle,'P_aim',.false.) call HDF5_read(P_aim,groupHandle,'P_aim',.false.)
call MPI_Bcast(P_aim,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr) call MPI_Bcast(P_aim,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
if(ierr /=0) error stop 'MPI error' if (ierr /=0) error stop 'MPI error'
call HDF5_read(F_aim,groupHandle,'F_aim',.false.) call HDF5_read(F_aim,groupHandle,'F_aim',.false.)
call MPI_Bcast(F_aim,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr) call MPI_Bcast(F_aim,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
if(ierr /=0) error stop 'MPI error' if (ierr /=0) error stop 'MPI error'
call HDF5_read(F_aim_lastInc,groupHandle,'F_aim_lastInc',.false.) call HDF5_read(F_aim_lastInc,groupHandle,'F_aim_lastInc',.false.)
call MPI_Bcast(F_aim_lastInc,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr) call MPI_Bcast(F_aim_lastInc,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
if(ierr /=0) error stop 'MPI error' if (ierr /=0) error stop 'MPI error'
call HDF5_read(F_aimDot,groupHandle,'F_aimDot',.false.) call HDF5_read(F_aimDot,groupHandle,'F_aimDot',.false.)
call MPI_Bcast(F_aimDot,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr) call MPI_Bcast(F_aimDot,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
if(ierr /=0) error stop 'MPI error' if (ierr /=0) error stop 'MPI error'
call HDF5_read(F,groupHandle,'F') call HDF5_read(F,groupHandle,'F')
call HDF5_read(F_lastInc,groupHandle,'F_lastInc') call HDF5_read(F_lastInc,groupHandle,'F_lastInc')
elseif (interface_restartInc == 0) then restartRead elseif (interface_restartInc == 0) then restartRead
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 restartRead end if restartRead
homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for homogenization_mechanical_response homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for homogenization_mechanical_response
call utilities_updateCoords(reshape(F,shape(F_lastInc))) call utilities_updateCoords(reshape(F,shape(F_lastInc)))
@ -219,13 +219,13 @@ subroutine grid_mechanical_spectral_basic_init
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! deassociate pointer call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! deassociate pointer
restartRead2: if (interface_restartInc > 0) then restartRead2: if (interface_restartInc > 0) then
print'(a,i0,a)', ' reading more restart data of increment ', interface_restartInc, ' from file' print'(1x,a,i0,a)', 'reading more restart data of increment ', interface_restartInc, ' from file'
call HDF5_read(C_volAvg,groupHandle,'C_volAvg',.false.) call HDF5_read(C_volAvg,groupHandle,'C_volAvg',.false.)
call MPI_Bcast(C_volAvg,81,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr) call MPI_Bcast(C_volAvg,81,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
if(ierr /=0) error stop 'MPI error' if (ierr /=0) error stop 'MPI error'
call HDF5_read(C_volAvgLastInc,groupHandle,'C_volAvgLastInc',.false.) call HDF5_read(C_volAvgLastInc,groupHandle,'C_volAvgLastInc',.false.)
call MPI_Bcast(C_volAvgLastInc,81,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr) call MPI_Bcast(C_volAvgLastInc,81,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
if(ierr /=0) error stop 'MPI error' if (ierr /=0) error stop 'MPI error'
call HDF5_closeGroup(groupHandle) call HDF5_closeGroup(groupHandle)
call HDF5_closeFile(fileHandle) call HDF5_closeFile(fileHandle)
@ -234,7 +234,7 @@ subroutine grid_mechanical_spectral_basic_init
MPI_MODE_RDONLY,MPI_INFO_NULL,fileUnit,ierr) MPI_MODE_RDONLY,MPI_INFO_NULL,fileUnit,ierr)
call MPI_File_read(fileUnit,C_minMaxAvg,81,MPI_DOUBLE,MPI_STATUS_IGNORE,ierr) call MPI_File_read(fileUnit,C_minMaxAvg,81,MPI_DOUBLE,MPI_STATUS_IGNORE,ierr)
call MPI_File_close(fileUnit,ierr) call MPI_File_close(fileUnit,ierr)
endif restartRead2 end if restartRead2
call utilities_updateGamma(C_minMaxAvg) call utilities_updateGamma(C_minMaxAvg)
call utilities_saveReferenceStiffness call utilities_saveReferenceStiffness
@ -263,7 +263,7 @@ function grid_mechanical_spectral_basic_solution(incInfoIn) result(solution)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator) ! update stiffness (and gamma operator)
S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg) S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg)
if(num%update_gamma) call utilities_updateGamma(C_minMaxAvg) if (num%update_gamma) call utilities_updateGamma(C_minMaxAvg)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! solve BVP ! solve BVP
@ -328,7 +328,7 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_
elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed
F_aimDot = F_aimDot & F_aimDot = F_aimDot &
+ merge(.0_pReal,(deformation_BC%values - F_aim_lastInc)/t_remaining,deformation_BC%mask) + merge(.0_pReal,(deformation_BC%values - F_aim_lastInc)/t_remaining,deformation_BC%mask)
endif end if
Fdot = utilities_calculateRate(guess, & Fdot = utilities_calculateRate(guess, &
F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),Delta_t_old, & F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),Delta_t_old, &
@ -336,7 +336,7 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_
F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid3]) F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid3])
homogenization_F0 = reshape(F,[3,3,product(grid(1:2))*grid3]) homogenization_F0 = reshape(F,[3,3,product(grid(1:2))*grid3])
endif end if
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update average and local deformation gradients ! update average and local deformation gradients
@ -385,7 +385,7 @@ subroutine grid_mechanical_spectral_basic_restartWrite
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
print*, 'writing solver data required for restart to file'; flush(IO_STDOUT) print'(1x,a)', 'writing solver data required for restart to file'; flush(IO_STDOUT)
fileHandle = HDF5_openFile(getSolverJobName()//'_restart.hdf5','w') fileHandle = HDF5_openFile(getSolverJobName()//'_restart.hdf5','w')
groupHandle = HDF5_addGroup(fileHandle,'solver') groupHandle = HDF5_addGroup(fileHandle,'solver')
@ -406,7 +406,7 @@ subroutine grid_mechanical_spectral_basic_restartWrite
call HDF5_write(C_minMaxAvg,groupHandle,'C_minMaxAvg',.false.) call HDF5_write(C_minMaxAvg,groupHandle,'C_minMaxAvg',.false.)
call HDF5_closeGroup(groupHandle) call HDF5_closeGroup(groupHandle)
call HDF5_closeFile(fileHandle) call HDF5_closeFile(fileHandle)
endif end if
if (num%update_gamma) call utilities_saveReferenceStiffness if (num%update_gamma) call utilities_saveReferenceStiffness
@ -443,14 +443,14 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
reason = -1 reason = -1
else else
reason = 0 reason = 0
endif end if
print'(1/,a)', ' ... reporting .............................................................' print'(/,1x,a)', '... reporting .............................................................'
print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', & print'(/,1x,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,')'
print'(a,f12.2,a,es8.2,a,es9.2,a)', ' error stress BC = ', & print'(1x,a,f12.2,a,es8.2,a,es9.2,a)', 'error stress BC = ', &
err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')'
print'(/,a)', ' ===========================================================================' print'(/,1x,a)', '==========================================================================='
flush(IO_STDOUT) flush(IO_STDOUT)
end subroutine converged end subroutine converged
@ -485,12 +485,12 @@ subroutine formResidual(in, F, &
newIteration: if (totalIter <= PETScIter) then newIteration: if (totalIter <= PETScIter) then
totalIter = totalIter + 1 totalIter = totalIter + 1
print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax
if (debugRotation) print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', & if (debugRotation) print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) 'deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', & print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
' deformation gradient aim =', transpose(F_aim) 'deformation gradient aim =', transpose(F_aim)
flush(IO_STDOUT) flush(IO_STDOUT)
endif newIteration end if newIteration
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate constitutive response ! evaluate constitutive response

74
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -124,10 +124,10 @@ subroutine grid_mechanical_spectral_polarisation_init
num_grid, & num_grid, &
debug_grid debug_grid
print'(/,a)', ' <<<+- grid_mechanical_spectral_polarization init -+>>>'; flush(IO_STDOUT) print'(/,1x,a)', '<<<+- grid_mechanical_spectral_polarization init -+>>>'; flush(IO_STDOUT)
print*, 'P. Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' print'(/,1x,a)', 'P. Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006' print'( 1x,a)', 'https://doi.org/10.1016/j.ijplas.2014.02.006'
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! debugging options ! debugging options
@ -208,23 +208,23 @@ subroutine grid_mechanical_spectral_polarisation_init
F_tau => FandF_tau(9:17,:,:,:) F_tau => FandF_tau(9:17,:,:,:)
restartRead: if (interface_restartInc > 0) then restartRead: if (interface_restartInc > 0) then
print'(/,a,i0,a)', ' reading restart data of increment ', interface_restartInc, ' from file' print'(/,1x,a,i0,a)', 'reading restart data of increment ', interface_restartInc, ' from file'
fileHandle = HDF5_openFile(getSolverJobName()//'_restart.hdf5','r') fileHandle = HDF5_openFile(getSolverJobName()//'_restart.hdf5','r')
groupHandle = HDF5_openGroup(fileHandle,'solver') groupHandle = HDF5_openGroup(fileHandle,'solver')
call HDF5_read(P_aim,groupHandle,'P_aim',.false.) call HDF5_read(P_aim,groupHandle,'P_aim',.false.)
call MPI_Bcast(P_aim,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr) call MPI_Bcast(P_aim,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
if(ierr /=0) error stop 'MPI error' if (ierr /=0) error stop 'MPI error'
call HDF5_read(F_aim,groupHandle,'F_aim',.false.) call HDF5_read(F_aim,groupHandle,'F_aim',.false.)
call MPI_Bcast(F_aim,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr) call MPI_Bcast(F_aim,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
if(ierr /=0) error stop 'MPI error' if (ierr /=0) error stop 'MPI error'
call HDF5_read(F_aim_lastInc,groupHandle,'F_aim_lastInc',.false.) call HDF5_read(F_aim_lastInc,groupHandle,'F_aim_lastInc',.false.)
call MPI_Bcast(F_aim_lastInc,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr) call MPI_Bcast(F_aim_lastInc,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
if(ierr /=0) error stop 'MPI error' if (ierr /=0) error stop 'MPI error'
call HDF5_read(F_aimDot,groupHandle,'F_aimDot',.false.) call HDF5_read(F_aimDot,groupHandle,'F_aimDot',.false.)
call MPI_Bcast(F_aimDot,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr) call MPI_Bcast(F_aimDot,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
if(ierr /=0) error stop 'MPI error' if (ierr /=0) error stop 'MPI error'
call HDF5_read(F,groupHandle,'F') call HDF5_read(F,groupHandle,'F')
call HDF5_read(F_lastInc,groupHandle,'F_lastInc') call HDF5_read(F_lastInc,groupHandle,'F_lastInc')
call HDF5_read(F_tau,groupHandle,'F_tau') call HDF5_read(F_tau,groupHandle,'F_tau')
@ -235,7 +235,7 @@ subroutine grid_mechanical_spectral_polarisation_init
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 restartRead end if restartRead
homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for homogenization_mechanical_response homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for homogenization_mechanical_response
call utilities_updateCoords(reshape(F,shape(F_lastInc))) call utilities_updateCoords(reshape(F,shape(F_lastInc)))
@ -245,13 +245,13 @@ subroutine grid_mechanical_spectral_polarisation_init
call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! deassociate pointer call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! deassociate pointer
restartRead2: if (interface_restartInc > 0) then restartRead2: if (interface_restartInc > 0) then
print'(a,i0,a)', ' reading more restart data of increment ', interface_restartInc, ' from file' print'(1x,a,i0,a)', 'reading more restart data of increment ', interface_restartInc, ' from file'
call HDF5_read(C_volAvg,groupHandle,'C_volAvg',.false.) call HDF5_read(C_volAvg,groupHandle,'C_volAvg',.false.)
call MPI_Bcast(C_volAvg,81,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr) call MPI_Bcast(C_volAvg,81,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
if(ierr /=0) error stop 'MPI error' if (ierr /=0) error stop 'MPI error'
call HDF5_read(C_volAvgLastInc,groupHandle,'C_volAvgLastInc',.false.) call HDF5_read(C_volAvgLastInc,groupHandle,'C_volAvgLastInc',.false.)
call MPI_Bcast(C_volAvgLastInc,81,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr) call MPI_Bcast(C_volAvgLastInc,81,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
if(ierr /=0) error stop 'MPI error' if (ierr /=0) error stop 'MPI error'
call HDF5_closeGroup(groupHandle) call HDF5_closeGroup(groupHandle)
call HDF5_closeFile(fileHandle) call HDF5_closeFile(fileHandle)
@ -260,7 +260,7 @@ subroutine grid_mechanical_spectral_polarisation_init
MPI_MODE_RDONLY,MPI_INFO_NULL,fileUnit,ierr) MPI_MODE_RDONLY,MPI_INFO_NULL,fileUnit,ierr)
call MPI_File_read(fileUnit,C_minMaxAvg,81,MPI_DOUBLE,MPI_STATUS_IGNORE,ierr) call MPI_File_read(fileUnit,C_minMaxAvg,81,MPI_DOUBLE,MPI_STATUS_IGNORE,ierr)
call MPI_File_close(fileUnit,ierr) call MPI_File_close(fileUnit,ierr)
endif restartRead2 end if restartRead2
call utilities_updateGamma(C_minMaxAvg) call utilities_updateGamma(C_minMaxAvg)
call utilities_saveReferenceStiffness call utilities_saveReferenceStiffness
@ -291,11 +291,11 @@ function grid_mechanical_spectral_polarisation_solution(incInfoIn) result(soluti
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator) ! update stiffness (and gamma operator)
S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg) S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg)
if(num%update_gamma) then if (num%update_gamma) then
call utilities_updateGamma(C_minMaxAvg) call utilities_updateGamma(C_minMaxAvg)
C_scale = C_minMaxAvg C_scale = C_minMaxAvg
S_scale = math_invSym3333(C_minMaxAvg) S_scale = math_invSym3333(C_minMaxAvg)
endif end if
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! solve BVP ! solve BVP
@ -364,7 +364,7 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D
elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed
F_aimDot = F_aimDot & F_aimDot = F_aimDot &
+ merge(.0_pReal,(deformation_BC%values - F_aim_lastInc)/t_remaining,deformation_BC%mask) + merge(.0_pReal,(deformation_BC%values - F_aim_lastInc)/t_remaining,deformation_BC%mask)
endif end if
Fdot = utilities_calculateRate(guess, & Fdot = utilities_calculateRate(guess, &
F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),Delta_t_old, & F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),Delta_t_old, &
@ -376,14 +376,14 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D
F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3]) F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3])
homogenization_F0 = reshape(F,[3,3,product(grid(1:2))*grid3]) homogenization_F0 = reshape(F,[3,3,product(grid(1:2))*grid3])
endif end if
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update average and local deformation gradients ! update average and local deformation gradients
F_aim = F_aim_lastInc + F_aimDot * Delta_t F_aim = F_aim_lastInc + F_aimDot * Delta_t
if(stress_BC%myType=='P') P_aim = P_aim & if (stress_BC%myType=='P') P_aim = P_aim &
+ merge(.0_pReal,(stress_BC%values - P_aim)/t_remaining,stress_BC%mask)*Delta_t + merge(.0_pReal,(stress_BC%values - P_aim)/t_remaining,stress_BC%mask)*Delta_t
if(stress_BC%myType=='dot_P') P_aim = P_aim & if (stress_BC%myType=='dot_P') P_aim = P_aim &
+ merge(.0_pReal,stress_BC%values,stress_BC%mask)*Delta_t + merge(.0_pReal,stress_BC%values,stress_BC%mask)*Delta_t
F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average
@ -399,8 +399,8 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D
+ math_mul3333xx33(S_scale,0.5_pReal*matmul(F_lambda33, & + math_mul3333xx33(S_scale,0.5_pReal*matmul(F_lambda33, &
math_mul3333xx33(C_scale,matmul(transpose(F_lambda33),F_lambda33)-math_I3))) math_mul3333xx33(C_scale,matmul(transpose(F_lambda33),F_lambda33)-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 end do; end do; end do
endif end if
call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
@ -441,7 +441,7 @@ subroutine grid_mechanical_spectral_polarisation_restartWrite
F => FandF_tau(0: 8,:,:,:) F => FandF_tau(0: 8,:,:,:)
F_tau => FandF_tau(9:17,:,:,:) F_tau => FandF_tau(9:17,:,:,:)
print*, 'writing solver data required for restart to file'; flush(IO_STDOUT) print'(1x,a)', 'writing solver data required for restart to file'; flush(IO_STDOUT)
fileHandle = HDF5_openFile(getSolverJobName()//'_restart.hdf5','w') fileHandle = HDF5_openFile(getSolverJobName()//'_restart.hdf5','w')
groupHandle = HDF5_addGroup(fileHandle,'solver') groupHandle = HDF5_addGroup(fileHandle,'solver')
@ -463,9 +463,9 @@ subroutine grid_mechanical_spectral_polarisation_restartWrite
call HDF5_write(C_volAvgLastInc,groupHandle,'C_volAvgLastInc',.false.) call HDF5_write(C_volAvgLastInc,groupHandle,'C_volAvgLastInc',.false.)
call HDF5_closeGroup(groupHandle) call HDF5_closeGroup(groupHandle)
call HDF5_closeFile(fileHandle) call HDF5_closeFile(fileHandle)
endif end if
if(num%update_gamma) call utilities_saveReferenceStiffness if (num%update_gamma) call utilities_saveReferenceStiffness
call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
@ -502,16 +502,16 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
reason = -1 reason = -1
else else
reason = 0 reason = 0
endif end if
print'(1/,a)', ' ... reporting .............................................................' print'(/,1x,a)', '... reporting .............................................................'
print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', & print'(/,1x,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,')'
print '(a,f12.2,a,es8.2,a,es9.2,a)', ' error curl = ', & print '(1x,a,f12.2,a,es8.2,a,es9.2,a)', 'error curl = ', &
err_curl/curlTol,' (',err_curl,' -, tol = ',curlTol,')' err_curl/curlTol,' (',err_curl,' -, tol = ',curlTol,')'
print '(a,f12.2,a,es8.2,a,es9.2,a)', ' error stress BC = ', & print '(1x,a,f12.2,a,es8.2,a,es9.2,a)', 'error stress BC = ', &
err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')'
print'(/,a)', ' ===========================================================================' print'(/,1x,a)', '==========================================================================='
flush(IO_STDOUT) flush(IO_STDOUT)
end subroutine converged end subroutine converged
@ -565,12 +565,12 @@ subroutine formResidual(in, FandF_tau, &
newIteration: if (totalIter <= PETScIter) then newIteration: if (totalIter <= PETScIter) then
totalIter = totalIter + 1 totalIter = totalIter + 1
print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax
if (debugRotation) print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', & if (debugRotation) print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) 'deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', & print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
' deformation gradient aim =', transpose(F_aim) 'deformation gradient aim =', transpose(F_aim)
flush(IO_STDOUT) flush(IO_STDOUT)
endif newIteration end if newIteration
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! !
@ -580,7 +580,7 @@ subroutine formResidual(in, FandF_tau, &
num%beta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& num%beta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -&
num%alpha*matmul(F(1:3,1:3,i,j,k), & num%alpha*matmul(F(1:3,1:3,i,j,k), &
math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3)) math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))
enddo; enddo; enddo end do; end do; end do
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! doing convolution in Fourier space ! doing convolution in Fourier space
@ -621,7 +621,7 @@ subroutine formResidual(in, FandF_tau, &
residual_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), & residual_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), &
math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) & math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) &
+ residual_F_tau(1:3,1:3,i,j,k) + residual_F_tau(1:3,1:3,i,j,k)
enddo; enddo; enddo end do; end do; end do
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculating curl ! calculating curl

28
src/grid/grid_thermal_spectral.f90
View File

@ -75,10 +75,10 @@ subroutine grid_thermal_spectral_init(T_0)
class(tNode), pointer :: & class(tNode), pointer :: &
num_grid num_grid
print'(/,a)', ' <<<+- grid_thermal_spectral init -+>>>' print'(/,1x,a)', '<<<+- grid_thermal_spectral init -+>>>'
print*, 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019' print'(/,1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019'
print*, 'https://doi.org/10.1007/978-981-10-6855-3_80' print'( 1x,a)', 'https://doi.org/10.1007/978-981-10-6855-3_80'
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! read numerical parameters and do sanity checks ! read numerical parameters and do sanity checks
@ -141,7 +141,7 @@ subroutine grid_thermal_spectral_init(T_0)
T_lastInc(i,j,k) = T_current(i,j,k) T_lastInc(i,j,k) = T_current(i,j,k)
T_stagInc(i,j,k) = T_current(i,j,k) T_stagInc(i,j,k) = T_current(i,j,k)
call homogenization_thermal_setField(T_0,0.0_pReal,ce) call homogenization_thermal_setField(T_0,0.0_pReal,ce)
enddo; enddo; enddo end do; end do; end do
call DMDAVecGetArrayF90(thermal_grid,solution_vec,T_PETSc,ierr); CHKERRQ(ierr) call DMDAVecGetArrayF90(thermal_grid,solution_vec,T_PETSc,ierr); CHKERRQ(ierr)
T_PETSc(xstart:xend,ystart:yend,zstart:zend) = T_current T_PETSc(xstart:xend,ystart:yend,zstart:zend) = T_current
@ -182,7 +182,7 @@ function grid_thermal_spectral_solution(Delta_t) result(solution)
else else
solution%converged = .true. solution%converged = .true.
solution%iterationsNeeded = totalIter solution%iterationsNeeded = totalIter
endif end if
stagNorm = maxval(abs(T_current - T_stagInc)) stagNorm = maxval(abs(T_current - T_stagInc))
solnNorm = maxval(abs(T_current)) solnNorm = maxval(abs(T_current))
call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD,ierr)
@ -196,14 +196,14 @@ function grid_thermal_spectral_solution(Delta_t) result(solution)
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1 ce = ce + 1
call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%Delta_t,ce) call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%Delta_t,ce)
enddo; enddo; enddo end do; end do; end do
call VecMin(solution_vec,devNull,T_min,ierr); CHKERRQ(ierr) call VecMin(solution_vec,devNull,T_min,ierr); CHKERRQ(ierr)
call VecMax(solution_vec,devNull,T_max,ierr); CHKERRQ(ierr) call VecMax(solution_vec,devNull,T_max,ierr); CHKERRQ(ierr)
if (solution%converged) & if (solution%converged) &
print'(/,a)', ' ... thermal conduction converged ..................................' print'(/,1x,a)', '... thermal conduction converged ..................................'
print'(/,a,f8.4,2x,f8.4,2x,f8.4)', ' Minimum|Maximum|Delta Temperature / K = ', T_min, T_max, stagNorm print'(/,1x,a,f8.4,2x,f8.4,2x,f8.4)', 'Minimum|Maximum|Delta Temperature / K = ', T_min, T_max, stagNorm
print'(/,a)', ' ===========================================================================' print'(/,1x,a)', '==========================================================================='
flush(IO_STDOUT) flush(IO_STDOUT)
end function grid_thermal_spectral_solution end function grid_thermal_spectral_solution
@ -234,11 +234,11 @@ subroutine grid_thermal_spectral_forward(cutBack)
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1 ce = ce + 1
call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%Delta_t,ce) call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%Delta_t,ce)
enddo; enddo; enddo end do; end do; end do
else else
T_lastInc = T_current T_lastInc = T_current
call updateReference call updateReference
endif end if
end subroutine grid_thermal_spectral_forward end subroutine grid_thermal_spectral_forward
@ -272,7 +272,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1 ce = ce + 1
vectorField_real(1:3,i,j,k) = matmul(homogenization_K_T(ce) - K_ref, vectorField_real(1:3,i,j,k)) vectorField_real(1:3,i,j,k) = matmul(homogenization_K_T(ce) - K_ref, vectorField_real(1:3,i,j,k))
enddo; enddo; enddo end do; end do; end do
call utilities_FFTvectorForward call utilities_FFTvectorForward
call utilities_fourierVectorDivergence !< calculate temperature divergence in fourier field call utilities_fourierVectorDivergence !< calculate temperature divergence in fourier field
call utilities_FFTscalarBackward call utilities_FFTscalarBackward
@ -282,7 +282,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
scalarField_real(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_T(ce)) & scalarField_real(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_T(ce)) &
+ homogenization_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) & + homogenization_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) &
+ mu_ref*T_current(i,j,k) + mu_ref*T_current(i,j,k)
enddo; enddo; enddo end do; end do; end do
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! convolution of temperature field with green operator ! convolution of temperature field with green operator
@ -310,7 +310,7 @@ subroutine updateReference()
do ce = 1, product(grid(1:2))*grid3 do ce = 1, product(grid(1:2))*grid3
K_ref = K_ref + homogenization_K_T(ce) K_ref = K_ref + homogenization_K_T(ce)
mu_ref = mu_ref + homogenization_mu_T(ce) mu_ref = mu_ref + homogenization_mu_T(ce)
enddo end do
K_ref = K_ref*wgt K_ref = K_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr)

108
src/grid/spectral_utilities.f90
View File

@ -177,19 +177,19 @@ subroutine spectral_utilities_init
num_grid, & num_grid, &
debug_grid ! pointer to grid debug options debug_grid ! pointer to grid debug options
print'(/,a)', ' <<<+- spectral_utilities init -+>>>' print'(/,1x,a)', '<<<+- spectral_utilities init -+>>>'
print*, 'M. Diehl, Diploma Thesis TU München, 2010' print'(/,1x,a)', 'M. Diehl, Diploma Thesis TU München, 2010'
print*, 'https://doi.org/10.13140/2.1.3234.3840'//IO_EOL print'( 1x,a)', 'https://doi.org/10.13140/2.1.3234.3840'//IO_EOL
print*, 'P. Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013' print'( 1x,a)', 'P. Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013'
print*, 'https://doi.org/10.1016/j.ijplas.2012.09.012'//IO_EOL print'( 1x,a)', 'https://doi.org/10.1016/j.ijplas.2012.09.012'//IO_EOL
print*, 'P. Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' print'( 1x,a)', 'P. Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006'//IO_EOL print'( 1x,a)', 'https://doi.org/10.1016/j.ijplas.2014.02.006'//IO_EOL
print*, 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019' print'( 1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019'
print*, 'https://doi.org/10.1007/978-981-10-6855-3_80' print'( 1x,a)', 'https://doi.org/10.1007/978-981-10-6855-3_80'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set debugging parameters ! set debugging parameters
@ -200,15 +200,15 @@ subroutine spectral_utilities_init
debugRotation = debug_grid%contains('rotation') debugRotation = debug_grid%contains('rotation')
debugPETSc = debug_grid%contains('PETSc') debugPETSc = debug_grid%contains('PETSc')
if(debugPETSc) print'(3(/,a),/)', & if (debugPETSc) print'(3(/,1x,a),/)', &
' Initializing PETSc with debug options: ', & 'Initializing PETSc with debug options: ', &
trim(PETScDebug), & trim(PETScDebug), &
' add more using the "PETSc_options" keyword in numerics.yaml' 'add more using the "PETSc_options" keyword in numerics.yaml'
flush(IO_STDOUT) flush(IO_STDOUT)
call PetscOptionsClear(PETSC_NULL_OPTIONS,ierr) call PetscOptionsClear(PETSC_NULL_OPTIONS,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
if(debugPETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),ierr) if (debugPETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,& call PetscOptionsInsertString(PETSC_NULL_OPTIONS,&
num_grid%get_asString('PETSc_options',defaultVal=''),ierr) num_grid%get_asString('PETSc_options',defaultVal=''),ierr)
@ -271,7 +271,7 @@ subroutine spectral_utilities_init
if (pReal /= C_DOUBLE .or. kind(1) /= C_INT) error stop 'C and Fortran datatypes do not match' if (pReal /= C_DOUBLE .or. kind(1) /= C_INT) error stop 'C and Fortran datatypes do not match'
call fftw_set_timelimit(num_grid%get_asFloat('fftw_timelimit',defaultVal=-1.0_pReal)) call fftw_set_timelimit(num_grid%get_asFloat('fftw_timelimit',defaultVal=-1.0_pReal))
print*, 'FFTW initialized'; flush(IO_STDOUT) print'(/,1x,a)', 'FFTW initialized'; flush(IO_STDOUT)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! MPI allocation ! MPI allocation
@ -342,10 +342,10 @@ subroutine spectral_utilities_init
! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around)
do k = grid3Offset+1, grid3Offset+grid3 do k = grid3Offset+1, grid3Offset+grid3
k_s(3) = k - 1 k_s(3) = k - 1
if(k > grid(3)/2 + 1) k_s(3) = k_s(3) - grid(3) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1 if (k > grid(3)/2 + 1) k_s(3) = k_s(3) - grid(3) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1
do j = 1, grid(2) do j = 1, grid(2)
k_s(2) = j - 1 k_s(2) = j - 1
if(j > grid(2)/2 + 1) k_s(2) = k_s(2) - grid(2) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1 if (j > grid(2)/2 + 1) k_s(2) = k_s(2) - grid(2) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1
do i = 1, grid1Red do i = 1, grid1Red
k_s(1) = i - 1 ! symmetry, junst running from 0,1,...,N/2,N/2+1 k_s(1) = i - 1 ! symmetry, junst running from 0,1,...,N/2,N/2+1
xi2nd(1:3,i,j,k-grid3Offset) = utilities_getFreqDerivative(k_s) xi2nd(1:3,i,j,k-grid3Offset) = utilities_getFreqDerivative(k_s)
@ -357,7 +357,7 @@ subroutine spectral_utilities_init
endwhere endwhere
enddo; enddo; enddo enddo; enddo; enddo
if(num%memory_efficient) then ! allocate just single fourth order tensor if (num%memory_efficient) then ! allocate just single fourth order tensor
allocate (gamma_hat(3,3,3,3,1,1,1), source = cmplx(0.0_pReal,0.0_pReal,pReal)) allocate (gamma_hat(3,3,3,3,1,1,1), source = cmplx(0.0_pReal,0.0_pReal,pReal))
else ! precalculation of gamma_hat field else ! precalculation of gamma_hat field
allocate (gamma_hat(3,3,3,3,grid1Red,grid(2),grid3), source = cmplx(0.0_pReal,0.0_pReal,pReal)) allocate (gamma_hat(3,3,3,3,grid1Red,grid(2),grid3), source = cmplx(0.0_pReal,0.0_pReal,pReal))
@ -384,7 +384,7 @@ subroutine utilities_updateGamma(C)
C_ref = C C_ref = C
if(.not. num%memory_efficient) then if (.not. num%memory_efficient) then
gamma_hat = cmplx(0.0_pReal,0.0_pReal,pReal) ! for the singular point and any non invertible A gamma_hat = cmplx(0.0_pReal,0.0_pReal,pReal) ! for the singular point and any non invertible A
do k = grid3Offset+1, grid3Offset+grid3; do j = 1, grid(2); do i = 1, grid1Red do k = grid3Offset+1, grid3Offset+grid3; do j = 1, grid(2); do i = 1, grid1Red
if (any([i,j,k] /= 1)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 if (any([i,j,k] /= 1)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
@ -497,12 +497,12 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
logical :: err logical :: err
print'(/,a)', ' ... doing gamma convolution ...............................................' print'(/,1x,a)', '... doing gamma convolution ...............................................'
flush(IO_STDOUT) flush(IO_STDOUT)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! do the actual spectral method calculation (mechanical equilibrium) ! do the actual spectral method calculation (mechanical equilibrium)
memoryEfficient: if(num%memory_efficient) then memoryEfficient: if (num%memory_efficient) then
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid1Red do k = 1, grid3; do j = 1, grid(2); do i = 1, grid1Red
if (any([i,j,k+grid3Offset] /= 1)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 if (any([i,j,k+grid3Offset] /= 1)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
forall(l = 1:3, m = 1:3) & forall(l = 1:3, m = 1:3) &
@ -567,7 +567,7 @@ real(pReal) function utilities_divergenceRMS()
integer :: i, j, k, ierr integer :: i, j, k, ierr
complex(pReal), dimension(3) :: rescaledGeom complex(pReal), dimension(3) :: rescaledGeom
print'(/,a)', ' ... calculating divergence ................................................' print'(/,1x,a)', '... calculating divergence ................................................'
flush(IO_STDOUT) flush(IO_STDOUT)
rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal) rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal)
@ -593,9 +593,9 @@ real(pReal) function utilities_divergenceRMS()
+ sum(aimag(matmul(tensorField_fourier(1:3,1:3,grid1Red,j,k), & + sum(aimag(matmul(tensorField_fourier(1:3,1:3,grid1Red,j,k), &
conjg(-xi1st(1:3,grid1Red,j,k))*rescaledGeom))**2.0_pReal) conjg(-xi1st(1:3,grid1Red,j,k))*rescaledGeom))**2.0_pReal)
enddo; enddo enddo; enddo
if(grid(1) == 1) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of grid(1) == 1 if (grid(1) == 1) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of grid(1) == 1
call MPI_Allreduce(MPI_IN_PLACE,utilities_divergenceRMS,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,utilities_divergenceRMS,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr)
if(ierr /=0) error stop 'MPI error' if (ierr /=0) error stop 'MPI error'
utilities_divergenceRMS = sqrt(utilities_divergenceRMS) * wgt ! RMS in real space calculated with Parsevals theorem from Fourier space utilities_divergenceRMS = sqrt(utilities_divergenceRMS) * wgt ! RMS in real space calculated with Parsevals theorem from Fourier space
end function utilities_divergenceRMS end function utilities_divergenceRMS
@ -610,7 +610,7 @@ real(pReal) function utilities_curlRMS()
complex(pReal), dimension(3,3) :: curl_fourier complex(pReal), dimension(3,3) :: curl_fourier
complex(pReal), dimension(3) :: rescaledGeom complex(pReal), dimension(3) :: rescaledGeom
print'(/,a)', ' ... calculating curl ......................................................' print'(/,1x,a)', '... calculating curl ......................................................'
flush(IO_STDOUT) flush(IO_STDOUT)
rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal) rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal)
@ -655,9 +655,9 @@ real(pReal) function utilities_curlRMS()
enddo; enddo enddo; enddo
call MPI_Allreduce(MPI_IN_PLACE,utilities_curlRMS,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,utilities_curlRMS,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr)
if(ierr /=0) error stop 'MPI error' if (ierr /=0) error stop 'MPI error'
utilities_curlRMS = sqrt(utilities_curlRMS) * wgt utilities_curlRMS = sqrt(utilities_curlRMS) * wgt
if(grid(1) == 1) utilities_curlRMS = utilities_curlRMS * 0.5_pReal ! counted twice in case of grid(1) == 1 if (grid(1) == 1) utilities_curlRMS = utilities_curlRMS * 0.5_pReal ! counted twice in case of grid(1) == 1
end function utilities_curlRMS end function utilities_curlRMS
@ -686,13 +686,13 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
mask_stressVector = .not. reshape(transpose(mask_stress), [9]) mask_stressVector = .not. reshape(transpose(mask_stress), [9])
size_reduced = count(mask_stressVector) size_reduced = count(mask_stressVector)
if(size_reduced > 0) then if (size_reduced > 0) then
temp99_real = math_3333to99(rot_BC%rotate(C)) temp99_real = math_3333to99(rot_BC%rotate(C))
if(debugGeneral) then if (debugGeneral) then
print'(/,a)', ' ... updating masked compliance ............................................' print'(/,1x,a)', '... updating masked compliance ............................................'
print'(/,a,/,8(9(2x,f12.7,1x)/),9(2x,f12.7,1x))', & print'(/,1x,a,/,8(9(2x,f12.7,1x)/),9(2x,f12.7,1x))', &
' Stiffness C (load) / GPa =', transpose(temp99_Real)*1.0e-9_pReal 'Stiffness C (load) / GPa =', transpose(temp99_Real)*1.0e-9_pReal
flush(IO_STDOUT) flush(IO_STDOUT)
endif endif
@ -711,10 +711,10 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
errmatinv = errmatinv .or. any(dNeq(sTimesC,math_eye(size_reduced),1.0e-12_pReal)) errmatinv = errmatinv .or. any(dNeq(sTimesC,math_eye(size_reduced),1.0e-12_pReal))
if (debugGeneral .or. errmatinv) then if (debugGeneral .or. errmatinv) then
write(formatString, '(i2)') size_reduced write(formatString, '(i2)') size_reduced
formatString = '(/,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))' formatString = '(/,1x,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))'
print trim(formatString), ' C * S (load) ', transpose(matmul(c_reduced,s_reduced)) print trim(formatString), 'C * S (load) ', transpose(matmul(c_reduced,s_reduced))
print trim(formatString), ' S (load) ', transpose(s_reduced) print trim(formatString), 'S (load) ', transpose(s_reduced)
if(errmatinv) error stop 'matrix inversion error' if (errmatinv) error stop 'matrix inversion error'
endif endif
temp99_real = reshape(unpack(reshape(s_reduced,[size_reduced**2]),reshape(mask,[81]),0.0_pReal),[9,9]) temp99_real = reshape(unpack(reshape(s_reduced,[size_reduced**2]),reshape(mask,[81]),0.0_pReal),[9,9])
else else
@ -723,9 +723,9 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
utilities_maskedCompliance = math_99to3333(temp99_Real) utilities_maskedCompliance = math_99to3333(temp99_Real)
if(debugGeneral) then if (debugGeneral) then
print'(/,a,/,9(9(2x,f10.5,1x)/),9(2x,f10.5,1x))', & print'(/,1x,a,/,9(9(2x,f10.5,1x)/),9(2x,f10.5,1x))', &
' Masked Compliance (load) * GPa =', transpose(temp99_Real)*1.0e9_pReal 'Masked Compliance (load) * GPa =', transpose(temp99_Real)*1.0e9_pReal
flush(IO_STDOUT) flush(IO_STDOUT)
endif endif
@ -810,7 +810,7 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
real(pReal) :: dPdF_norm_max, dPdF_norm_min real(pReal) :: dPdF_norm_max, dPdF_norm_min
real(pReal), dimension(2) :: valueAndRank !< pair of min/max norm of dPdF to synchronize min/max of dPdF real(pReal), dimension(2) :: valueAndRank !< pair of min/max norm of dPdF to synchronize min/max of dPdF
print'(/,a)', ' ... evaluating constitutive response ......................................' print'(/,1x,a)', '... evaluating constitutive response ......................................'
flush(IO_STDOUT) flush(IO_STDOUT)
homogenization_F = reshape(F,[3,3,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field homogenization_F = reshape(F,[3,3,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field
@ -824,11 +824,11 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
P = reshape(homogenization_P, [3,3,grid(1),grid(2),grid3]) P = reshape(homogenization_P, [3,3,grid(1),grid(2),grid3])
P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt
call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr)
if (debugRotation) print'(/,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', & if (debugRotation) print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &
' Piola--Kirchhoff stress (lab) / MPa =', transpose(P_av)*1.e-6_pReal 'Piola--Kirchhoff stress (lab) / MPa =', transpose(P_av)*1.e-6_pReal
if(present(rotation_BC)) P_av = rotation_BC%rotate(P_av) if (present(rotation_BC)) P_av = rotation_BC%rotate(P_av)
print'(/,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', & print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &
' Piola--Kirchhoff stress / MPa =', transpose(P_av)*1.e-6_pReal 'Piola--Kirchhoff stress / MPa =', transpose(P_av)*1.e-6_pReal
flush(IO_STDOUT) flush(IO_STDOUT)
dPdF_max = 0.0_pReal dPdF_max = 0.0_pReal
@ -1017,7 +1017,7 @@ subroutine utilities_updateCoords(F)
call utilities_FFTtensorForward() call utilities_FFTtensorForward()
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid1Red do k = 1, grid3; do j = 1, grid(2); do i = 1, grid1Red
if(any([i,j,k+grid3Offset] /= 1)) then if (any([i,j,k+grid3Offset] /= 1)) then
vectorField_fourier(1:3,i,j,k) = matmul(tensorField_fourier(1:3,1:3,i,j,k),xi2nd(1:3,i,j,k)) & vectorField_fourier(1:3,i,j,k) = matmul(tensorField_fourier(1:3,1:3,i,j,k),xi2nd(1:3,i,j,k)) &
/ sum(conjg(-xi2nd(1:3,i,j,k))*xi2nd(1:3,i,j,k)) * cmplx(wgt,0.0,pReal) / sum(conjg(-xi2nd(1:3,i,j,k))*xi2nd(1:3,i,j,k)) * cmplx(wgt,0.0,pReal)
else else
@ -1031,7 +1031,7 @@ subroutine utilities_updateCoords(F)
! average F ! average F
if (grid3Offset == 0) Favg = real(tensorField_fourier(1:3,1:3,1,1,1),pReal)*wgt if (grid3Offset == 0) Favg = real(tensorField_fourier(1:3,1:3,1,1,1),pReal)*wgt
call MPI_Bcast(Favg,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr) call MPI_Bcast(Favg,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
if(ierr /=0) error stop 'MPI error' if (ierr /=0) error stop 'MPI error'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! pad cell center fluctuations along z-direction (needed when running MPI simulation) ! pad cell center fluctuations along z-direction (needed when running MPI simulation)
@ -1042,22 +1042,22 @@ subroutine utilities_updateCoords(F)
! send bottom layer to process below ! send bottom layer to process below
call MPI_Isend(IPfluct_padded(:,:,:,2), c,MPI_DOUBLE,rank_b,0,MPI_COMM_WORLD,request(1),ierr) call MPI_Isend(IPfluct_padded(:,:,:,2), c,MPI_DOUBLE,rank_b,0,MPI_COMM_WORLD,request(1),ierr)
if(ierr /=0) error stop 'MPI error' if (ierr /=0) error stop 'MPI error'
call MPI_Irecv(IPfluct_padded(:,:,:,grid3+2),c,MPI_DOUBLE,rank_t,0,MPI_COMM_WORLD,request(2),ierr) call MPI_Irecv(IPfluct_padded(:,:,:,grid3+2),c,MPI_DOUBLE,rank_t,0,MPI_COMM_WORLD,request(2),ierr)
if(ierr /=0) error stop 'MPI error' if (ierr /=0) error stop 'MPI error'
! send top layer to process above ! send top layer to process above
call MPI_Isend(IPfluct_padded(:,:,:,grid3+1),c,MPI_DOUBLE,rank_t,1,MPI_COMM_WORLD,request(3),ierr) call MPI_Isend(IPfluct_padded(:,:,:,grid3+1),c,MPI_DOUBLE,rank_t,1,MPI_COMM_WORLD,request(3),ierr)
if(ierr /=0) error stop 'MPI error' if (ierr /=0) error stop 'MPI error'
call MPI_Irecv(IPfluct_padded(:,:,:,1), c,MPI_DOUBLE,rank_b,1,MPI_COMM_WORLD,request(4),ierr) call MPI_Irecv(IPfluct_padded(:,:,:,1), c,MPI_DOUBLE,rank_b,1,MPI_COMM_WORLD,request(4),ierr)
if(ierr /=0) error stop 'MPI error' if (ierr /=0) error stop 'MPI error'
call MPI_Waitall(4,request,status,ierr) call MPI_Waitall(4,request,status,ierr)
if(ierr /=0) error stop 'MPI error' if (ierr /=0) error stop 'MPI error'
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
! ToDo ! ToDo
#else #else
if(any(status(MPI_ERROR,:) /= 0)) error stop 'MPI error' if (any(status(MPI_ERROR,:) /= 0)) error stop 'MPI error'
#endif #endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -1094,10 +1094,10 @@ subroutine utilities_saveReferenceStiffness
fileUnit,ierr fileUnit,ierr
if (worldrank == 0) then if (worldrank == 0) then
print'(a)', ' writing reference stiffness data required for restart to file'; flush(IO_STDOUT) print'(1x,a)', 'writing reference stiffness data required for restart to file'; flush(IO_STDOUT)
open(newunit=fileUnit, file=getSolverJobName()//'.C_ref',& open(newunit=fileUnit, file=getSolverJobName()//'.C_ref',&
status='replace',access='stream',action='write',iostat=ierr) status='replace',access='stream',action='write',iostat=ierr)
if(ierr /=0) call IO_error(100,ext_msg='could not open file '//getSolverJobName()//'.C_ref') if (ierr /=0) call IO_error(100,ext_msg='could not open file '//getSolverJobName()//'.C_ref')
write(fileUnit) C_ref write(fileUnit) C_ref
close(fileUnit) close(fileUnit)
endif endif

2
src/homogenization.f90
View File

@ -199,7 +199,7 @@ subroutine homogenization_init()
num_homog, & num_homog, &
num_homogGeneric num_homogGeneric
print'(/,a)', ' <<<+- homogenization init -+>>>'; flush(IO_STDOUT) print'(/,1x,a)', '<<<+- homogenization init -+>>>'; flush(IO_STDOUT)
allocate(homogState (size(material_name_homogenization))) allocate(homogState (size(material_name_homogenization)))

2
src/homogenization_damage.f90
View File

@ -41,7 +41,7 @@ module subroutine damage_init()
integer :: ho,Nmembers integer :: ho,Nmembers
print'(/,a)', ' <<<+- homogenization:damage init -+>>>' print'(/,1x,a)', '<<<+- homogenization:damage init -+>>>'
configHomogenizations => config_material%get('homogenization') configHomogenizations => config_material%get('homogenization')

2
src/homogenization_damage_pass.f90
View File

@ -8,7 +8,7 @@ contains
module subroutine pass_init() module subroutine pass_init()
print'(/,a)', ' <<<+- homogenization:damage:pass init -+>>>' print'(/,1x,a)', '<<<+- homogenization:damage:pass init -+>>>'
end subroutine pass_init end subroutine pass_init

12
src/homogenization_mechanical.f90
View File

@ -68,7 +68,7 @@ module subroutine mechanical_init(num_homog)
class(tNode), pointer :: & class(tNode), pointer :: &
num_homogMech num_homogMech
print'(/,a)', ' <<<+- homogenization:mechanical init -+>>>' print'(/,1x,a)', '<<<+- homogenization:mechanical init -+>>>'
call material_parseHomogenization2() call material_parseHomogenization2()
@ -114,7 +114,7 @@ module subroutine mechanical_partition(subF,ce)
do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
call phase_set_F(Fs(1:3,1:3,co),co,ce) call phase_set_F(Fs(1:3,1:3,co),co,ce)
enddo end do
end subroutine mechanical_partition end subroutine mechanical_partition
@ -138,7 +138,7 @@ module subroutine mechanical_homogenize(Delta_t,ce)
+ phase_P(co,ce) + phase_P(co,ce)
homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = homogenization_dPdF(1:3,1:3,1:3,1:3,ce) & homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = homogenization_dPdF(1:3,1:3,1:3,1:3,ce) &
+ phase_mechanical_dPdF(Delta_t,co,ce) + phase_mechanical_dPdF(Delta_t,co,ce)
enddo end do
homogenization_P(1:3,1:3,ce) = homogenization_P(1:3,1:3,ce) & homogenization_P(1:3,1:3,ce) = homogenization_P(1:3,1:3,ce) &
/ real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal) / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
@ -173,11 +173,11 @@ module function mechanical_updateState(subdt,subF,ce) result(doneAndHappy)
dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(subdt,co,ce) dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(subdt,co,ce)
Fs(:,:,co) = phase_F(co,ce) Fs(:,:,co) = phase_F(co,ce)
Ps(:,:,co) = phase_P(co,ce) Ps(:,:,co) = phase_P(co,ce)
enddo end do
doneAndHappy = RGC_updateState(Ps,Fs,subF,subdt,dPdFs,ce) doneAndHappy = RGC_updateState(Ps,Fs,subF,subdt,dPdFs,ce)
else else
doneAndHappy = .true. doneAndHappy = .true.
endif end if
end function mechanical_updateState end function mechanical_updateState
@ -241,7 +241,7 @@ subroutine material_parseHomogenization2()
case default case default
call IO_error(500,ext_msg=homogMech%get_asString('type')) call IO_error(500,ext_msg=homogMech%get_asString('type'))
end select end select
enddo end do
end subroutine material_parseHomogenization2 end subroutine material_parseHomogenization2

90
src/homogenization_mechanical_RGC.f90
View File

@ -87,16 +87,16 @@ module subroutine RGC_init(num_homogMech)
homog, & homog, &
homogMech homogMech
print'(/,a)', ' <<<+- homogenization:mechanical:RGC init -+>>>' print'(/,1x,a)', '<<<+- homogenization:mechanical:RGC init -+>>>'
print'(a,i0)', ' # homogenizations: ',count(homogenization_type == HOMOGENIZATION_RGC_ID) print'(/,a,i0)', ' # homogenizations: ',count(homogenization_type == HOMOGENIZATION_RGC_ID)
flush(IO_STDOUT) flush(IO_STDOUT)
print*, 'D.D. Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009' print'(/,1x,a)', 'D.D. Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009'
print*, 'https://doi.org/10.1007/s12289-009-0619-1'//IO_EOL print'( 1x,a)', 'https://doi.org/10.1007/s12289-009-0619-1'//IO_EOL
print*, 'D.D. Tjahjanto et al., Modelling and Simulation in Materials Science and Engineering 18:015006, 2010' print'(/,1x,a)', 'D.D. Tjahjanto et al., Modelling and Simulation in Materials Science and Engineering 18:015006, 2010'
print*, 'https://doi.org/10.1088/0965-0393/18/1/015006'//IO_EOL print'( 1x,a)', 'https://doi.org/10.1088/0965-0393/18/1/015006'//IO_EOL
material_homogenization => config_material%get('homogenization') material_homogenization => config_material%get('homogenization')
@ -186,7 +186,7 @@ module subroutine RGC_init(num_homogMech)
end associate end associate
enddo end do
end subroutine RGC_init end subroutine RGC_init
@ -222,9 +222,9 @@ module subroutine RGC_partitionDeformation(F,avgF,ce)
nVect = interfaceNormal(intFace,ho,en) nVect = interfaceNormal(intFace,ho,en)
forall (i=1:3,j=1:3) & forall (i=1:3,j=1:3) &
F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! calculating deformation relaxations due to interface relaxation F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! calculating deformation relaxations due to interface relaxation
enddo end do
F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF ! resulting relaxed deformation gradient F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF ! resulting relaxed deformation gradient
enddo end do
end associate end associate
@ -257,10 +257,10 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
real(pReal), dimension(:,:), allocatable :: tract,jmatrix,jnverse,smatrix,pmatrix,rmatrix real(pReal), dimension(:,:), allocatable :: tract,jmatrix,jnverse,smatrix,pmatrix,rmatrix
real(pReal), dimension(:), allocatable :: resid,relax,p_relax,p_resid,drelax real(pReal), dimension(:), allocatable :: resid,relax,p_relax,p_resid,drelax
zeroTimeStep: if(dEq0(dt)) then zeroTimeStep: if (dEq0(dt)) then
doneAndHappy = .true. ! pretend everything is fine and return doneAndHappy = .true. ! pretend everything is fine and return
return return
endif zeroTimeStep end if zeroTimeStep
ho = material_homogenizationID(ce) ho = material_homogenizationID(ce)
en = material_homogenizationEntry(ce) en = material_homogenizationEntry(ce)
@ -319,10 +319,10 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
tract(iNum,i) = tract(iNum,i) + (P(i,j,iGrP) + R(i,j,iGrP) + D(i,j,iGrP))*normP(j) & ! contribution from material stress P, mismatch penalty R, and volume penalty D projected into the interface tract(iNum,i) = tract(iNum,i) + (P(i,j,iGrP) + R(i,j,iGrP) + D(i,j,iGrP))*normP(j) & ! contribution from material stress P, mismatch penalty R, and volume penalty D projected into the interface
+ (P(i,j,iGrN) + R(i,j,iGrN) + D(i,j,iGrN))*normN(j) + (P(i,j,iGrN) + R(i,j,iGrN) + D(i,j,iGrN))*normN(j)
resid(i+3*(iNum-1)) = tract(iNum,i) ! translate the local residual into global 1-dimensional residual array resid(i+3*(iNum-1)) = tract(iNum,i) ! translate the local residual into global 1-dimensional residual array
enddo end do
enddo end do
enddo end do
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! convergence check for stress residual ! convergence check for stress residual
@ -347,7 +347,7 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
elseif (residMax > num%relMax*stresMax .or. residMax > num%absMax) then ! try to restart when residual blows up exceeding maximum bound elseif (residMax > num%relMax*stresMax .or. residMax > num%absMax) then ! try to restart when residual blows up exceeding maximum bound
doneAndHappy = [.true.,.false.] ! with direct cut-back doneAndHappy = [.true.,.false.] ! with direct cut-back
return return
endif end if
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
! construct the global Jacobian matrix for updating the global relaxation vector array when ! construct the global Jacobian matrix for updating the global relaxation vector array when
@ -373,11 +373,11 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
do i=1,3; do j=1,3; do k=1,3; do l=1,3 do i=1,3; do j=1,3; do k=1,3; do l=1,3
smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) & smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) &
+ dPdF(i,k,j,l,iGrN)*normN(k)*mornN(l) + dPdF(i,k,j,l,iGrN)*normN(k)*mornN(l)
enddo;enddo;enddo;enddo end do;enddo;enddo;enddo
! projecting the material tangent dPdF into the interface ! projecting the material tangent dPdF into the interface
! to obtain the Jacobian matrix contribution of dPdF ! to obtain the Jacobian matrix contribution of dPdF
endif end if
enddo end do
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! identify the right/up/front grain (+|P) ! identify the right/up/front grain (+|P)
@ -394,10 +394,10 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
do i=1,3; do j=1,3; do k=1,3; do l=1,3 do i=1,3; do j=1,3; do k=1,3; do l=1,3
smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) & smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) &
+ dPdF(i,k,j,l,iGrP)*normP(k)*mornP(l) + dPdF(i,k,j,l,iGrP)*normP(k)*mornP(l)
enddo;enddo;enddo;enddo end do;enddo;enddo;enddo
endif end if
enddo end do
enddo end do
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! ... of the stress penalty tangent (mismatch penalty and volume penalty, computed using numerical ! ... of the stress penalty tangent (mismatch penalty and volume penalty, computed using numerical
@ -443,10 +443,10 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
+ (pR(i,j,iGrN) - R(i,j,iGrN))*normN(j) & + (pR(i,j,iGrN) - R(i,j,iGrN))*normN(j) &
+ (pD(i,j,iGrP) - D(i,j,iGrP))*normP(j) & + (pD(i,j,iGrP) - D(i,j,iGrP))*normP(j) &
+ (pD(i,j,iGrN) - D(i,j,iGrN))*normN(j) + (pD(i,j,iGrN) - D(i,j,iGrN))*normN(j)
enddo; enddo end do; end do
enddo end do
pmatrix(:,ipert) = p_resid/num%pPert pmatrix(:,ipert) = p_resid/num%pPert
enddo end do
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -455,7 +455,7 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
do i=1,3*nIntFaceTot do i=1,3*nIntFaceTot
rmatrix(i,i) = num%viscModus*num%viscPower/(num%refRelaxRate*dt)* & ! tangent due to numerical viscosity traction appears rmatrix(i,i) = num%viscModus*num%viscPower/(num%refRelaxRate*dt)* & ! tangent due to numerical viscosity traction appears
(abs(drelax(i))/(num%refRelaxRate*dt))**(num%viscPower - 1.0_pReal) ! only in the main diagonal term (abs(drelax(i))/(num%refRelaxRate*dt))**(num%viscPower - 1.0_pReal) ! only in the main diagonal term
enddo end do
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -472,7 +472,7 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
drelax = 0.0_pReal drelax = 0.0_pReal
do i = 1,3*nIntFaceTot;do j = 1,3*nIntFaceTot do i = 1,3*nIntFaceTot;do j = 1,3*nIntFaceTot
drelax(i) = drelax(i) - jnverse(i,j)*resid(j) ! Calculate the correction for the state variable drelax(i) = drelax(i) - jnverse(i,j)*resid(j) ! Calculate the correction for the state variable
enddo; enddo end do; end do
stt%relaxationVector(:,en) = relax + drelax ! Updateing the state variable for the next iteration stt%relaxationVector(:,en) = relax + drelax ! Updateing the state variable for the next iteration
if (any(abs(drelax) > num%maxdRelax)) then ! Forcing cutback when the incremental change of relaxation vector becomes too large if (any(abs(drelax) > num%maxdRelax)) then ! Forcing cutback when the incremental change of relaxation vector becomes too large
doneAndHappy = [.true.,.false.] doneAndHappy = [.true.,.false.]
@ -481,7 +481,7 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
print'(a,e15.8)',' due to large relaxation change = ',maxval(abs(drelax)) print'(a,e15.8)',' due to large relaxation change = ',maxval(abs(drelax))
flush(IO_STDOUT) flush(IO_STDOUT)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif end if
end associate end associate
@ -545,9 +545,9 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
do i = 1,3; do j = 1,3 do i = 1,3; do j = 1,3
do k = 1,3; do l = 1,3 do k = 1,3; do l = 1,3
nDef(i,j) = nDef(i,j) - nVect(k)*gDef(i,l)*math_LeviCivita(j,k,l) ! compute the interface mismatch tensor from the jump of deformation gradient nDef(i,j) = nDef(i,j) - nVect(k)*gDef(i,l)*math_LeviCivita(j,k,l) ! compute the interface mismatch tensor from the jump of deformation gradient
enddo; enddo end do; end do
nDefNorm = nDefNorm + nDef(i,j)**2.0_pReal ! compute the norm of the mismatch tensor nDefNorm = nDefNorm + nDef(i,j)**2.0_pReal ! compute the norm of the mismatch tensor
enddo; enddo end do; end do
nDefNorm = max(nDefToler,sqrt(nDefNorm)) ! approximation to zero mismatch if mismatch is zero (singularity) nDefNorm = max(nDefToler,sqrt(nDefNorm)) ! approximation to zero mismatch if mismatch is zero (singularity)
nMis(abs(intFace(1)),iGrain) = nMis(abs(intFace(1)),iGrain) + nDefNorm ! total amount of mismatch experienced by the grain (at all six interfaces) nMis(abs(intFace(1)),iGrain) = nMis(abs(intFace(1)),iGrain) + nDefNorm ! total amount of mismatch experienced by the grain (at all six interfaces)
@ -560,11 +560,11 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
*cosh(prm%c_alpha*nDefNorm) & *cosh(prm%c_alpha*nDefNorm) &
*0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_LeviCivita(k,l,j) & *0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_LeviCivita(k,l,j) &
*tanh(nDefNorm/num%xSmoo) *tanh(nDefNorm/num%xSmoo)
enddo; enddo;enddo; enddo end do; end do;enddo; end do
enddo interfaceLoop end do interfaceLoop
enddo grainLoop end do grainLoop
end associate end associate
@ -594,7 +594,7 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
gVol(i) = math_det33(fDef(1:3,1:3,i)) ! compute the volume of individual grains gVol(i) = math_det33(fDef(1:3,1:3,i)) ! compute the volume of individual grains
vDiscrep = vDiscrep - gVol(i)/real(nGrain,pReal) ! calculate the difference/dicrepancy between vDiscrep = vDiscrep - gVol(i)/real(nGrain,pReal) ! calculate the difference/dicrepancy between
! the volume of the cluster and the the total volume of grains ! the volume of the cluster and the the total volume of grains
enddo end do
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
! calculate the stress and penalty due to volume discrepancy ! calculate the stress and penalty due to volume discrepancy
@ -603,7 +603,7 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
vPen(:,:,i) = -1.0_pReal/real(nGrain,pReal)*num%volDiscrMod*num%volDiscrPow/num%maxVolDiscr* & vPen(:,:,i) = -1.0_pReal/real(nGrain,pReal)*num%volDiscrMod*num%volDiscrPow/num%maxVolDiscr* &
sign((abs(vDiscrep)/num%maxVolDiscr)**(num%volDiscrPow - 1.0),vDiscrep)* & sign((abs(vDiscrep)/num%maxVolDiscr)**(num%volDiscrPow - 1.0),vDiscrep)* &
gVol(i)*transpose(math_inv33(fDef(:,:,i))) gVol(i)*transpose(math_inv33(fDef(:,:,i)))
enddo end do
end subroutine volumePenalty end subroutine volumePenalty
@ -633,9 +633,9 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
nVect = interfaceNormal([iBase,1,1,1],ho,en) nVect = interfaceNormal([iBase,1,1,1],ho,en)
do i = 1,3; do j = 1,3 do i = 1,3; do j = 1,3
surfaceCorrection(iBase) = surfaceCorrection(iBase) + invC(i,j)*nVect(i)*nVect(j) ! compute the component of (the inverse of) the stretch in the direction of the normal surfaceCorrection(iBase) = surfaceCorrection(iBase) + invC(i,j)*nVect(i)*nVect(j) ! compute the component of (the inverse of) the stretch in the direction of the normal
enddo; enddo end do; end do
surfaceCorrection(iBase) = sqrt(surfaceCorrection(iBase))*detF ! get the surface correction factor (area contraction/enlargement) surfaceCorrection(iBase) = sqrt(surfaceCorrection(iBase))*detF ! get the surface correction factor (area contraction/enlargement)
enddo end do
end function surfaceCorrection end function surfaceCorrection
@ -690,9 +690,9 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
nVect = interfaceNormal(intFace,ho,en) nVect = interfaceNormal(intFace,ho,en)
forall (i=1:3,j=1:3) & forall (i=1:3,j=1:3) &
F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations
enddo end do
F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF ! relaxed deformation gradient F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF ! relaxed deformation gradient
enddo end do
end associate end associate
@ -727,7 +727,7 @@ module subroutine RGC_results(ho,group)
call results_writeDataset(dst%relaxationrate_avg,group,trim(prm%output(o)), & call results_writeDataset(dst%relaxationrate_avg,group,trim(prm%output(o)), &
'average relaxation rate','m/s') 'average relaxation rate','m/s')
end select end select
enddo outputsLoop end do outputsLoop
end associate end associate
end subroutine RGC_results end subroutine RGC_results
@ -756,7 +756,7 @@ pure function relaxationVector(intFace,ho,en)
relaxationVector = stt%relaxationVector((3*iNum-2):(3*iNum),en) relaxationVector = stt%relaxationVector((3*iNum-2):(3*iNum),en)
else else
relaxationVector = 0.0_pReal relaxationVector = 0.0_pReal
endif end if
end associate end associate
@ -855,7 +855,7 @@ integer pure function interface4to1(iFace4D, nGDim)
else else
interface4to1 = iFace4D(3) + nGDim(2)*(iFace4D(4)-1) & interface4to1 = iFace4D(3) + nGDim(2)*(iFace4D(4)-1) &
+ nGDim(2)*nGDim(3)*(iFace4D(2)-1) + nGDim(2)*nGDim(3)*(iFace4D(2)-1)
endif end if
case(2) case(2)
if ((iFace4D(3) == 0) .or. (iFace4D(3) == nGDim(2))) then if ((iFace4D(3) == 0) .or. (iFace4D(3) == nGDim(2))) then
@ -864,7 +864,7 @@ integer pure function interface4to1(iFace4D, nGDim)
interface4to1 = iFace4D(4) + nGDim(3)*(iFace4D(2)-1) & interface4to1 = iFace4D(4) + nGDim(3)*(iFace4D(2)-1) &
+ nGDim(3)*nGDim(1)*(iFace4D(3)-1) & + nGDim(3)*nGDim(1)*(iFace4D(3)-1) &
+ (nGDim(1)-1)*nGDim(2)*nGDim(3) ! total # of interfaces normal || e1 + (nGDim(1)-1)*nGDim(2)*nGDim(3) ! total # of interfaces normal || e1
endif end if
case(3) case(3)
if ((iFace4D(4) == 0) .or. (iFace4D(4) == nGDim(3))) then if ((iFace4D(4) == 0) .or. (iFace4D(4) == nGDim(3))) then
@ -874,7 +874,7 @@ integer pure function interface4to1(iFace4D, nGDim)
+ nGDim(1)*nGDim(2)*(iFace4D(4)-1) & + nGDim(1)*nGDim(2)*(iFace4D(4)-1) &
+ (nGDim(1)-1)*nGDim(2)*nGDim(3) & ! total # of interfaces normal || e1 + (nGDim(1)-1)*nGDim(2)*nGDim(3) & ! total # of interfaces normal || e1
+ nGDim(1)*(nGDim(2)-1)*nGDim(3) ! total # of interfaces normal || e2 + nGDim(1)*(nGDim(2)-1)*nGDim(3) ! total # of interfaces normal || e2
endif end if
case default case default
interface4to1 = -1 interface4to1 = -1
@ -918,7 +918,7 @@ pure function interface1to4(iFace1D, nGDim)
interface1to4(2) = mod((iFace1D-nIntFace(2)-nIntFace(1)-1),nGDim(1))+1 interface1to4(2) = mod((iFace1D-nIntFace(2)-nIntFace(1)-1),nGDim(1))+1
interface1to4(3) = mod(int(real(iFace1D-nIntFace(2)-nIntFace(1)-1,pReal)/real(nGDim(1),pReal)),nGDim(2))+1 interface1to4(3) = mod(int(real(iFace1D-nIntFace(2)-nIntFace(1)-1,pReal)/real(nGDim(1),pReal)),nGDim(2))+1
interface1to4(4) = int(real(iFace1D-nIntFace(2)-nIntFace(1)-1,pReal)/real(nGDim(1),pReal)/real(nGDim(2),pReal))+1 interface1to4(4) = int(real(iFace1D-nIntFace(2)-nIntFace(1)-1,pReal)/real(nGDim(1),pReal)/real(nGDim(2),pReal))+1
endif end if
end function interface1to4 end function interface1to4

6
src/homogenization_mechanical_isostrain.f90
View File

@ -17,9 +17,9 @@ module subroutine isostrain_init
ho, & ho, &
Nmembers Nmembers
print'(/,a)', ' <<<+- homogenization:mechanical:isostrain init -+>>>' print'(/,1x,a)', '<<<+- homogenization:mechanical:isostrain init -+>>>'
print'(a,i0)', ' # homogenizations: ',count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID) print'(/,a,i0)', ' # homogenizations: ',count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)
flush(IO_STDOUT) flush(IO_STDOUT)
do ho = 1, size(homogenization_type) do ho = 1, size(homogenization_type)
@ -30,7 +30,7 @@ module subroutine isostrain_init
allocate(homogState(ho)%state0(0,Nmembers)) allocate(homogState(ho)%state0(0,Nmembers))
allocate(homogState(ho)%state (0,Nmembers)) allocate(homogState(ho)%state (0,Nmembers))
enddo end do
end subroutine isostrain_init end subroutine isostrain_init

10
src/homogenization_mechanical_pass.f90
View File

@ -17,15 +17,15 @@ module subroutine pass_init
ho, & ho, &
Nmembers Nmembers
print'(/,a)', ' <<<+- homogenization:mechanical:pass init -+>>>' print'(/,1x,a)', '<<<+- homogenization:mechanical:pass init -+>>>'
print'(a,i0)', ' # homogenizations: ',count(homogenization_type == HOMOGENIZATION_NONE_ID) print'(/,a,i0)', ' # homogenizations: ',count(homogenization_type == HOMOGENIZATION_NONE_ID)
flush(IO_STDOUT) flush(IO_STDOUT)
do ho = 1, size(homogenization_type) do ho = 1, size(homogenization_type)
if(homogenization_type(ho) /= HOMOGENIZATION_NONE_ID) cycle if (homogenization_type(ho) /= HOMOGENIZATION_NONE_ID) cycle
if(homogenization_Nconstituents(ho) /= 1) & if (homogenization_Nconstituents(ho) /= 1) &
call IO_error(211,ext_msg='N_constituents (pass)') call IO_error(211,ext_msg='N_constituents (pass)')
Nmembers = count(material_homogenizationID == ho) Nmembers = count(material_homogenizationID == ho)
@ -33,7 +33,7 @@ module subroutine pass_init
allocate(homogState(ho)%state0(0,Nmembers)) allocate(homogState(ho)%state0(0,Nmembers))
allocate(homogState(ho)%state (0,Nmembers)) allocate(homogState(ho)%state (0,Nmembers))
enddo end do
end subroutine pass_init end subroutine pass_init

16
src/homogenization_thermal.f90
View File

@ -44,7 +44,7 @@ module subroutine thermal_init()
integer :: ho integer :: ho
print'(/,a)', ' <<<+- homogenization:thermal init -+>>>' print'(/,1x,a)', '<<<+- homogenization:thermal init -+>>>'
configHomogenizations => config_material%get('homogenization') configHomogenizations => config_material%get('homogenization')
@ -65,9 +65,9 @@ module subroutine thermal_init()
#endif #endif
else else
prm%output = emptyStringArray prm%output = emptyStringArray
endif end if
end associate end associate
enddo end do
call pass_init() call pass_init()
@ -89,7 +89,7 @@ module subroutine thermal_partition(ce)
dot_T = current(material_homogenizationID(ce))%dot_T(material_homogenizationEntry(ce)) dot_T = current(material_homogenizationID(ce))%dot_T(material_homogenizationEntry(ce))
do co = 1, homogenization_Nconstituents(material_homogenizationID(ce)) do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
call phase_thermal_setField(T,dot_T,co,ce) call phase_thermal_setField(T,dot_T,co,ce)
enddo end do
end subroutine thermal_partition end subroutine thermal_partition
@ -108,7 +108,7 @@ module function homogenization_mu_T(ce) result(mu)
mu = phase_mu_T(1,ce) mu = phase_mu_T(1,ce)
do co = 2, homogenization_Nconstituents(material_homogenizationID(ce)) do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
mu = mu + phase_mu_T(co,ce) mu = mu + phase_mu_T(co,ce)
enddo end do
mu = mu / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal) mu = mu / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
@ -129,7 +129,7 @@ module function homogenization_K_T(ce) result(K)
K = phase_K_T(1,ce) K = phase_K_T(1,ce)
do co = 2, homogenization_Nconstituents(material_homogenizationID(ce)) do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
K = K + phase_K_T(co,ce) K = K + phase_K_T(co,ce)
enddo end do
K = K / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal) K = K / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
@ -150,7 +150,7 @@ module function homogenization_f_T(ce) result(f)
f = phase_f_T(material_phaseID(1,ce),material_phaseEntry(1,ce)) f = phase_f_T(material_phaseID(1,ce),material_phaseEntry(1,ce))
do co = 2, homogenization_Nconstituents(material_homogenizationID(ce)) do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
f = f + phase_f_T(material_phaseID(co,ce),material_phaseEntry(co,ce)) f = f + phase_f_T(material_phaseID(co,ce),material_phaseEntry(co,ce))
enddo end do
f = f/real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal) f = f/real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
@ -189,7 +189,7 @@ module subroutine thermal_results(ho,group)
case('T') case('T')
call results_writeDataset(current(ho)%T,group,'T','temperature','K') call results_writeDataset(current(ho)%T,group,'T','temperature','K')
end select end select
enddo outputsLoop end do outputsLoop
end associate end associate
end subroutine thermal_results end subroutine thermal_results

2
src/homogenization_thermal_pass.f90
View File

@ -8,7 +8,7 @@ contains
module subroutine pass_init() module subroutine pass_init()
print'(/,a)', ' <<<+- homogenization:thermal:pass init -+>>>' print'(/,1x,a)', '<<<+- homogenization:thermal:pass init -+>>>'
end subroutine pass_init end subroutine pass_init

2
src/lattice.f90
View File

@ -409,7 +409,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine lattice_init subroutine lattice_init
print'(/,a)', ' <<<+- lattice init -+>>>'; flush(IO_STDOUT) print'(/,1x,a)', '<<<+- lattice init -+>>>'; flush(IO_STDOUT)
call selfTest call selfTest

30
src/material.f90
View File

@ -61,11 +61,11 @@ subroutine material_init(restart)
logical, intent(in) :: restart logical, intent(in) :: restart
print'(/,a)', ' <<<+- material init -+>>>'; flush(IO_STDOUT) print'(/,1x,a)', '<<<+- material init -+>>>'; flush(IO_STDOUT)
call parse call parse
print*, 'parsed material.yaml' print'(/,1x,a)', 'parsed material.yaml'
if (.not. restart) then if (.not. restart) then
@ -73,7 +73,7 @@ subroutine material_init(restart)
call results_mapping_phase(material_phaseID,material_phaseEntry,material_name_phase) call results_mapping_phase(material_phaseID,material_phaseEntry,material_name_phase)
call results_mapping_homogenization(material_homogenizationID,material_homogenizationEntry,material_name_homogenization) call results_mapping_homogenization(material_homogenizationID,material_homogenizationEntry,material_name_homogenization)
call results_closeJobFile call results_closeJobFile
endif end if
end subroutine material_init end subroutine material_init
@ -113,7 +113,7 @@ subroutine parse()
do h=1, homogenizations%length do h=1, homogenizations%length
homogenization => homogenizations%get(h) homogenization => homogenizations%get(h)
homogenization_Nconstituents(h) = homogenization%get_asInt('N_constituents') homogenization_Nconstituents(h) = homogenization%get_asInt('N_constituents')
enddo end do
homogenization_maxNconstituents = maxval(homogenization_Nconstituents) homogenization_maxNconstituents = maxval(homogenization_Nconstituents)
allocate(counterPhase(phases%length),source=0) allocate(counterPhase(phases%length),source=0)
@ -137,7 +137,7 @@ subroutine parse()
material_homogenizationID(ce) = homogenizations%getIndex(material%get_asString('homogenization')) material_homogenizationID(ce) = homogenizations%getIndex(material%get_asString('homogenization'))
counterHomogenization(material_homogenizationID(ce)) = counterHomogenization(material_homogenizationID(ce)) + 1 counterHomogenization(material_homogenizationID(ce)) = counterHomogenization(material_homogenizationID(ce)) + 1
material_homogenizationEntry(ce) = counterHomogenization(material_homogenizationID(ce)) material_homogenizationEntry(ce) = counterHomogenization(material_homogenizationID(ce))
enddo end do
frac = 0.0_pReal frac = 0.0_pReal
do co = 1, constituents%length do co = 1, constituents%length
@ -151,12 +151,12 @@ subroutine parse()
material_phaseMemberAt(co,ip,el) = counterPhase(material_phaseAt(co,el)) material_phaseMemberAt(co,ip,el) = counterPhase(material_phaseAt(co,el))
material_phaseEntry(co,ce) = counterPhase(material_phaseAt(co,el)) material_phaseEntry(co,ce) = counterPhase(material_phaseAt(co,el))
material_phaseID(co,ce) = material_phaseAt(co,el) material_phaseID(co,ce) = material_phaseAt(co,el)
enddo end do
enddo end do
if (dNeq(frac,1.0_pReal,1.e-12_pReal)) call IO_error(153,ext_msg='constituent') if (dNeq(frac,1.0_pReal,1.e-12_pReal)) call IO_error(153,ext_msg='constituent')
enddo end do
allocate(material_O_0(materials%length)) allocate(material_O_0(materials%length))
@ -167,8 +167,8 @@ subroutine parse()
do co = 1, constituents%length do co = 1, constituents%length
constituent => constituents%get(co) constituent => constituents%get(co)
call material_O_0(ma)%data(co)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4)) call material_O_0(ma)%data(co)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4))
enddo end do
enddo end do
end subroutine parse end subroutine parse
@ -186,15 +186,15 @@ subroutine sanityCheck(materials,homogenizations)
constituents constituents
integer :: m integer :: m
if(maxval(discretization_materialAt) > materials%length) & if (maxval(discretization_materialAt) > materials%length) &
call IO_error(155,ext_msg='More materials requested than found in material.yaml') call IO_error(155,ext_msg='More materials requested than found in material.yaml')
do m = 1, materials%length do m = 1, materials%length
material => materials%get(m) material => materials%get(m)
constituents => material%get('constituents') constituents => material%get('constituents')
homogenization => homogenizations%get(material%get_asString('homogenization')) homogenization => homogenizations%get(material%get_asString('homogenization'))
if(constituents%length /= homogenization%get_asInt('N_constituents')) call IO_error(148) if (constituents%length /= homogenization%get_asInt('N_constituents')) call IO_error(148)
enddo end do
end subroutine sanityCheck end subroutine sanityCheck
@ -215,12 +215,12 @@ function getKeys(dict)
do i=1, dict%length do i=1, dict%length
temp(i) = dict%getKey(i) temp(i) = dict%getKey(i)
l = max(len_trim(temp(i)),l) l = max(len_trim(temp(i)),l)
enddo end do
allocate(character(l)::getKeys(dict%length)) allocate(character(l)::getKeys(dict%length))
do i=1, dict%length do i=1, dict%length
getKeys(i) = trim(temp(i)) getKeys(i) = trim(temp(i))
enddo end do
end function getKeys end function getKeys

8
src/math.f90
View File

@ -91,7 +91,7 @@ subroutine math_init
class(tNode), pointer :: & class(tNode), pointer :: &
num_generic num_generic
print'(/,a)', ' <<<+- math init -+>>>'; flush(IO_STDOUT) print'(/,1x,a)', '<<<+- math init -+>>>'; flush(IO_STDOUT)
num_generic => config_numerics%get('generic',defaultVal=emptyDict) num_generic => config_numerics%get('generic',defaultVal=emptyDict)
randomSeed = num_generic%get_asInt('random_seed', defaultVal = 0) randomSeed = num_generic%get_asInt('random_seed', defaultVal = 0)
@ -109,9 +109,9 @@ subroutine math_init
call random_seed(put = randInit) call random_seed(put = randInit)
call random_number(randTest) call random_number(randTest)
print'(a,i2)', ' size of random seed: ', randSize print'(/,a,i2)', ' size of random seed: ', randSize
print'(a,i0)', ' value of random seed: ', randInit(1) print'( a,i0)', ' value of random seed: ', randInit(1)
print'(a,4(/,26x,f17.14),/)', ' start of random sequence: ', randTest print'( a,4(/,26x,f17.14),/)', ' start of random sequence: ', randTest
call random_seed(put = randInit) call random_seed(put = randInit)

88
src/mesh/DAMASK_mesh.f90
View File

@ -86,7 +86,7 @@ program DAMASK_mesh
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init DAMASK (all modules) ! init DAMASK (all modules)
call CPFEM_initAll call CPFEM_initAll
print'(/,a)', ' <<<+- DAMASK_mesh init -+>>>'; flush(IO_STDOUT) print'(/,1x,a)', '<<<+- DAMASK_mesh init -+>>>'; flush(IO_STDOUT)
!--------------------------------------------------------------------- !---------------------------------------------------------------------
! reading field information from numerics file and do sanity checks ! reading field information from numerics file and do sanity checks
@ -115,16 +115,16 @@ program DAMASK_mesh
case('$Loadcase') case('$Loadcase')
N_def = N_def + 1 N_def = N_def + 1
end select end select
enddo ! count all identifiers to allocate memory and do sanity check end do ! count all identifiers to allocate memory and do sanity check
enddo end do
if(N_def < 1) call IO_error(error_ID = 837) if (N_def < 1) call IO_error(error_ID = 837)
allocate(loadCases(N_def)) allocate(loadCases(N_def))
do i = 1, size(loadCases) do i = 1, size(loadCases)
allocate(loadCases(i)%fieldBC(1)) allocate(loadCases(i)%fieldBC(1))
loadCases(i)%fieldBC(1)%ID = FIELD_MECH_ID loadCases(i)%fieldBC(1)%ID = FIELD_MECH_ID
enddo end do
do i = 1, size(loadCases) do i = 1, size(loadCases)
loadCases(i)%fieldBC(1)%nComponents = dimPlex !< X, Y (, Z) displacements loadCases(i)%fieldBC(1)%nComponents = dimPlex !< X, Y (, Z) displacements
@ -138,12 +138,12 @@ program DAMASK_mesh
case (3) case (3)
loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_Z_ID loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_Z_ID
end select end select
enddo end do
do component = 1, loadCases(i)%fieldBC(1)%nComponents do component = 1, loadCases(i)%fieldBC(1)%nComponents
allocate(loadCases(i)%fieldBC(1)%componentBC(component)%Value(mesh_Nboundaries), source = 0.0_pReal) allocate(loadCases(i)%fieldBC(1)%componentBC(component)%Value(mesh_Nboundaries), source = 0.0_pReal)
allocate(loadCases(i)%fieldBC(1)%componentBC(component)%Mask (mesh_Nboundaries), source = .false.) allocate(loadCases(i)%fieldBC(1)%componentBC(component)%Mask (mesh_Nboundaries), source = .false.)
enddo end do
enddo end do
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! reading the load case and assign values to the allocated data structure ! reading the load case and assign values to the allocated data structure
@ -163,7 +163,7 @@ program DAMASK_mesh
currentFaceSet = -1 currentFaceSet = -1
do faceSet = 1, mesh_Nboundaries do faceSet = 1, mesh_Nboundaries
if (mesh_boundaries(faceSet) == currentFace) currentFaceSet = faceSet if (mesh_boundaries(faceSet) == currentFace) currentFaceSet = faceSet
enddo end do
if (currentFaceSet < 0) call IO_error(error_ID = 837, ext_msg = 'invalid BC') if (currentFaceSet < 0) call IO_error(error_ID = 837, ext_msg = 'invalid BC')
case('t') case('t')
loadCases(currentLoadCase)%time = IO_floatValue(line,chunkPos,i+1) loadCases(currentLoadCase)%time = IO_floatValue(line,chunkPos,i+1)
@ -192,11 +192,11 @@ program DAMASK_mesh
.true. .true.
loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Value(currentFaceSet) = & loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1) IO_floatValue(line,chunkPos,i+1)
endif end if
enddo end do
end select end select
enddo end do
enddo end do
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! consistency checks and output of load case ! consistency checks and output of load case
@ -204,28 +204,28 @@ program DAMASK_mesh
errorID = 0 errorID = 0
checkLoadcases: do currentLoadCase = 1, size(loadCases) checkLoadcases: do currentLoadCase = 1, size(loadCases)
write (loadcase_string, '(i0)' ) currentLoadCase write (loadcase_string, '(i0)' ) currentLoadCase
print'(a,i0)', ' load case: ', currentLoadCase print'(/,1x,a,i0)', 'load case: ', currentLoadCase
if (.not. loadCases(currentLoadCase)%followFormerTrajectory) & if (.not. loadCases(currentLoadCase)%followFormerTrajectory) &
print'(a)', ' drop guessing along trajectory' print'(2x,a)', 'drop guessing along trajectory'
print'(a)', ' Field '//trim(FIELD_MECH_label) print'(2x,a)', 'Field '//trim(FIELD_MECH_label)
do faceSet = 1, mesh_Nboundaries do faceSet = 1, mesh_Nboundaries
do component = 1, loadCases(currentLoadCase)%fieldBC(1)%nComponents do component = 1, loadCases(currentLoadCase)%fieldBC(1)%nComponents
if (loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Mask(faceSet)) & if (loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Mask(faceSet)) &
print'(a,i2,a,i2,a,f12.7)', ' Face ', mesh_boundaries(faceSet), & print'(a,i2,a,i2,a,f12.7)', &
' Component ', component, & ' Face ', mesh_boundaries(faceSet), &
' Value ', loadCases(currentLoadCase)%fieldBC(1)% & ' Component ', component, &
componentBC(component)%Value(faceSet) ' Value ', loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Value(faceSet)
enddo end do
enddo end do
print'(a,f12.6)', ' time: ', loadCases(currentLoadCase)%time print'(2x,a,f12.6)', 'time: ', loadCases(currentLoadCase)%time
if (loadCases(currentLoadCase)%incs < 1) errorID = 835 ! non-positive incs count if (loadCases(currentLoadCase)%incs < 1) errorID = 835 ! non-positive incs count
print'(a,i5)', ' increments: ', loadCases(currentLoadCase)%incs print'(2x,a,i5)', 'increments: ', loadCases(currentLoadCase)%incs
if (loadCases(currentLoadCase)%outputfrequency < 1) errorID = 836 ! non-positive result frequency if (loadCases(currentLoadCase)%outputfrequency < 1) errorID = 836 ! non-positive result frequency
print'(a,i5)', ' output frequency: ', & print'(2x,a,i5)', 'output frequency: ', &
loadCases(currentLoadCase)%outputfrequency loadCases(currentLoadCase)%outputfrequency
if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message
enddo checkLoadcases end do checkLoadcases
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! doing initialization depending on active solvers ! doing initialization depending on active solvers
@ -235,9 +235,9 @@ program DAMASK_mesh
if (worldrank == 0) then if (worldrank == 0) then
open(newunit=statUnit,file=trim(getSolverJobName())//'.sta',form='FORMATTED',status='REPLACE') open(newunit=statUnit,file=trim(getSolverJobName())//'.sta',form='FORMATTED',status='REPLACE')
write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded' ! statistics file write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded' ! statistics file
endif end if
print'(/,a)', ' ... writing initial configuration to file ........................' print'(/,1x,a)', '... writing initial configuration to file ........................'
flush(IO_STDOUT) flush(IO_STDOUT)
call CPFEM_results(0,0.0_pReal) call CPFEM_results(0,0.0_pReal)
@ -262,7 +262,7 @@ program DAMASK_mesh
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report begin of new step ! report begin of new step
print'(/,a)', ' ###########################################################################' print'(/,1x,a)', '###########################################################################'
print'(1x,a,es12.5,6(a,i0))',& print'(1x,a,es12.5,6(a,i0))',&
'Time', time, & 'Time', time, &
's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,& 's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,&
@ -281,60 +281,60 @@ program DAMASK_mesh
stagIterate = .true. stagIterate = .true.
do while (stagIterate) do while (stagIterate)
solres(1) = FEM_mechanical_solution(incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(1)) solres(1) = FEM_mechanical_solution(incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(1))
if(.not. solres(1)%converged) exit if (.not. solres(1)%converged) exit
stagIter = stagIter + 1 stagIter = stagIter + 1
stagIterate = stagIter < stagItMax & stagIterate = stagIter < stagItMax &
.and. all(solres(:)%converged) & .and. all(solres(:)%converged) &
.and. .not. all(solres(:)%stagConverged) ! stationary with respect to staggered iteration .and. .not. all(solres(:)%stagConverged) ! stationary with respect to staggered iteration
enddo end do
! check solution ! check solution
cutBack = .False. cutBack = .False.
if(.not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found if (.not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found
if (cutBackLevel < maxCutBack) then ! do cut back if (cutBackLevel < maxCutBack) then ! do cut back
cutBack = .True. cutBack = .True.
stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator
cutBackLevel = cutBackLevel + 1 cutBackLevel = cutBackLevel + 1
time = time - timeinc ! rewind time time = time - timeinc ! rewind time
timeinc = timeinc/2.0_pReal timeinc = timeinc/2.0_pReal
print'(/,a)', ' cutting back' print'(/,1x,a)', 'cutting back'
else ! default behavior, exit if spectral solver does not converge else ! default behavior, exit if spectral solver does not converge
if (worldrank == 0) close(statUnit) if (worldrank == 0) close(statUnit)
call IO_error(950) call IO_error(950)
endif end if
else else
guess = .true. ! start guessing after first converged (sub)inc guess = .true. ! start guessing after first converged (sub)inc
timeIncOld = timeinc timeIncOld = timeinc
endif end if
if (.not. cutBack .and. worldrank == 0) & if (.not. cutBack .and. worldrank == 0) &
write(statUnit,*) totalIncsCounter, time, cutBackLevel, & write(statUnit,*) totalIncsCounter, time, cutBackLevel, &
solres%converged, solres%iterationsNeeded ! write statistics about accepted solution solres%converged, solres%iterationsNeeded ! write statistics about accepted solution
enddo subStepLooping end do subStepLooping
cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc
if (all(solres(:)%converged)) then if (all(solres(:)%converged)) then
print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' converged' print'(/,1x,a,i0,a)', 'increment ', totalIncsCounter, ' converged'
else else
print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' NOT converged' print'(/,1x,a,i0,a)', 'increment ', totalIncsCounter, ' NOT converged'
endif; flush(IO_STDOUT) end if; flush(IO_STDOUT)
if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then ! at output frequency if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then ! at output frequency
print'(/,a)', ' ... writing results to file ......................................' print'(/,1x,a)', '... writing results to file ......................................'
call FEM_mechanical_updateCoords call FEM_mechanical_updateCoords
call CPFEM_results(totalIncsCounter,time) call CPFEM_results(totalIncsCounter,time)
endif end if
enddo incLooping end do incLooping
enddo loadCaseLooping end do loadCaseLooping
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report summary of whole calculation ! report summary of whole calculation
print'(/,a)', ' ###########################################################################' print'(/,1x,a)', '###########################################################################'
if (worldrank == 0) close(statUnit) if (worldrank == 0) close(statUnit)
call quit(0) ! no complains ;) call quit(0) ! no complains ;)

6
src/mesh/FEM_quadrature.f90
View File

@ -41,10 +41,10 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_quadrature_init() subroutine FEM_quadrature_init()
print'(/,a)', ' <<<+- FEM_quadrature init -+>>>'; flush(6) print'(/,1x,a)', '<<<+- FEM_quadrature init -+>>>'; flush(6)
print*, 'L. Zhang et al., Journal of Computational Mathematics 27(1):89-96, 2009' print'(/,1x,a)', 'L. Zhang et al., Journal of Computational Mathematics 27(1):89-96, 2009'
print*, 'https://www.jstor.org/stable/43693493' print'( 1x,a)', 'https://www.jstor.org/stable/43693493'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 2D linear ! 2D linear

10
src/mesh/FEM_utilities.f90
View File

@ -96,7 +96,7 @@ subroutine FEM_utilities_init
logical :: debugPETSc !< use some in debug defined options for more verbose PETSc solution logical :: debugPETSc !< use some in debug defined options for more verbose PETSc solution
print'(/,a)', ' <<<+- FEM_utilities init -+>>>' print'(/,1x,a)', '<<<+- FEM_utilities init -+>>>'
num_mesh => config_numerics%get('mesh',defaultVal=emptyDict) num_mesh => config_numerics%get('mesh',defaultVal=emptyDict)
@ -111,10 +111,10 @@ subroutine FEM_utilities_init
debug_mesh => config_debug%get('mesh',defaultVal=emptyList) debug_mesh => config_debug%get('mesh',defaultVal=emptyList)
debugPETSc = debug_mesh%contains('PETSc') debugPETSc = debug_mesh%contains('PETSc')
if(debugPETSc) print'(3(/,a),/)', & if(debugPETSc) print'(3(/,1x,a),/)', &
' Initializing PETSc with debug options: ', & 'Initializing PETSc with debug options: ', &
trim(PETScDebug), & trim(PETScDebug), &
' add more using the "PETSc_options" keyword in numerics.yaml' 'add more using the "PETSc_options" keyword in numerics.yaml'
flush(IO_STDOUT) flush(IO_STDOUT)
call PetscOptionsClear(PETSC_NULL_OPTIONS,ierr) call PetscOptionsClear(PETSC_NULL_OPTIONS,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
@ -151,7 +151,7 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
PetscErrorCode :: ierr PetscErrorCode :: ierr
print'(/,a)', ' ... evaluating constitutive response ......................................' print'(/,1x,a)', '... evaluating constitutive response ......................................'
call homogenization_mechanical_response(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) ! calculate P field call homogenization_mechanical_response(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) ! calculate P field
if (.not. terminallyIll) & if (.not. terminallyIll) &

3
src/mesh/discretization_mesh.f90
View File

@ -88,7 +88,7 @@ subroutine discretization_mesh_init(restart)
integer :: p_i !< integration order (quadrature rule) integer :: p_i !< integration order (quadrature rule)
type(tvec) :: coords_node0 type(tvec) :: coords_node0
print'(/,a)', ' <<<+- discretization_mesh init -+>>>' print'(/,1x,a)', '<<<+- discretization_mesh init -+>>>'
!-------------------------------------------------------------------------------- !--------------------------------------------------------------------------------
! read numerics parameter ! read numerics parameter
@ -106,6 +106,7 @@ subroutine discretization_mesh_init(restart)
CHKERRQ(ierr) CHKERRQ(ierr)
call DMGetStratumSize(globalMesh,'depth',dimPlex,mesh_NcpElemsGlobal,ierr) call DMGetStratumSize(globalMesh,'depth',dimPlex,mesh_NcpElemsGlobal,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
print'()'
call DMView(globalMesh, PETSC_VIEWER_STDOUT_WORLD,ierr) call DMView(globalMesh, PETSC_VIEWER_STDOUT_WORLD,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)

8
src/mesh/mesh_mech_FEM.f90
View File

@ -113,7 +113,7 @@ subroutine FEM_mechanical_init(fieldBC)
class(tNode), pointer :: & class(tNode), pointer :: &
num_mesh num_mesh
print'(/,a)', ' <<<+- FEM_mech init -+>>>'; flush(IO_STDOUT) print'(/,1x,a)', '<<<+- FEM_mech init -+>>>'; flush(IO_STDOUT)
!----------------------------------------------------------------------------- !-----------------------------------------------------------------------------
! read numerical parametes and do sanity checks ! read numerical parametes and do sanity checks
@ -302,7 +302,7 @@ type(tSolutionState) function FEM_mechanical_solution( &
CHKERRQ(ierr) CHKERRQ(ierr)
endif endif
print'(/,a)', ' ===========================================================================' print'(/,1x,a)', '==========================================================================='
flush(IO_STDOUT) flush(IO_STDOUT)
end function FEM_mechanical_solution end function FEM_mechanical_solution
@ -663,8 +663,8 @@ subroutine FEM_mechanical_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reaso
print'(/,1x,a,a,i0,a,i0,f0.3)', trim(incInfo), & print'(/,1x,a,a,i0,a,i0,f0.3)', trim(incInfo), &
' @ Iteration ',PETScIter,' mechanical residual norm = ', & ' @ Iteration ',PETScIter,' mechanical residual norm = ', &
int(fnorm/divTol),fnorm/divTol-int(fnorm/divTol) int(fnorm/divTol),fnorm/divTol-int(fnorm/divTol)
print'(/,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', & print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &
' Piola--Kirchhoff stress / MPa =',transpose(P_av)*1.e-6_pReal 'Piola--Kirchhoff stress / MPa =',transpose(P_av)*1.e-6_pReal
flush(IO_STDOUT) flush(IO_STDOUT)
end subroutine FEM_mechanical_converged end subroutine FEM_mechanical_converged

10
src/parallelization.f90
View File

@ -76,11 +76,11 @@ subroutine parallelization_init
call MPI_Comm_rank(MPI_COMM_WORLD,worldrank,err) call MPI_Comm_rank(MPI_COMM_WORLD,worldrank,err)
if (err /= 0) error stop 'Could not determine worldrank' if (err /= 0) error stop 'Could not determine worldrank'
if (worldrank == 0) print'(/,a)', ' <<<+- parallelization init -+>>>' if (worldrank == 0) print'(/,1x,a)', '<<<+- parallelization init -+>>>'
call MPI_Comm_size(MPI_COMM_WORLD,worldsize,err) call MPI_Comm_size(MPI_COMM_WORLD,worldsize,err)
if (err /= 0) error stop 'Could not determine worldsize' if (err /= 0) error stop 'Could not determine worldsize'
if (worldrank == 0) print'(a,i3)', ' MPI processes: ',worldsize if (worldrank == 0) print'(/,1x,a,i3)', 'MPI processes: ',worldsize
call MPI_Type_size(MPI_INTEGER,typeSize,err) call MPI_Type_size(MPI_INTEGER,typeSize,err)
if (err /= 0) error stop 'Could not determine MPI integer size' if (err /= 0) error stop 'Could not determine MPI integer size'
@ -97,16 +97,16 @@ subroutine parallelization_init
!$ call get_environment_variable(name='OMP_NUM_THREADS',value=NumThreadsString,STATUS=got_env) !$ call get_environment_variable(name='OMP_NUM_THREADS',value=NumThreadsString,STATUS=got_env)
!$ if(got_env /= 0) then !$ if(got_env /= 0) then
!$ print*, 'Could not get $OMP_NUM_THREADS, using default' !$ print'(1x,a)', 'Could not get $OMP_NUM_THREADS, using default'
!$ OMP_NUM_THREADS = 4_pI32 !$ OMP_NUM_THREADS = 4_pI32
!$ else !$ else
!$ read(NumThreadsString,'(i6)') OMP_NUM_THREADS !$ read(NumThreadsString,'(i6)') OMP_NUM_THREADS
!$ if (OMP_NUM_THREADS < 1_pI32) then !$ if (OMP_NUM_THREADS < 1_pI32) then
!$ print*, 'Invalid OMP_NUM_THREADS: "'//trim(NumThreadsString)//'", using default' !$ print'(1x,a)', 'Invalid OMP_NUM_THREADS: "'//trim(NumThreadsString)//'", using default'
!$ OMP_NUM_THREADS = 4_pI32 !$ OMP_NUM_THREADS = 4_pI32
!$ endif !$ endif
!$ endif !$ endif
!$ print'(a,i2)', ' OMP_NUM_THREADS: ',OMP_NUM_THREADS !$ print'(1x,a,1x,i2)', 'OMP_NUM_THREADS:',OMP_NUM_THREADS
!$ call omp_set_num_threads(OMP_NUM_THREADS) !$ call omp_set_num_threads(OMP_NUM_THREADS)
end subroutine parallelization_init end subroutine parallelization_init

52
src/phase.f90
View File

@ -343,7 +343,7 @@ subroutine phase_init
phase phase
print'(/,a)', ' <<<+- phase init -+>>>'; flush(IO_STDOUT) print'(/,1x,a)', '<<<+- phase init -+>>>'; flush(IO_STDOUT)
debug_constitutive => config_debug%get('phase', defaultVal=emptyList) debug_constitutive => config_debug%get('phase', defaultVal=emptyList)
debugConstitutive%basic = debug_constitutive%contains('basic') debugConstitutive%basic = debug_constitutive%contains('basic')
@ -371,20 +371,20 @@ subroutine phase_init
phase_cOverA(ph) = phase%get_asFloat('c/a') phase_cOverA(ph) = phase%get_asFloat('c/a')
phase_rho(ph) = phase%get_asFloat('rho',defaultVal=0.0_pReal) phase_rho(ph) = phase%get_asFloat('rho',defaultVal=0.0_pReal)
allocate(phase_O_0(ph)%data(count(material_phaseID==ph))) allocate(phase_O_0(ph)%data(count(material_phaseID==ph)))
enddo end do
do ce = 1, size(material_phaseID,2) do ce = 1, size(material_phaseID,2)
ma = discretization_materialAt((ce-1)/discretization_nIPs+1) ma = discretization_materialAt((ce-1)/discretization_nIPs+1)
do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
ph = material_phaseID(co,ce) ph = material_phaseID(co,ce)
phase_O_0(ph)%data(material_phaseEntry(co,ce)) = material_O_0(ma)%data(co) phase_O_0(ph)%data(material_phaseEntry(co,ce)) = material_O_0(ma)%data(co)
enddo end do
enddo end do
allocate(phase_O(phases%length)) allocate(phase_O(phases%length))
do ph = 1,phases%length do ph = 1,phases%length
phase_O(ph)%data = phase_O_0(ph)%data phase_O(ph)%data = phase_O_0(ph)%data
enddo end do
call mechanical_init(phases) call mechanical_init(phases)
call damage_init call damage_init
@ -471,7 +471,7 @@ subroutine phase_results()
call mechanical_results(group,ph) call mechanical_results(group,ph)
call damage_results(group,ph) call damage_results(group,ph)
enddo end do
end subroutine phase_results end subroutine phase_results
@ -495,7 +495,7 @@ subroutine crystallite_init()
phases phases
print'(/,a)', ' <<<+- crystallite init -+>>>' print'(/,1x,a)', '<<<+- crystallite init -+>>>'
cMax = homogenization_maxNconstituents cMax = homogenization_maxNconstituents
iMax = discretization_nIPs iMax = discretization_nIPs
@ -515,28 +515,28 @@ subroutine crystallite_init()
num%nState = num_crystallite%get_asInt ('nState', defaultVal=20) num%nState = num_crystallite%get_asInt ('nState', defaultVal=20)
num%nStress = num_crystallite%get_asInt ('nStress', defaultVal=40) num%nStress = num_crystallite%get_asInt ('nStress', defaultVal=40)
if(num%subStepMinCryst <= 0.0_pReal) call IO_error(301,ext_msg='subStepMinCryst') if (num%subStepMinCryst <= 0.0_pReal) call IO_error(301,ext_msg='subStepMinCryst')
if(num%subStepSizeCryst <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeCryst') if (num%subStepSizeCryst <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeCryst')
if(num%stepIncreaseCryst <= 0.0_pReal) call IO_error(301,ext_msg='stepIncreaseCryst') if (num%stepIncreaseCryst <= 0.0_pReal) call IO_error(301,ext_msg='stepIncreaseCryst')
if(num%subStepSizeLp <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeLp') if (num%subStepSizeLp <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeLp')
if(num%subStepSizeLi <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeLi') if (num%subStepSizeLi <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeLi')
if(num%rtol_crystalliteState <= 0.0_pReal) call IO_error(301,ext_msg='rtol_crystalliteState') if (num%rtol_crystalliteState <= 0.0_pReal) call IO_error(301,ext_msg='rtol_crystalliteState')
if(num%rtol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='rtol_crystalliteStress') if (num%rtol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='rtol_crystalliteStress')
if(num%atol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='atol_crystalliteStress') if (num%atol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='atol_crystalliteStress')
if(num%iJacoLpresiduum < 1) call IO_error(301,ext_msg='iJacoLpresiduum') if (num%iJacoLpresiduum < 1) call IO_error(301,ext_msg='iJacoLpresiduum')
if(num%nState < 1) call IO_error(301,ext_msg='nState') if (num%nState < 1) call IO_error(301,ext_msg='nState')
if(num%nStress< 1) call IO_error(301,ext_msg='nStress') if (num%nStress< 1) call IO_error(301,ext_msg='nStress')
phases => config_material%get('phase') phases => config_material%get('phase')
print'(a42,1x,i10)', ' # of elements: ', eMax print'(/,a42,1x,i10)', ' # of elements: ', eMax
print'(a42,1x,i10)', ' # of integration points/element: ', iMax print'( a42,1x,i10)', ' # of integration points/element: ', iMax
print'(a42,1x,i10)', 'max # of constituents/integration point: ', cMax print'( a42,1x,i10)', 'max # of constituents/integration point: ', cMax
flush(IO_STDOUT) flush(IO_STDOUT)
@ -547,9 +547,9 @@ subroutine crystallite_init()
do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
call crystallite_orientations(co,ip,el) call crystallite_orientations(co,ip,el)
call plastic_dependentState(co,ip,el) ! update dependent state variables to be consistent with basic states call plastic_dependentState(co,ip,el) ! update dependent state variables to be consistent with basic states
enddo end do
enddo end do
enddo end do
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
@ -642,7 +642,7 @@ subroutine phase_restartWrite(fileHandle)
call HDF5_closeGroup(groupHandle(2)) call HDF5_closeGroup(groupHandle(2))
enddo end do
call HDF5_closeGroup(groupHandle(1)) call HDF5_closeGroup(groupHandle(1))
@ -670,7 +670,7 @@ subroutine phase_restartRead(fileHandle)
call HDF5_closeGroup(groupHandle(2)) call HDF5_closeGroup(groupHandle(2))
enddo end do
call HDF5_closeGroup(groupHandle(1)) call HDF5_closeGroup(groupHandle(1))

36
src/phase_damage.f90
View File

@ -83,7 +83,7 @@ module subroutine damage_init
source source
logical:: damage_active logical:: damage_active
print'(/,a)', ' <<<+- phase:damage init -+>>>' print'(/,1x,a)', '<<<+- phase:damage init -+>>>'
phases => config_material%get('phase') phases => config_material%get('phase')
@ -108,16 +108,16 @@ module subroutine damage_init
param(ph)%D(1,1) = source%get_asFloat('D_11') param(ph)%D(1,1) = source%get_asFloat('D_11')
if (any(phase_lattice(ph) == ['hP','tI'])) param(ph)%D(3,3) = source%get_asFloat('D_33') if (any(phase_lattice(ph) == ['hP','tI'])) param(ph)%D(3,3) = source%get_asFloat('D_33')
param(ph)%D = lattice_symmetrize_33(param(ph)%D,phase_lattice(ph)) param(ph)%D = lattice_symmetrize_33(param(ph)%D,phase_lattice(ph))
endif end if
enddo end do
allocate(phase_damage(phases%length), source = DAMAGE_UNDEFINED_ID) allocate(phase_damage(phases%length), source = DAMAGE_UNDEFINED_ID)
if (damage_active) then if (damage_active) then
where(isobrittle_init() ) phase_damage = DAMAGE_ISOBRITTLE_ID where(isobrittle_init() ) phase_damage = DAMAGE_ISOBRITTLE_ID
where(anisobrittle_init()) phase_damage = DAMAGE_ANISOBRITTLE_ID where(anisobrittle_init()) phase_damage = DAMAGE_ANISOBRITTLE_ID
endif end if
phase_damage_maxSizeDotState = maxval(damageState%sizeDotState) phase_damage_maxSizeDotState = maxval(damageState%sizeDotState)
@ -181,7 +181,7 @@ module subroutine damage_restore(ce)
if (damageState(material_phaseID(co,ce))%sizeState > 0) & if (damageState(material_phaseID(co,ce))%sizeState > 0) &
damageState(material_phaseID(co,ce))%state( :,material_phaseEntry(co,ce)) = & damageState(material_phaseID(co,ce))%state( :,material_phaseEntry(co,ce)) = &
damageState(material_phaseID(co,ce))%state0(:,material_phaseEntry(co,ce)) damageState(material_phaseID(co,ce))%state0(:,material_phaseEntry(co,ce))
enddo end do
end subroutine damage_restore end subroutine damage_restore
@ -242,11 +242,11 @@ function integrateDamageState(Delta_t,ph,en) result(broken)
if (damageState(ph)%sizeState == 0) then if (damageState(ph)%sizeState == 0) then
broken = .false. broken = .false.
return return
endif end if
converged_ = .true. converged_ = .true.
broken = phase_damage_collectDotState(ph,en) broken = phase_damage_collectDotState(ph,en)
if(broken) return if (broken) return
size_so = damageState(ph)%sizeDotState size_so = damageState(ph)%sizeDotState
damageState(ph)%state(1:size_so,en) = damageState(ph)%state0 (1:size_so,en) & damageState(ph)%state(1:size_so,en) = damageState(ph)%state0 (1:size_so,en) &
@ -255,11 +255,11 @@ function integrateDamageState(Delta_t,ph,en) result(broken)
iteration: do NiterationState = 1, num%nState iteration: do NiterationState = 1, num%nState
if(nIterationState > 1) source_dotState(1:size_so,2) = source_dotState(1:size_so,1) if (nIterationState > 1) source_dotState(1:size_so,2) = source_dotState(1:size_so,1)
source_dotState(1:size_so,1) = damageState(ph)%dotState(:,en) source_dotState(1:size_so,1) = damageState(ph)%dotState(:,en)
broken = phase_damage_collectDotState(ph,en) broken = phase_damage_collectDotState(ph,en)
if(broken) exit iteration if (broken) exit iteration
zeta = damper(damageState(ph)%dotState(:,en),source_dotState(1:size_so,1),source_dotState(1:size_so,2)) zeta = damper(damageState(ph)%dotState(:,en),source_dotState(1:size_so,1),source_dotState(1:size_so,2))
@ -274,12 +274,12 @@ function integrateDamageState(Delta_t,ph,en) result(broken)
damageState(ph)%atol(1:size_so)) damageState(ph)%atol(1:size_so))
if(converged_) then if (converged_) then
broken = phase_damage_deltaState(mechanical_F_e(ph,en),ph,en) broken = phase_damage_deltaState(mechanical_F_e(ph,en),ph,en)
exit iteration exit iteration
endif end if
enddo iteration end do iteration
broken = broken .or. .not. converged_ broken = broken .or. .not. converged_
@ -302,7 +302,7 @@ function integrateDamageState(Delta_t,ph,en) result(broken)
damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22)
else else
damper = 1.0_pReal damper = 1.0_pReal
endif end if
end function damper end function damper
@ -358,7 +358,7 @@ function phase_damage_collectDotState(ph,en) result(broken)
broken = broken .or. any(IEEE_is_NaN(damageState(ph)%dotState(:,en))) broken = broken .or. any(IEEE_is_NaN(damageState(ph)%dotState(:,en)))
endif end if
end function phase_damage_collectDotState end function phase_damage_collectDotState
@ -419,12 +419,12 @@ function phase_damage_deltaState(Fe, ph, en) result(broken)
case (DAMAGE_ISOBRITTLE_ID) sourceType case (DAMAGE_ISOBRITTLE_ID) sourceType
call isobrittle_deltaState(phase_homogenizedC(ph,en), Fe, ph,en) call isobrittle_deltaState(phase_homogenizedC(ph,en), Fe, ph,en)
broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,en))) broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,en)))
if(.not. broken) then if (.not. broken) then
myOffset = damageState(ph)%offsetDeltaState myOffset = damageState(ph)%offsetDeltaState
mySize = damageState(ph)%sizeDeltaState mySize = damageState(ph)%sizeDeltaState
damageState(ph)%state(myOffset + 1: myOffset + mySize,en) = & damageState(ph)%state(myOffset + 1: myOffset + mySize,en) = &
damageState(ph)%state(myOffset + 1: myOffset + mySize,en) + damageState(ph)%deltaState(1:mySize,en) damageState(ph)%state(myOffset + 1: myOffset + mySize,en) + damageState(ph)%deltaState(1:mySize,en)
endif end if
end select sourceType end select sourceType
@ -454,7 +454,7 @@ function source_active(source_label) result(active_source)
sources => phase%get('damage',defaultVal=emptyList) sources => phase%get('damage',defaultVal=emptyList)
src => sources%get(1) src => sources%get(1)
active_source(ph) = src%get_asString('type',defaultVal = 'x') == source_label active_source(ph) = src%get_asString('type',defaultVal = 'x') == source_label
enddo end do
end function source_active end function source_active
@ -497,7 +497,7 @@ module subroutine damage_forward()
do ph = 1, size(damageState) do ph = 1, size(damageState)
if (damageState(ph)%sizeState > 0) & if (damageState(ph)%sizeState > 0) &
damageState(ph)%state0 = damageState(ph)%state damageState(ph)%state0 = damageState(ph)%state
enddo end do
end subroutine damage_forward end subroutine damage_forward

26
src/phase_damage_anisobrittle.f90
View File

@ -46,10 +46,10 @@ module function anisobrittle_init() result(mySources)
mySources = source_active('anisobrittle') mySources = source_active('anisobrittle')
if(count(mySources) == 0) return if (count(mySources) == 0) return
print'(/,a)', ' <<<+- phase:damage:anisobrittle init -+>>>' print'(/,1x,a)', '<<<+- phase:damage:anisobrittle init -+>>>'
print'(a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT) print'(/,a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT)
phases => config_material%get('phase') phases => config_material%get('phase')
@ -57,7 +57,7 @@ module function anisobrittle_init() result(mySources)
do ph = 1, phases%length do ph = 1, phases%length
if(mySources(ph)) then if (mySources(ph)) then
phase => phases%get(ph) phase => phases%get(ph)
sources => phase%get('damage') sources => phase%get('damage')
@ -94,14 +94,14 @@ module function anisobrittle_init() result(mySources)
Nmembers = count(material_phaseID==ph) Nmembers = count(material_phaseID==ph)
call phase_allocateState(damageState(ph),Nmembers,1,1,0) call phase_allocateState(damageState(ph),Nmembers,1,1,0)
damageState(ph)%atol = src%get_asFloat('atol_phi',defaultVal=1.0e-9_pReal) damageState(ph)%atol = src%get_asFloat('atol_phi',defaultVal=1.0e-9_pReal)
if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi' if (any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi'
end associate end associate
if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_anisoBrittle)') if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_anisoBrittle)')
endif end if
enddo end do
end function anisobrittle_init end function anisobrittle_init
@ -136,7 +136,7 @@ module subroutine anisobrittle_dotState(S, ph,en)
* ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%q + & * ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%q + &
(max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%q + & (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%q + &
(max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**prm%q) (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**prm%q)
enddo end do
end associate end associate
end subroutine anisobrittle_dotState end subroutine anisobrittle_dotState
@ -159,7 +159,7 @@ module subroutine anisobrittle_results(phase,group)
case ('f_phi') case ('f_phi')
call results_writeDataset(stt,group,trim(prm%output(o)),'driving force','J/m³') call results_writeDataset(stt,group,trim(prm%output(o)),'driving force','J/m³')
end select end select
enddo outputsLoop end do outputsLoop
end associate end associate
end subroutine anisobrittle_results end subroutine anisobrittle_results
@ -200,7 +200,7 @@ module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,en)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) & dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) &
+ dudotd_dt*prm%cleavage_systems(k,l,1,i) * prm%cleavage_systems(m,n,1,i) + dudotd_dt*prm%cleavage_systems(k,l,1,i) * prm%cleavage_systems(m,n,1,i)
endif end if
traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i)) traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i))
if (abs(traction_t) > traction_crit + tol_math_check) then if (abs(traction_t) > traction_crit + tol_math_check) then
@ -210,7 +210,7 @@ module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,en)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) & dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) &
+ dudott_dt*prm%cleavage_systems(k,l,2,i) * prm%cleavage_systems(m,n,2,i) + dudott_dt*prm%cleavage_systems(k,l,2,i) * prm%cleavage_systems(m,n,2,i)
endif end if
traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i)) traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i))
if (abs(traction_n) > traction_crit + tol_math_check) then if (abs(traction_n) > traction_crit + tol_math_check) then
@ -220,8 +220,8 @@ module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,en)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) & dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) &
+ dudotn_dt*prm%cleavage_systems(k,l,3,i) * prm%cleavage_systems(m,n,3,i) + dudotn_dt*prm%cleavage_systems(k,l,3,i) * prm%cleavage_systems(m,n,3,i)
endif end if
enddo end do
end associate end associate
end subroutine damage_anisobrittle_LiAndItsTangent end subroutine damage_anisobrittle_LiAndItsTangent

16
src/phase_damage_isobrittle.f90
View File

@ -44,10 +44,10 @@ module function isobrittle_init() result(mySources)
mySources = source_active('isobrittle') mySources = source_active('isobrittle')
if(count(mySources) == 0) return if (count(mySources) == 0) return
print'(/,a)', ' <<<+- phase:damage:isobrittle init -+>>>' print'(/,1x,a)', '<<<+- phase:damage:isobrittle init -+>>>'
print'(a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT) print'(/,a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT)
phases => config_material%get('phase') phases => config_material%get('phase')
@ -56,7 +56,7 @@ module function isobrittle_init() result(mySources)
allocate(deltaState(phases%length)) allocate(deltaState(phases%length))
do ph = 1, phases%length do ph = 1, phases%length
if(mySources(ph)) then if (mySources(ph)) then
phase => phases%get(ph) phase => phases%get(ph)
sources => phase%get('damage') sources => phase%get('damage')
@ -77,7 +77,7 @@ module function isobrittle_init() result(mySources)
Nmembers = count(material_phaseID==ph) Nmembers = count(material_phaseID==ph)
call phase_allocateState(damageState(ph),Nmembers,1,1,1) call phase_allocateState(damageState(ph),Nmembers,1,1,1)
damageState(ph)%atol = src%get_asFloat('atol_phi',defaultVal=1.0e-9_pReal) damageState(ph)%atol = src%get_asFloat('atol_phi',defaultVal=1.0e-9_pReal)
if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi' if (any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi'
stt%r_W => damageState(ph)%state(1,:) stt%r_W => damageState(ph)%state(1,:)
dlt%r_W => damageState(ph)%deltaState(1,:) dlt%r_W => damageState(ph)%deltaState(1,:)
@ -86,9 +86,9 @@ module function isobrittle_init() result(mySources)
if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isobrittle)') if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isobrittle)')
endif end if
enddo end do
end function isobrittle_init end function isobrittle_init
@ -141,7 +141,7 @@ module subroutine isobrittle_results(phase,group)
case ('f_phi') case ('f_phi')
call results_writeDataset(stt,group,trim(prm%output(o)),'driving force','J/m³') ! Wrong, this is dimensionless call results_writeDataset(stt,group,trim(prm%output(o)),'driving force','J/m³') ! Wrong, this is dimensionless
end select end select
enddo outputsLoop end do outputsLoop
end associate end associate

2
src/phase_mechanical.f90
View File

@ -214,7 +214,7 @@ module subroutine mechanical_init(phases)
phase, & phase, &
mech mech
print'(/,a)', ' <<<+- phase:mechanical init -+>>>' print'(/,1x,a)', '<<<+- phase:mechanical init -+>>>'
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
allocate(output_constituent(phases%length)) allocate(output_constituent(phases%length))

20
src/phase_mechanical_eigen.f90
View File

@ -44,7 +44,7 @@ module subroutine eigen_init(phases)
kinematics, & kinematics, &
mechanics mechanics
print'(/,a)', ' <<<+- phase:mechanical:eigen init -+>>>' print'(/,1x,a)', '<<<+- phase:mechanical:eigen init -+>>>'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! explicit eigen mechanisms ! explicit eigen mechanisms
@ -55,11 +55,11 @@ module subroutine eigen_init(phases)
mechanics => phase%get('mechanical') mechanics => phase%get('mechanical')
kinematics => mechanics%get('eigen',defaultVal=emptyList) kinematics => mechanics%get('eigen',defaultVal=emptyList)
Nmodels(ph) = kinematics%length Nmodels(ph) = kinematics%length
enddo end do
allocate(model(maxval(Nmodels),phases%length), source = KINEMATICS_undefined_ID) allocate(model(maxval(Nmodels),phases%length), source = KINEMATICS_undefined_ID)
if(maxval(Nmodels) /= 0) then if (maxval(Nmodels) /= 0) then
where(thermalexpansion_init(maxval(Nmodels))) model = KINEMATICS_thermal_expansion_ID where(thermalexpansion_init(maxval(Nmodels))) model = KINEMATICS_thermal_expansion_ID
endif endif
@ -97,8 +97,8 @@ function kinematics_active(kinematics_label,kinematics_length) result(active_ki
do k = 1, kinematics%length do k = 1, kinematics%length
kinematics_type => kinematics%get(k) kinematics_type => kinematics%get(k)
active_kinematics(k,ph) = kinematics_type%get_asString('type') == kinematics_label active_kinematics(k,ph) = kinematics_type%get_asString('type') == kinematics_label
enddo end do
enddo end do
end function kinematics_active end function kinematics_active
@ -125,11 +125,11 @@ function kinematics_active2(kinematics_label) result(active_kinematics)
do ph = 1, phases%length do ph = 1, phases%length
phase => phases%get(ph) phase => phases%get(ph)
kinematics => phase%get('damage',defaultVal=emptyList) kinematics => phase%get('damage',defaultVal=emptyList)
if(kinematics%length < 1) return if (kinematics%length < 1) return
kinematics_type => kinematics%get(1) kinematics_type => kinematics%get(1)
if (.not. kinematics_type%contains('type')) continue if (.not. kinematics_type%contains('type')) continue
active_kinematics(ph) = kinematics_type%get_asString('type',defaultVal='n/a') == kinematics_label active_kinematics(ph) = kinematics_type%get_asString('type',defaultVal='n/a') == kinematics_label
enddo end do
end function kinematics_active2 end function kinematics_active2
@ -188,7 +188,7 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
dLi_dS = dLi_dS + my_dLi_dS dLi_dS = dLi_dS + my_dLi_dS
active = .true. active = .true.
end select kinematicsType end select kinematicsType
enddo KinematicsLoop end do KinematicsLoop
select case (model_damage(ph)) select case (model_damage(ph))
case (KINEMATICS_cleavage_opening_ID) case (KINEMATICS_cleavage_opening_ID)
@ -198,7 +198,7 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
active = .true. active = .true.
end select end select
if(.not. active) return if (.not. active) return
FiInv = math_inv33(Fi) FiInv = math_inv33(Fi)
detFi = math_det33(Fi) detFi = math_det33(Fi)
@ -209,7 +209,7 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
dLi_dS(1:3,1:3,i,j) = matmul(matmul(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi dLi_dS(1:3,1:3,i,j) = matmul(matmul(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi
dLi_dFi(1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i) dLi_dFi(1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i)
dLi_dFi(1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i) dLi_dFi(1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i)
enddo; enddo end do; end do
end subroutine phase_LiAndItsTangents end subroutine phase_LiAndItsTangents

4
src/phase_mechanical_eigen_cleavageopening.f90
View File

@ -21,8 +21,8 @@ module function damage_anisobrittle_init() result(myKinematics)
myKinematics = kinematics_active2('anisobrittle') myKinematics = kinematics_active2('anisobrittle')
if(count(myKinematics) == 0) return if(count(myKinematics) == 0) return
print'(/,a)', ' <<<+- phase:mechanical:eigen:cleavageopening init -+>>>' print'(/,1x,a)', '<<<+- phase:mechanical:eigen:cleavageopening init -+>>>'
print'(a,i2)', ' # phases: ',count(myKinematics); flush(IO_STDOUT) print'(/,a,i2)', ' # phases: ',count(myKinematics); flush(IO_STDOUT)
end function damage_anisobrittle_init end function damage_anisobrittle_init

22
src/phase_mechanical_eigen_thermalexpansion.f90
View File

@ -36,25 +36,25 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics)
kinematics, & kinematics, &
kinematic_type kinematic_type
print'(/,a)', ' <<<+- phase:mechanical:eigen:thermalexpansion init -+>>>' print'(/,1x,a)', '<<<+- phase:mechanical:eigen:thermalexpansion init -+>>>'
myKinematics = kinematics_active('thermalexpansion',kinematics_length) myKinematics = kinematics_active('thermalexpansion',kinematics_length)
Ninstances = count(myKinematics) Ninstances = count(myKinematics)
print'(a,i2)', ' # phases: ',Ninstances; flush(IO_STDOUT) print'(/,a,i2)', ' # phases: ',Ninstances; flush(IO_STDOUT)
if(Ninstances == 0) return if (Ninstances == 0) return
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(param(Ninstances)) allocate(param(Ninstances))
allocate(kinematics_thermal_expansion_instance(phases%length), source=0) allocate(kinematics_thermal_expansion_instance(phases%length), source=0)
do p = 1, phases%length do p = 1, phases%length
if(any(myKinematics(:,p))) kinematics_thermal_expansion_instance(p) = count(myKinematics(:,1:p)) if (any(myKinematics(:,p))) kinematics_thermal_expansion_instance(p) = count(myKinematics(:,1:p))
phase => phases%get(p) phase => phases%get(p)
if(count(myKinematics(:,p)) == 0) cycle if (count(myKinematics(:,p)) == 0) cycle
mech => phase%get('mechanical') mech => phase%get('mechanical')
kinematics => mech%get('eigen') kinematics => mech%get('eigen')
do k = 1, kinematics%length do k = 1, kinematics%length
if(myKinematics(k,p)) then if (myKinematics(k,p)) then
associate(prm => param(kinematics_thermal_expansion_instance(p))) associate(prm => param(kinematics_thermal_expansion_instance(p)))
kinematic_type => kinematics%get(k) kinematic_type => kinematics%get(k)
@ -67,14 +67,14 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics)
prm%A(3,3,1) = kinematic_type%get_asFloat('A_33') prm%A(3,3,1) = kinematic_type%get_asFloat('A_33')
prm%A(3,3,2) = kinematic_type%get_asFloat('A_33,T',defaultVal=0.0_pReal) prm%A(3,3,2) = kinematic_type%get_asFloat('A_33,T',defaultVal=0.0_pReal)
prm%A(3,3,3) = kinematic_type%get_asFloat('A_33,T^2',defaultVal=0.0_pReal) prm%A(3,3,3) = kinematic_type%get_asFloat('A_33,T^2',defaultVal=0.0_pReal)
endif end if
do i=1, size(prm%A,3) do i=1, size(prm%A,3)
prm%A(1:3,1:3,i) = lattice_symmetrize_33(prm%A(1:3,1:3,i),phase_lattice(p)) prm%A(1:3,1:3,i) = lattice_symmetrize_33(prm%A(1:3,1:3,i),phase_lattice(p))
enddo end do
end associate end associate
endif end if
enddo end do
enddo end do
end function thermalexpansion_init end function thermalexpansion_init

12
src/phase_mechanical_elastic.f90
View File

@ -26,10 +26,10 @@ module subroutine elastic_init(phases)
elastic elastic
print'(/,a)', ' <<<+- phase:mechanical:elastic init -+>>>' print'(/,1x,a)', '<<<+- phase:mechanical:elastic init -+>>>'
print'(/,a)', ' <<<+- phase:mechanical:elastic:Hooke init -+>>>' print'(/,1x,a)', '<<<+- phase:mechanical:elastic:Hooke init -+>>>'
print'(a,i0)', ' # phases: ',phases%length; flush(IO_STDOUT) print'(/,a,i0)', ' # phases: ',phases%length; flush(IO_STDOUT)
allocate(param(phases%length)) allocate(param(phases%length))
@ -48,7 +48,7 @@ module subroutine elastic_init(phases)
if (any(phase_lattice(ph) == ['hP','tI'])) then if (any(phase_lattice(ph) == ['hP','tI'])) then
prm%C66(1,3) = elastic%get_asFloat('C_13') prm%C66(1,3) = elastic%get_asFloat('C_13')
prm%C66(3,3) = elastic%get_asFloat('C_33') prm%C66(3,3) = elastic%get_asFloat('C_33')
endif end if
if (phase_lattice(ph) == 'tI') prm%C66(6,6) = elastic%get_asFloat('C_66') if (phase_lattice(ph) == 'tI') prm%C66(6,6) = elastic%get_asFloat('C_66')
prm%C66 = lattice_symmetrize_C66(prm%C66,phase_lattice(ph)) prm%C66 = lattice_symmetrize_C66(prm%C66,phase_lattice(ph))
@ -59,7 +59,7 @@ module subroutine elastic_init(phases)
prm%C66 = math_sym3333to66(math_Voigt66to3333(prm%C66)) ! Literature data is in Voigt notation prm%C66 = math_sym3333to66(math_Voigt66to3333(prm%C66)) ! Literature data is in Voigt notation
end associate end associate
enddo end do
end subroutine elastic_init end subroutine elastic_init
@ -98,7 +98,7 @@ module subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
do i =1, 3;do j=1,3 do i =1, 3;do j=1,3
dS_dFe(i,j,1:3,1:3) = matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko dS_dFe(i,j,1:3,1:3) = matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko
dS_dFi(i,j,1:3,1:3) = 2.0_pReal*matmul(matmul(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn dS_dFi(i,j,1:3,1:3) = 2.0_pReal*matmul(matmul(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn
enddo; enddo end do; end do
end subroutine phase_hooke_SandItsTangents end subroutine phase_hooke_SandItsTangents

12
src/phase_mechanical_plastic.f90
View File

@ -222,7 +222,7 @@ contains
module subroutine plastic_init module subroutine plastic_init
print'(/,a)', ' <<<+- phase:mechanical:plastic init -+>>>' print'(/,1x,a)', '<<<+- phase:mechanical:plastic init -+>>>'
where(plastic_none_init()) phase_plasticity = PLASTICITY_NONE_ID where(plastic_none_init()) phase_plasticity = PLASTICITY_NONE_ID
where(plastic_isotropic_init()) phase_plasticity = PLASTICITY_ISOTROPIC_ID where(plastic_isotropic_init()) phase_plasticity = PLASTICITY_ISOTROPIC_ID
@ -295,7 +295,7 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + & dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + &
matmul(matmul(Fi,dLp_dMp(i,j,1:3,1:3)),S) matmul(matmul(Fi,dLp_dMp(i,j,1:3,1:3)),S)
dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi) dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi)
enddo; enddo end do; end do
end subroutine plastic_LpAndItsTangents end subroutine plastic_LpAndItsTangents
@ -418,7 +418,7 @@ module function plastic_deltaState(ph, en) result(broken)
end select plasticType end select plasticType
if(.not. broken) then if (.not. broken) then
select case(phase_plasticity(ph)) select case(phase_plasticity(ph))
case (PLASTICITY_NONLOCAL_ID,PLASTICITY_KINEHARDENING_ID) case (PLASTICITY_NONLOCAL_ID,PLASTICITY_KINEHARDENING_ID)
@ -427,7 +427,7 @@ module function plastic_deltaState(ph, en) result(broken)
plasticState(ph)%state(myOffset + 1:myOffset + mySize,en) = & plasticState(ph)%state(myOffset + 1:myOffset + mySize,en) = &
plasticState(ph)%state(myOffset + 1:myOffset + mySize,en) + plasticState(ph)%deltaState(1:mySize,en) plasticState(ph)%state(myOffset + 1:myOffset + mySize,en) + plasticState(ph)%deltaState(1:mySize,en)
end select end select
endif end if
end function plastic_deltaState end function plastic_deltaState
@ -453,8 +453,8 @@ function plastic_active(plastic_label) result(active_plastic)
phase => phases%get(ph) phase => phases%get(ph)
mech => phase%get('mechanical') mech => phase%get('mechanical')
pl => mech%get('plastic',defaultVal = emptyDict) pl => mech%get('plastic',defaultVal = emptyDict)
if(pl%get_asString('type',defaultVal='none') == plastic_label) active_plastic(ph) = .true. if (pl%get_asString('type',defaultVal='none') == plastic_label) active_plastic(ph) = .true.
enddo end do
end function plastic_active end function plastic_active

14
src/phase_mechanical_plastic_dislotungsten.f90
View File

@ -99,13 +99,13 @@ module function plastic_dislotungsten_init() result(myPlasticity)
myPlasticity = plastic_active('dislotungsten') myPlasticity = plastic_active('dislotungsten')
if(count(myPlasticity) == 0) return if (count(myPlasticity) == 0) return
print'(/,a)', ' <<<+- phase:mechanical:plastic:dislotungsten init -+>>>' print'(/,1x,a)', '<<<+- phase:mechanical:plastic:dislotungsten init -+>>>'
print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
print*, 'D. Cereceda et al., International Journal of Plasticity 78:242–256, 2016' print'(/,1x,a)', 'D. Cereceda et al., International Journal of Plasticity 78:242–256, 2016'
print*, 'https://doi.org/10.1016/j.ijplas.2015.09.002' print'( 1x,a)', 'https://doi.org/10.1016/j.ijplas.2015.09.002'
phases => config_material%get('phase') phases => config_material%get('phase')
@ -116,7 +116,7 @@ module function plastic_dislotungsten_init() result(myPlasticity)
do ph = 1, phases%length do ph = 1, phases%length
if(.not. myPlasticity(ph)) cycle if (.not. myPlasticity(ph)) cycle
associate(prm => param(ph), dot => dotState(ph), stt => state(ph), dst => dependentState(ph)) associate(prm => param(ph), dot => dotState(ph), stt => state(ph), dst => dependentState(ph))
@ -243,7 +243,7 @@ module function plastic_dislotungsten_init() result(myPlasticity)
stt%gamma_sl => plasticState(ph)%state(startIndex:endIndex,:) stt%gamma_sl => plasticState(ph)%state(startIndex:endIndex,:)
dot%gamma_sl => plasticState(ph)%dotState(startIndex:endIndex,:) dot%gamma_sl => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
allocate(dst%Lambda_sl(prm%sum_N_sl,Nmembers), source=0.0_pReal) allocate(dst%Lambda_sl(prm%sum_N_sl,Nmembers), source=0.0_pReal)
allocate(dst%tau_pass(prm%sum_N_sl,Nmembers), source=0.0_pReal) allocate(dst%tau_pass(prm%sum_N_sl,Nmembers), source=0.0_pReal)

72
src/phase_mechanical_plastic_dislotwin.f90
View File

@ -150,17 +150,17 @@ module function plastic_dislotwin_init() result(myPlasticity)
myPlasticity = plastic_active('dislotwin') myPlasticity = plastic_active('dislotwin')
if(count(myPlasticity) == 0) return if(count(myPlasticity) == 0) return
print'(/,a)', ' <<<+- phase:mechanical:plastic:dislotwin init -+>>>' print'(/,1x,a)', '<<<+- phase:mechanical:plastic:dislotwin init -+>>>'
print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
print*, 'A. Ma and F. Roters, Acta Materialia 52(12):3603–3612, 2004' print'(/,1x,a)', 'A. Ma and F. Roters, Acta Materialia 52(12):3603–3612, 2004'
print*, 'https://doi.org/10.1016/j.actamat.2004.04.012'//IO_EOL print'( 1x,a)', 'https://doi.org/10.1016/j.actamat.2004.04.012'//IO_EOL
print*, 'F. Roters et al., Computational Materials Science 39:91–95, 2007' print'(/,1x,a)', 'F. Roters et al., Computational Materials Science 39:91–95, 2007'
print*, 'https://doi.org/10.1016/j.commatsci.2006.04.014'//IO_EOL print'( 1x,a)', 'https://doi.org/10.1016/j.commatsci.2006.04.014'//IO_EOL
print*, 'S.L. Wong et al., Acta Materialia 118:140–151, 2016' print'(/,1x,a)', 'S.L. Wong et al., Acta Materialia 118:140–151, 2016'
print*, 'https://doi.org/10.1016/j.actamat.2016.07.032' print'( 1x,a)', 'https://doi.org/10.1016/j.actamat.2016.07.032'
phases => config_material%get('phase') phases => config_material%get('phase')
@ -256,7 +256,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
rho_mob_0 = emptyRealArray; rho_dip_0 = emptyRealArray rho_mob_0 = emptyRealArray; rho_dip_0 = emptyRealArray
allocate(prm%b_sl,prm%Q_sl,prm%v_0,prm%i_sl,prm%p,prm%q,prm%B,source=emptyRealArray) allocate(prm%b_sl,prm%Q_sl,prm%v_0,prm%i_sl,prm%p,prm%q,prm%B,source=emptyRealArray)
allocate(prm%forestProjection(0,0),prm%h_sl_sl(0,0)) allocate(prm%forestProjection(0,0),prm%h_sl_sl(0,0))
endif slipActive end if slipActive
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! twin related parameters ! twin related parameters
@ -283,7 +283,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
if (.not. prm%fccTwinTransNucleation) then if (.not. prm%fccTwinTransNucleation) then
prm%dot_N_0_tw = pl%get_as1dFloat('dot_N_0_tw') prm%dot_N_0_tw = pl%get_as1dFloat('dot_N_0_tw')
prm%dot_N_0_tw = math_expand(prm%dot_N_0_tw,N_tw) prm%dot_N_0_tw = math_expand(prm%dot_N_0_tw,N_tw)
endif end if
! expand: family => system ! expand: family => system
prm%b_tw = math_expand(prm%b_tw,N_tw) prm%b_tw = math_expand(prm%b_tw,N_tw)
@ -299,11 +299,11 @@ module function plastic_dislotwin_init() result(myPlasticity)
if (any(prm%r < 0.0_pReal)) extmsg = trim(extmsg)//' p_tw' if (any(prm%r < 0.0_pReal)) extmsg = trim(extmsg)//' p_tw'
if (.not. prm%fccTwinTransNucleation) then if (.not. prm%fccTwinTransNucleation) then
if (any(prm%dot_N_0_tw < 0.0_pReal)) extmsg = trim(extmsg)//' dot_N_0_tw' if (any(prm%dot_N_0_tw < 0.0_pReal)) extmsg = trim(extmsg)//' dot_N_0_tw'
endif end if
else twinActive else twinActive
allocate(prm%gamma_char,prm%b_tw,prm%dot_N_0_tw,prm%t_tw,prm%r,source=emptyRealArray) allocate(prm%gamma_char,prm%b_tw,prm%dot_N_0_tw,prm%t_tw,prm%r,source=emptyRealArray)
allocate(prm%h_tw_tw(0,0)) allocate(prm%h_tw_tw(0,0))
endif twinActive end if twinActive
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! transformation related parameters ! transformation related parameters
@ -335,7 +335,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
if (phase_lattice(ph) /= 'cF') then if (phase_lattice(ph) /= 'cF') then
prm%dot_N_0_tr = pl%get_as1dFloat('dot_N_0_tr') prm%dot_N_0_tr = pl%get_as1dFloat('dot_N_0_tr')
prm%dot_N_0_tr = math_expand(prm%dot_N_0_tr,N_tr) prm%dot_N_0_tr = math_expand(prm%dot_N_0_tr,N_tr)
endif end if
prm%t_tr = pl%get_as1dFloat('t_tr') prm%t_tr = pl%get_as1dFloat('t_tr')
prm%t_tr = math_expand(prm%t_tr,N_tr) prm%t_tr = math_expand(prm%t_tr,N_tr)
prm%s = pl%get_as1dFloat('p_tr',defaultVal=[0.0_pReal]) prm%s = pl%get_as1dFloat('p_tr',defaultVal=[0.0_pReal])
@ -349,11 +349,11 @@ module function plastic_dislotwin_init() result(myPlasticity)
if (any(prm%s < 0.0_pReal)) extmsg = trim(extmsg)//' p_tr' if (any(prm%s < 0.0_pReal)) extmsg = trim(extmsg)//' p_tr'
if (phase_lattice(ph) /= 'cF') then if (phase_lattice(ph) /= 'cF') then
if (any(prm%dot_N_0_tr < 0.0_pReal)) extmsg = trim(extmsg)//' dot_N_0_tr' if (any(prm%dot_N_0_tr < 0.0_pReal)) extmsg = trim(extmsg)//' dot_N_0_tr'
endif end if
else transActive else transActive
allocate(prm%s,prm%b_tr,prm%t_tr,prm%dot_N_0_tr,source=emptyRealArray) allocate(prm%s,prm%b_tr,prm%t_tr,prm%dot_N_0_tr,source=emptyRealArray)
allocate(prm%h_tr_tr(0,0)) allocate(prm%h_tr_tr(0,0))
endif transActive end if transActive
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! shearband related parameters ! shearband related parameters
@ -369,7 +369,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
if (prm%E_sb < 0.0_pReal) extmsg = trim(extmsg)//' Q_sb' if (prm%E_sb < 0.0_pReal) extmsg = trim(extmsg)//' Q_sb'
if (prm%p_sb <= 0.0_pReal) extmsg = trim(extmsg)//' p_sb' if (prm%p_sb <= 0.0_pReal) extmsg = trim(extmsg)//' p_sb'
if (prm%q_sb <= 0.0_pReal) extmsg = trim(extmsg)//' q_sb' if (prm%q_sb <= 0.0_pReal) extmsg = trim(extmsg)//' q_sb'
endif end if
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! parameters required for several mechanisms and their interactions ! parameters required for several mechanisms and their interactions
@ -383,19 +383,19 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%T_ref = pl%get_asFloat('T_ref') prm%T_ref = pl%get_asFloat('T_ref')
prm%Gamma_sf(1) = pl%get_asFloat('Gamma_sf') prm%Gamma_sf(1) = pl%get_asFloat('Gamma_sf')
prm%Gamma_sf(2) = pl%get_asFloat('Gamma_sf,T',defaultVal=0.0_pReal) prm%Gamma_sf(2) = pl%get_asFloat('Gamma_sf,T',defaultVal=0.0_pReal)
endif end if
slipAndTwinActive: if (prm%sum_N_sl * prm%sum_N_tw > 0) then slipAndTwinActive: if (prm%sum_N_sl * prm%sum_N_tw > 0) then
prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,N_tw,pl%get_as1dFloat('h_sl-tw'), & prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,N_tw,pl%get_as1dFloat('h_sl-tw'), &
phase_lattice(ph)) phase_lattice(ph))
if (prm%fccTwinTransNucleation .and. size(N_tw) /= 1) extmsg = trim(extmsg)//' N_tw: nucleation' if (prm%fccTwinTransNucleation .and. size(N_tw) /= 1) extmsg = trim(extmsg)//' N_tw: nucleation'
endif slipAndTwinActive end if slipAndTwinActive
slipAndTransActive: if (prm%sum_N_sl * prm%sum_N_tr > 0) then slipAndTransActive: if (prm%sum_N_sl * prm%sum_N_tr > 0) then
prm%h_sl_tr = lattice_interaction_SlipByTrans(N_sl,N_tr,pl%get_as1dFloat('h_sl-tr'), & prm%h_sl_tr = lattice_interaction_SlipByTrans(N_sl,N_tr,pl%get_as1dFloat('h_sl-tr'), &
phase_lattice(ph)) phase_lattice(ph))
if (prm%fccTwinTransNucleation .and. size(N_tr) /= 1) extmsg = trim(extmsg)//' N_tr: nucleation' if (prm%fccTwinTransNucleation .and. size(N_tr) /= 1) extmsg = trim(extmsg)//' N_tr: nucleation'
endif slipAndTransActive end if slipAndTransActive
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
@ -465,7 +465,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
! exit if any parameter is out of range ! exit if any parameter is out of range
if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(dislotwin)') if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(dislotwin)')
enddo end do
end function plastic_dislotwin_init end function plastic_dislotwin_init
@ -494,11 +494,11 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC)
do i=1,prm%sum_N_tw do i=1,prm%sum_N_tw
homogenizedC = homogenizedC & homogenizedC = homogenizedC &
+ stt%f_tw(i,en)*prm%C66_tw(1:6,1:6,i) + stt%f_tw(i,en)*prm%C66_tw(1:6,1:6,i)
enddo end do
do i=1,prm%sum_N_tr do i=1,prm%sum_N_tr
homogenizedC = homogenizedC & homogenizedC = homogenizedC &
+ stt%f_tr(i,en)*prm%C66_tr(1:6,1:6,i) + stt%f_tr(i,en)*prm%C66_tr(1:6,1:6,i)
enddo end do
end associate end associate
@ -566,7 +566,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ ddot_gamma_dtau_sl(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i) + ddot_gamma_dtau_sl(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i)
enddo slipContribution end do slipContribution
call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw) call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw)
twinContibution: do i = 1, prm%sum_N_tw twinContibution: do i = 1, prm%sum_N_tw
@ -574,7 +574,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ ddot_gamma_dtau_tw(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) + ddot_gamma_dtau_tw(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i)
enddo twinContibution end do twinContibution
call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr,ddot_gamma_dtau_tr) call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr,ddot_gamma_dtau_tr)
transContibution: do i = 1, prm%sum_N_tr transContibution: do i = 1, prm%sum_N_tr
@ -582,7 +582,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ ddot_gamma_dtau_tr(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i) + ddot_gamma_dtau_tr(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i)
enddo transContibution end do transContibution
Lp = Lp * f_unrotated Lp = Lp * f_unrotated
dLp_dMp = dLp_dMp * f_unrotated dLp_dMp = dLp_dMp * f_unrotated
@ -608,10 +608,10 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ ddot_gamma_dtau * P_sb(k,l) * P_sb(m,n) + ddot_gamma_dtau * P_sb(k,l) * P_sb(m,n)
endif significantShearBandStress end if significantShearBandStress
enddo end do
endif shearBandingContribution end if shearBandingContribution
end associate end associate
@ -686,9 +686,9 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,en) & dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,en) &
/ (d_hat-prm%d_caron(i)) / (d_hat-prm%d_caron(i))
endif end if
endif significantSlipStress end if significantSlipStress
enddo slipState end do slipState
dot%rho_mob(:,en) = abs(dot_gamma_sl)/(prm%b_sl*dst%Lambda_sl(:,en)) & dot%rho_mob(:,en) = abs(dot_gamma_sl)/(prm%b_sl*dst%Lambda_sl(:,en)) &
- dot_rho_dip_formation & - dot_rho_dip_formation &
@ -833,7 +833,7 @@ module subroutine plastic_dislotwin_results(ph,group)
end select end select
enddo end do
end associate end associate
@ -968,8 +968,8 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
end if end if
else isFCC else isFCC
Ndot0=prm%dot_N_0_tw(i) Ndot0=prm%dot_N_0_tw(i)
endif isFCC end if isFCC
enddo end do
significantStress: where(tau > tol_math_check) significantStress: where(tau > tol_math_check)
StressRatio_r = (dst%tau_hat_tw(:,en)/tau)**prm%r StressRatio_r = (dst%tau_hat_tw(:,en)/tau)**prm%r
@ -1037,8 +1037,8 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
end if end if
else isFCC else isFCC
Ndot0=prm%dot_N_0_tr(i) Ndot0=prm%dot_N_0_tr(i)
endif isFCC end if isFCC
enddo end do
significantStress: where(tau > tol_math_check) significantStress: where(tau > tol_math_check)
StressRatio_s = (dst%tau_hat_tr(:,en)/tau)**prm%s StressRatio_s = (dst%tau_hat_tr(:,en)/tau)**prm%s

20
src/phase_mechanical_plastic_isotropic.f90
View File

@ -68,11 +68,11 @@ module function plastic_isotropic_init() result(myPlasticity)
myPlasticity = plastic_active('isotropic') myPlasticity = plastic_active('isotropic')
if(count(myPlasticity) == 0) return if(count(myPlasticity) == 0) return
print'(/,a)', ' <<<+- phase:mechanical:plastic:isotropic init -+>>>' print'(/,1x,a)', '<<<+- phase:mechanical:plastic:isotropic init -+>>>'
print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
print*, 'T. Maiti and P. Eisenlohr, Scripta Materialia 145:37–40, 2018' print'(/,a)', 'T. Maiti and P. Eisenlohr, Scripta Materialia 145:37–40, 2018'
print*, 'https://doi.org/10.1016/j.scriptamat.2017.09.047' print'(/,a)', 'https://doi.org/10.1016/j.scriptamat.2017.09.047'
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(param(phases%length)) allocate(param(phases%length))
@ -140,7 +140,7 @@ module function plastic_isotropic_init() result(myPlasticity)
! exit if any parameter is out of range ! exit if any parameter is out of range
if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(isotropic)') if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(isotropic)')
enddo end do
end function plastic_isotropic_init end function plastic_isotropic_init
@ -232,7 +232,7 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,en)
else else
Li = 0.0_pReal Li = 0.0_pReal
dLi_dMi = 0.0_pReal dLi_dMi = 0.0_pReal
endif end if
end associate end associate
@ -262,7 +262,7 @@ module subroutine isotropic_dotState(Mp,ph,en)
norm_Mp = sqrt(math_tensordot(Mp,Mp)) norm_Mp = sqrt(math_tensordot(Mp,Mp))
else else
norm_Mp = sqrt(math_tensordot(math_deviatoric33(Mp),math_deviatoric33(Mp))) norm_Mp = sqrt(math_tensordot(math_deviatoric33(Mp),math_deviatoric33(Mp)))
endif end if
dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp /(prm%M*stt%xi(en))) **prm%n dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp /(prm%M*stt%xi(en))) **prm%n
@ -273,13 +273,13 @@ module subroutine isotropic_dotState(Mp,ph,en)
xi_inf_star = prm%xi_inf & xi_inf_star = prm%xi_inf &
+ asinh( (dot_gamma / prm%c_1)**(1.0_pReal / prm%c_2))**(1.0_pReal / prm%c_3) & + asinh( (dot_gamma / prm%c_1)**(1.0_pReal / prm%c_2))**(1.0_pReal / prm%c_3) &
/ prm%c_4 * (dot_gamma / prm%dot_gamma_0)**(1.0_pReal / prm%n) / prm%c_4 * (dot_gamma / prm%dot_gamma_0)**(1.0_pReal / prm%n)
endif end if
dot%xi(en) = dot_gamma & dot%xi(en) = dot_gamma &
* ( prm%h_0 + prm%h_ln * log(dot_gamma) ) & * ( prm%h_0 + prm%h_ln * log(dot_gamma) ) &
* sign(abs(1.0_pReal - stt%xi(en)/xi_inf_star)**prm%a *prm%h, 1.0_pReal-stt%xi(en)/xi_inf_star) * sign(abs(1.0_pReal - stt%xi(en)/xi_inf_star)**prm%a *prm%h, 1.0_pReal-stt%xi(en)/xi_inf_star)
else else
dot%xi(en) = 0.0_pReal dot%xi(en) = 0.0_pReal
endif end if
end associate end associate
@ -303,7 +303,7 @@ module subroutine plastic_isotropic_results(ph,group)
call results_writeDataset(stt%xi,group,trim(prm%output(o)), & call results_writeDataset(stt%xi,group,trim(prm%output(o)), &
'resistance against plastic flow','Pa') 'resistance against plastic flow','Pa')
end select end select
enddo outputsLoop end do outputsLoop
end associate end associate
end subroutine plastic_isotropic_results end subroutine plastic_isotropic_results

22
src/phase_mechanical_plastic_kinehardening.f90
View File

@ -83,8 +83,8 @@ module function plastic_kinehardening_init() result(myPlasticity)
myPlasticity = plastic_active('kinehardening') myPlasticity = plastic_active('kinehardening')
if(count(myPlasticity) == 0) return if(count(myPlasticity) == 0) return
print'(/,a)', ' <<<+- phase:mechanical:plastic:kinehardening init -+>>>' print'(/,1x,a)', '<<<+- phase:mechanical:plastic:kinehardening init -+>>>'
print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
phases => config_material%get('phase') phases => config_material%get('phase')
@ -95,7 +95,7 @@ module function plastic_kinehardening_init() result(myPlasticity)
do ph = 1, phases%length do ph = 1, phases%length
if(.not. myPlasticity(ph)) cycle if (.not. myPlasticity(ph)) cycle
associate(prm => param(ph), dot => dotState(ph), dlt => deltaState(ph), stt => state(ph)) associate(prm => param(ph), dot => dotState(ph), dlt => deltaState(ph), stt => state(ph))
@ -125,7 +125,7 @@ module function plastic_kinehardening_init() result(myPlasticity)
else else
prm%P_nS_pos = prm%P prm%P_nS_pos = prm%P
prm%P_nS_neg = prm%P prm%P_nS_neg = prm%P
endif end if
prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dFloat('h_sl-sl'), & prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dFloat('h_sl-sl'), &
phase_lattice(ph)) phase_lattice(ph))
@ -161,7 +161,7 @@ module function plastic_kinehardening_init() result(myPlasticity)
xi_0 = emptyRealArray xi_0 = emptyRealArray
allocate(prm%xi_inf_f,prm%xi_inf_b,prm%h_0_f,prm%h_inf_f,prm%h_0_b,prm%h_inf_b,source=emptyRealArray) allocate(prm%xi_inf_f,prm%xi_inf_b,prm%h_0_f,prm%h_inf_f,prm%h_0_b,prm%h_inf_b,source=emptyRealArray)
allocate(prm%h_sl_sl(0,0)) allocate(prm%h_sl_sl(0,0))
endif slipActive end if slipActive
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
@ -217,7 +217,7 @@ module function plastic_kinehardening_init() result(myPlasticity)
! exit if any parameter is out of range ! exit if any parameter is out of range
if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(kinehardening)') if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(kinehardening)')
enddo end do
end function plastic_kinehardening_init end function plastic_kinehardening_init
@ -258,7 +258,7 @@ pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ ddot_gamma_dtau_pos(i) * prm%P(k,l,i) * prm%P_nS_pos(m,n,i) & + ddot_gamma_dtau_pos(i) * prm%P(k,l,i) * prm%P_nS_pos(m,n,i) &
+ ddot_gamma_dtau_neg(i) * prm%P(k,l,i) * prm%P_nS_neg(m,n,i) + ddot_gamma_dtau_neg(i) * prm%P(k,l,i) * prm%P_nS_neg(m,n,i)
enddo end do
end associate end associate
@ -382,7 +382,7 @@ module subroutine plastic_kinehardening_results(ph,group)
'plastic shear','1',prm%systems_sl) 'plastic shear','1',prm%systems_sl)
end select end select
enddo end do
end associate end associate
@ -424,7 +424,7 @@ pure subroutine kinetics(Mp,ph,en, &
tau_pos(i) = math_tensordot(Mp,prm%P_nS_pos(1:3,1:3,i)) - stt%chi(i,en) tau_pos(i) = math_tensordot(Mp,prm%P_nS_pos(1:3,1:3,i)) - stt%chi(i,en)
tau_neg(i) = merge(math_tensordot(Mp,prm%P_nS_neg(1:3,1:3,i)) - stt%chi(i,en), & tau_neg(i) = merge(math_tensordot(Mp,prm%P_nS_neg(1:3,1:3,i)) - stt%chi(i,en), &
0.0_pReal, prm%nonSchmidActive) 0.0_pReal, prm%nonSchmidActive)
enddo end do
where(dNeq0(tau_pos)) where(dNeq0(tau_pos))
dot_gamma_pos = prm%dot_gamma_0 * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active dot_gamma_pos = prm%dot_gamma_0 * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active
@ -446,14 +446,14 @@ pure subroutine kinetics(Mp,ph,en, &
else where else where
ddot_gamma_dtau_pos = 0.0_pReal ddot_gamma_dtau_pos = 0.0_pReal
end where end where
endif end if
if (present(ddot_gamma_dtau_neg)) then if (present(ddot_gamma_dtau_neg)) then
where(dNeq0(dot_gamma_neg)) where(dNeq0(dot_gamma_neg))
ddot_gamma_dtau_neg = dot_gamma_neg*prm%n/tau_neg ddot_gamma_dtau_neg = dot_gamma_neg*prm%n/tau_neg
else where else where
ddot_gamma_dtau_neg = 0.0_pReal ddot_gamma_dtau_neg = 0.0_pReal
end where end where
endif end if
end associate end associate

10
src/phase_mechanical_plastic_none.f90
View File

@ -22,16 +22,16 @@ module function plastic_none_init() result(myPlasticity)
myPlasticity = plastic_active('none') myPlasticity = plastic_active('none')
if(count(myPlasticity) == 0) return if (count(myPlasticity) == 0) return
print'(/,a)', ' <<<+- phase:mechanical:plastic:none init -+>>>' print'(/,1x,a)', '<<<+- phase:mechanical:plastic:none init -+>>>'
print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
phases => config_material%get('phase') phases => config_material%get('phase')
do ph = 1, phases%length do ph = 1, phases%length
if(.not. myPlasticity(ph)) cycle if (.not. myPlasticity(ph)) cycle
call phase_allocateState(plasticState(ph),count(material_phaseID == ph),0,0,0) call phase_allocateState(plasticState(ph),count(material_phaseID == ph),0,0,0)
enddo end do
end function plastic_none_init end function plastic_none_init

140
src/phase_mechanical_plastic_nonlocal.f90
View File

@ -197,19 +197,19 @@ module function plastic_nonlocal_init() result(myPlasticity)
myPlasticity = plastic_active('nonlocal') myPlasticity = plastic_active('nonlocal')
Ninstances = count(myPlasticity) Ninstances = count(myPlasticity)
if(Ninstances == 0) then if (Ninstances == 0) then
call geometry_plastic_nonlocal_disable call geometry_plastic_nonlocal_disable
return return
endif end if
print'(/,a)', ' <<<+- phase:mechanical:plastic:nonlocal init -+>>>' print'(/,1x,a)', '<<<+- phase:mechanical:plastic:nonlocal init -+>>>'
print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT) print'(/,a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT)
print*, 'C. Reuber et al., Acta Materialia 71:333–348, 2014' print'(/,1x,a)', 'C. Reuber et al., Acta Materialia 71:333–348, 2014'
print*, 'https://doi.org/10.1016/j.actamat.2014.03.012'//IO_EOL print'( 1x,a)', 'https://doi.org/10.1016/j.actamat.2014.03.012'//IO_EOL
print*, 'C. Kords, Dissertation RWTH Aachen, 2014' print'(/,1x,a)', 'C. Kords, Dissertation RWTH Aachen, 2014'
print*, 'http://publications.rwth-aachen.de/record/229993' print'( 1x,a)', 'http://publications.rwth-aachen.de/record/229993'
phases => config_material%get('phase') phases => config_material%get('phase')
@ -224,7 +224,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
allocate(dependentState(phases%length)) allocate(dependentState(phases%length))
do ph = 1, phases%length do ph = 1, phases%length
if(.not. myPlasticity(ph)) cycle if (.not. myPlasticity(ph)) cycle
associate(prm => param(ph), dot => dotState(ph), stt => state(ph), & associate(prm => param(ph), dot => dotState(ph), stt => state(ph), &
st0 => state0(ph), del => deltaState(ph), dst => dependentState(ph)) st0 => state0(ph), del => deltaState(ph), dst => dependentState(ph))
@ -259,7 +259,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
else else
prm%P_nS_pos = prm%P_sl prm%P_nS_pos = prm%P_sl
prm%P_nS_neg = prm%P_sl prm%P_nS_neg = prm%P_sl
endif end if
prm%h_sl_sl = lattice_interaction_SlipBySlip(ini%N_sl,pl%get_as1dFloat('h_sl-sl'), & prm%h_sl_sl = lattice_interaction_SlipBySlip(ini%N_sl,pl%get_as1dFloat('h_sl-sl'), &
phase_lattice(ph)) phase_lattice(ph))
@ -280,8 +280,8 @@ module function plastic_nonlocal_init() result(myPlasticity)
if (all(dEq0 (math_cross(prm%slip_direction(1:3,s1),prm%slip_direction(1:3,s2)))) .and. & if (all(dEq0 (math_cross(prm%slip_direction(1:3,s1),prm%slip_direction(1:3,s2)))) .and. &
any(dNeq0(math_cross(prm%slip_normal (1:3,s1),prm%slip_normal (1:3,s2))))) & any(dNeq0(math_cross(prm%slip_normal (1:3,s1),prm%slip_normal (1:3,s2))))) &
prm%colinearSystem(s1) = s2 prm%colinearSystem(s1) = s2
enddo end do
enddo end do
ini%rho_u_ed_pos_0 = pl%get_as1dFloat('rho_u_ed_pos_0', requiredSize=size(ini%N_sl)) ini%rho_u_ed_pos_0 = pl%get_as1dFloat('rho_u_ed_pos_0', requiredSize=size(ini%N_sl))
ini%rho_u_ed_neg_0 = pl%get_as1dFloat('rho_u_ed_neg_0', requiredSize=size(ini%N_sl)) ini%rho_u_ed_neg_0 = pl%get_as1dFloat('rho_u_ed_neg_0', requiredSize=size(ini%N_sl))
@ -391,7 +391,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
if (prm%f_ed_mult < 0.0_pReal .or. prm%f_ed_mult > 1.0_pReal) & if (prm%f_ed_mult < 0.0_pReal .or. prm%f_ed_mult > 1.0_pReal) &
extmsg = trim(extmsg)//' f_ed_mult' extmsg = trim(extmsg)//' f_ed_mult'
endif slipActive end if slipActive
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
@ -506,7 +506,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
! exit if any parameter is out of range ! exit if any parameter is out of range
if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(nonlocal)') if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(nonlocal)')
enddo end do
allocate(compatibility(2,maxval(param%sum_N_sl),maxval(param%sum_N_sl),nIPneighbors,& allocate(compatibility(2,maxval(param%sum_N_sl),maxval(param%sum_N_sl),nIPneighbors,&
discretization_nIPs,discretization_Nelems), source=0.0_pReal) discretization_nIPs,discretization_Nelems), source=0.0_pReal)
@ -527,24 +527,24 @@ module function plastic_nonlocal_init() result(myPlasticity)
do s = 1,param(ph)%sum_N_sl do s = 1,param(ph)%sum_N_sl
l = l + 1 l = l + 1
iRhoU(s,t,ph) = l iRhoU(s,t,ph) = l
enddo end do
enddo end do
l = l + (4+2+1+1)*param(ph)%sum_N_sl ! immobile(4), dipole(2), shear, forest l = l + (4+2+1+1)*param(ph)%sum_N_sl ! immobile(4), dipole(2), shear, forest
do t = 1,4 do t = 1,4
do s = 1,param(ph)%sum_N_sl do s = 1,param(ph)%sum_N_sl
l = l + 1 l = l + 1
iV(s,t,ph) = l iV(s,t,ph) = l
enddo end do
enddo end do
do t = 1,2 do t = 1,2
do s = 1,param(ph)%sum_N_sl do s = 1,param(ph)%sum_N_sl
l = l + 1 l = l + 1
iD(s,t,ph) = l iD(s,t,ph) = l
enddo end do
enddo end do
if (iD(param(ph)%sum_N_sl,2,ph) /= plasticState(ph)%sizeState) & if (iD(param(ph)%sum_N_sl,2,ph) /= plasticState(ph)%sizeState) &
error stop 'state indices not properly set (nonlocal)' error stop 'state indices not properly set (nonlocal)'
enddo end do
end function plastic_nonlocal_init end function plastic_nonlocal_init
@ -625,7 +625,7 @@ module subroutine nonlocal_dependentState(ph, en, ip, el)
/ log(0.35_pReal * prm%b_sl * 1e6_pReal))** 2.0_pReal,2,prm%sum_N_sl) / log(0.35_pReal * prm%b_sl * 1e6_pReal))** 2.0_pReal,2,prm%sum_N_sl)
else else
myInteractionMatrix = prm%h_sl_sl myInteractionMatrix = prm%h_sl_sl
endif end if
dst%tau_pass(:,en) = prm%mu * prm%b_sl & dst%tau_pass(:,en) = prm%mu * prm%b_sl &
* sqrt(matmul(myInteractionMatrix,sum(abs(rho),2))) * sqrt(matmul(myInteractionMatrix,sum(abs(rho),2)))
@ -680,14 +680,14 @@ module subroutine nonlocal_dependentState(ph, en, ip, el)
connection_latticeConf(1:3,n) = 0.0_pReal connection_latticeConf(1:3,n) = 0.0_pReal
rho_edg_delta_neighbor(:,n) = rho_edg_delta rho_edg_delta_neighbor(:,n) = rho_edg_delta
rho_scr_delta_neighbor(:,n) = rho_scr_delta rho_scr_delta_neighbor(:,n) = rho_scr_delta
endif end if
else else
! free surface -> use central values instead ! free surface -> use central values instead
connection_latticeConf(1:3,n) = 0.0_pReal connection_latticeConf(1:3,n) = 0.0_pReal
rho_edg_delta_neighbor(:,n) = rho_edg_delta rho_edg_delta_neighbor(:,n) = rho_edg_delta
rho_scr_delta_neighbor(:,n) = rho_scr_delta rho_scr_delta_neighbor(:,n) = rho_scr_delta
endif end if
enddo end do
neighbor_rhoExcess(1,:,:) = rho_edg_delta_neighbor neighbor_rhoExcess(1,:,:) = rho_edg_delta_neighbor
neighbor_rhoExcess(2,:,:) = rho_scr_delta_neighbor neighbor_rhoExcess(2,:,:) = rho_scr_delta_neighbor
@ -709,12 +709,12 @@ module subroutine nonlocal_dependentState(ph, en, ip, el)
- connection_latticeConf(1:3,neighbors(2)) - connection_latticeConf(1:3,neighbors(2))
rhoExcessDifferences(dir) = neighbor_rhoExcess(c,s,neighbors(1)) & rhoExcessDifferences(dir) = neighbor_rhoExcess(c,s,neighbors(1)) &
- neighbor_rhoExcess(c,s,neighbors(2)) - neighbor_rhoExcess(c,s,neighbors(2))
enddo end do
invConnections = math_inv33(connections) invConnections = math_inv33(connections)
if (all(dEq0(invConnections))) call IO_error(-1,ext_msg='back stress calculation: inversion error') if (all(dEq0(invConnections))) call IO_error(-1,ext_msg='back stress calculation: inversion error')
rhoExcessGradient(c) = math_inner(m(1:3,s,c), matmul(invConnections,rhoExcessDifferences)) rhoExcessGradient(c) = math_inner(m(1:3,s,c), matmul(invConnections,rhoExcessDifferences))
enddo end do
! ... plus gradient from deads ... ! ... plus gradient from deads ...
rhoExcessGradient(1) = rhoExcessGradient(1) + sum(rho(s,imm_edg)) / FVsize rhoExcessGradient(1) = rhoExcessGradient(1) + sum(rho(s,imm_edg)) / FVsize
@ -731,8 +731,8 @@ module subroutine nonlocal_dependentState(ph, en, ip, el)
dst%tau_back(s,en) = - prm%mu * prm%b_sl(s) / (2.0_pReal * PI) & dst%tau_back(s,en) = - prm%mu * prm%b_sl(s) / (2.0_pReal * PI) &
* ( rhoExcessGradient_over_rho(1) / (1.0_pReal - prm%nu) & * ( rhoExcessGradient_over_rho(1) / (1.0_pReal - prm%nu) &
+ rhoExcessGradient_over_rho(2)) + rhoExcessGradient_over_rho(2))
enddo end do
endif end if
end associate end associate
@ -790,8 +790,8 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
else else
tauNS(s,3) = math_tensordot(Mp, +prm%P_nS_neg(1:3,1:3,s)) tauNS(s,3) = math_tensordot(Mp, +prm%P_nS_neg(1:3,1:3,s))
tauNS(s,4) = math_tensordot(Mp, -prm%P_nS_pos(1:3,1:3,s)) tauNS(s,4) = math_tensordot(Mp, -prm%P_nS_pos(1:3,1:3,s))
endif end if
enddo end do
tauNS = tauNS + spread(dst%tau_back(:,en),2,4) tauNS = tauNS + spread(dst%tau_back(:,en),2,4)
tau = tau + dst%tau_back(:,en) tau = tau + dst%tau_back(:,en)
@ -807,12 +807,12 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
do t = 3,4 do t = 3,4
call kinetics(v(:,t), dv_dtau(:,t), dv_dtauNS(:,t), & call kinetics(v(:,t), dv_dtau(:,t), dv_dtauNS(:,t), &
tau, tauNS(:,t), dst%tau_pass(:,en),2,Temperature, ph) tau, tauNS(:,t), dst%tau_pass(:,en),2,Temperature, ph)
enddo end do
else else
v(:,3:4) = spread(v(:,1),2,2) v(:,3:4) = spread(v(:,1),2,2)
dv_dtau(:,3:4) = spread(dv_dtau(:,1),2,2) dv_dtau(:,3:4) = spread(dv_dtau(:,1),2,2)
dv_dtauNS(:,3:4) = spread(dv_dtauNS(:,1),2,2) dv_dtauNS(:,3:4) = spread(dv_dtauNS(:,1),2,2)
endif end if
stt%v(:,en) = pack(v,.true.) stt%v(:,en) = pack(v,.true.)
@ -833,7 +833,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
+ prm%P_sl(i,j,s) & + prm%P_sl(i,j,s) &
* (+ prm%P_nS_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) & * (+ prm%P_nS_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) &
- prm%P_nS_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4)) * prm%b_sl(s) - prm%P_nS_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4)) * prm%b_sl(s)
enddo end do
end associate end associate
@ -899,7 +899,7 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
do s = 1,prm%sum_N_sl do s = 1,prm%sum_N_sl
tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s)) +dst%tau_back(s,en) tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s)) +dst%tau_back(s,en)
if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal
enddo end do
dUpper(:,1) = prm%mu * prm%b_sl/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau)) dUpper(:,1) = prm%mu * prm%b_sl/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau))
dUpper(:,2) = prm%mu * prm%b_sl/(4.0_pReal * PI * abs(tau)) dUpper(:,2) = prm%mu * prm%b_sl/(4.0_pReal * PI * abs(tau))
@ -980,7 +980,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
if (timestep <= 0.0_pReal) then if (timestep <= 0.0_pReal) then
plasticState(ph)%dotState = 0.0_pReal plasticState(ph)%dotState = 0.0_pReal
return return
endif end if
associate(prm => param(ph), dst => dependentState(ph), dot => dotState(ph), stt => state(ph)) associate(prm => param(ph), dst => dependentState(ph), dot => dotState(ph), stt => state(ph))
@ -1002,7 +1002,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
do s = 1,prm%sum_N_sl do s = 1,prm%sum_N_sl
tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s)) + dst%tau_back(s,en) tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s)) + dst%tau_back(s,en)
if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal
enddo end do
dLower = prm%minDipoleHeight dLower = prm%minDipoleHeight
dUpper(:,1) = prm%mu * prm%b_sl/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau)) dUpper(:,1) = prm%mu * prm%b_sl/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau))
@ -1032,7 +1032,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
rhoDotMultiplication(:,1:4) = spread( & rhoDotMultiplication(:,1:4) = spread( &
(sum(abs(dot_gamma(:,1:2)),2) * prm%f_ed_mult + sum(abs(dot_gamma(:,3:4)),2)) & (sum(abs(dot_gamma(:,1:2)),2) * prm%f_ed_mult + sum(abs(dot_gamma(:,3:4)),2)) &
* sqrt(stt%rho_forest(:,en)) / prm%i_sl / prm%b_sl, 2, 4) ! eq. 3.26 * sqrt(stt%rho_forest(:,en)) / prm%i_sl / prm%b_sl, 2, 4) ! eq. 3.26
endif isBCC end if isBCC
forall (s = 1:prm%sum_N_sl, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en) forall (s = 1:prm%sum_N_sl, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en)
@ -1062,7 +1062,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
+ abs(rhoDotSingle2DipoleGlide(:,2*c+4)) & + abs(rhoDotSingle2DipoleGlide(:,2*c+4)) &
- rhoDotSingle2DipoleGlide(:,2*c-1) & - rhoDotSingle2DipoleGlide(:,2*c-1) &
- rhoDotSingle2DipoleGlide(:,2*c) - rhoDotSingle2DipoleGlide(:,2*c)
enddo end do
! athermal annihilation ! athermal annihilation
@ -1103,13 +1103,13 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
if (debugConstitutive%extensive) then if (debugConstitutive%extensive) then
print'(a,i5,a,i2)', '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip print'(a,i5,a,i2)', '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip
print'(a)', '<< CONST >> enforcing cutback !!!' print'(a)', '<< CONST >> enforcing cutback !!!'
endif end if
#endif #endif
plasticState(ph)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN) plasticState(ph)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN)
else else
dot%rho(:,en) = pack(rhoDot,.true.) dot%rho(:,en) = pack(rhoDot,.true.)
dot%gamma(:,en) = sum(dot_gamma,2) dot%gamma(:,en) = sum(dot_gamma,2)
endif end if
end associate end associate
@ -1217,11 +1217,11 @@ function rhoDotFlux(timestep,ph,en,ip,el)
> IPvolume(ip,el) / maxval(IParea(:,ip,el))), & > IPvolume(ip,el) / maxval(IParea(:,ip,el))), &
' at a timestep of ',timestep ' at a timestep of ',timestep
print*, '<< CONST >> enforcing cutback !!!' print*, '<< CONST >> enforcing cutback !!!'
endif end if
#endif #endif
rhoDotFlux = IEEE_value(1.0_pReal,IEEE_quiet_NaN) ! enforce cutback rhoDotFlux = IEEE_value(1.0_pReal,IEEE_quiet_NaN) ! enforce cutback
return return
endif end if
!*** be aware of the definition of slip_transverse = slip_direction x slip_normal !!! !*** be aware of the definition of slip_transverse = slip_direction x slip_normal !!!
@ -1255,7 +1255,7 @@ function rhoDotFlux(timestep,ph,en,ip,el)
Favg = 0.5_pReal * (my_F + neighbor_F) Favg = 0.5_pReal * (my_F + neighbor_F)
else ! if no neighbor, take my value as average else ! if no neighbor, take my value as average
Favg = my_F Favg = my_F
endif end if
neighbor_v0 = 0.0_pReal ! needed for check of sign change in flux density below neighbor_v0 = 0.0_pReal ! needed for check of sign change in flux density below
@ -1300,10 +1300,10 @@ function rhoDotFlux(timestep,ph,en,ip,el)
rhoDotFlux(:,topp) = rhoDotFlux(:,topp) & rhoDotFlux(:,topp) = rhoDotFlux(:,topp) &
+ lineLength/IPvolume(ip,el)*compatibility(c,:,s,n,ip,el)**2.0_pReal ! transferring to opposite signed mobile dislocation type + lineLength/IPvolume(ip,el)*compatibility(c,:,s,n,ip,el)**2.0_pReal ! transferring to opposite signed mobile dislocation type
endif end if
enddo end do
enddo end do
endif; endif end if; end if
!* FLUX FROM ME TO MY NEIGHBOR !* FLUX FROM ME TO MY NEIGHBOR
@ -1330,20 +1330,20 @@ function rhoDotFlux(timestep,ph,en,ip,el)
transmissivity = sum(compatibility(c,:,s,n,ip,el)**2.0_pReal) ! overall transmissivity from this slip system to my neighbor transmissivity = sum(compatibility(c,:,s,n,ip,el)**2.0_pReal) ! overall transmissivity from this slip system to my neighbor
else ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor else ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor
transmissivity = 0.0_pReal transmissivity = 0.0_pReal
endif end if
lineLength = my_rhoSgl0(s,t) * v0(s,t) & lineLength = my_rhoSgl0(s,t) * v0(s,t) &
* math_inner(m(1:3,s,t), normal_me2neighbor) * area ! positive line length of mobiles that wants to leave through this interface * math_inner(m(1:3,s,t), normal_me2neighbor) * area ! positive line length of mobiles that wants to leave through this interface
rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / IPvolume(ip,el) ! subtract dislocation flux from current type rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / IPvolume(ip,el) ! subtract dislocation flux from current type
rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) & rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) &
+ lineLength / IPvolume(ip,el) * (1.0_pReal - transmissivity) & + lineLength / IPvolume(ip,el) * (1.0_pReal - transmissivity) &
* sign(1.0_pReal, v0(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point * sign(1.0_pReal, v0(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point
endif end if
enddo end do
enddo end do
endif; endif end if; end if
enddo neighbors end do neighbors
endif end if
end associate end associate
@ -1433,7 +1433,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
mis%rotate(prm%slip_normal(1:3,s2)))) & mis%rotate(prm%slip_normal(1:3,s2)))) &
* abs(math_inner(prm%slip_direction(1:3,s1), & * abs(math_inner(prm%slip_direction(1:3,s1), &
mis%rotate(prm%slip_direction(1:3,s2)))) mis%rotate(prm%slip_direction(1:3,s2))))
enddo neighborSlipSystems end do neighborSlipSystems
my_compatibilitySum = 0.0_pReal my_compatibilitySum = 0.0_pReal
belowThreshold = .true. belowThreshold = .true.
@ -1446,15 +1446,15 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
my_compatibility(:,:,s1,n) = sign((1.0_pReal - my_compatibilitySum)/nThresholdValues,& my_compatibility(:,:,s1,n) = sign((1.0_pReal - my_compatibilitySum)/nThresholdValues,&
my_compatibility(:,:,s1,n)) my_compatibility(:,:,s1,n))
my_compatibilitySum = my_compatibilitySum + nThresholdValues * thresholdValue my_compatibilitySum = my_compatibilitySum + nThresholdValues * thresholdValue
enddo end do
where(belowThreshold) my_compatibility(1,:,s1,n) = 0.0_pReal where(belowThreshold) my_compatibility(1,:,s1,n) = 0.0_pReal
where(belowThreshold) my_compatibility(2,:,s1,n) = 0.0_pReal where(belowThreshold) my_compatibility(2,:,s1,n) = 0.0_pReal
enddo mySlipSystems end do mySlipSystems
endif end if
enddo neighbors end do neighbors
compatibility(:,:,:,:,i,e) = my_compatibility compatibility(:,:,:,:,i,e) = my_compatibility
@ -1532,7 +1532,7 @@ module subroutine plastic_nonlocal_results(ph,group)
'passing stress for slip','Pa', prm%systems_sl) 'passing stress for slip','Pa', prm%systems_sl)
end select end select
enddo end do
end associate end associate
@ -1581,7 +1581,7 @@ subroutine stateInit(ini,phase,Nentries)
s = nint(rnd(2)*real(sum(ini%N_sl),pReal)*4.0_pReal + 0.5_pReal) s = nint(rnd(2)*real(sum(ini%N_sl),pReal)*4.0_pReal + 0.5_pReal)
meanDensity = meanDensity + densityBinning * geom(phase)%V_0(e) / totalVolume meanDensity = meanDensity + densityBinning * geom(phase)%V_0(e) / totalVolume
stt%rhoSglMobile(s,e) = densityBinning stt%rhoSglMobile(s,e) = densityBinning
enddo end do
else ! homogeneous distribution with noise else ! homogeneous distribution with noise
do e = 1, Nentries do e = 1, Nentries
do f = 1,size(ini%N_sl,1) do f = 1,size(ini%N_sl,1)
@ -1594,12 +1594,12 @@ subroutine stateInit(ini,phase,Nentries)
stt%rho_sgl_mob_edg_neg(s,e) = ini%rho_u_ed_neg_0(f) + noise(1) stt%rho_sgl_mob_edg_neg(s,e) = ini%rho_u_ed_neg_0(f) + noise(1)
stt%rho_sgl_mob_scr_pos(s,e) = ini%rho_u_sc_pos_0(f) + noise(2) stt%rho_sgl_mob_scr_pos(s,e) = ini%rho_u_sc_pos_0(f) + noise(2)
stt%rho_sgl_mob_scr_neg(s,e) = ini%rho_u_sc_neg_0(f) + noise(2) stt%rho_sgl_mob_scr_neg(s,e) = ini%rho_u_sc_neg_0(f) + noise(2)
enddo end do
stt%rho_dip_edg(from:upto,e) = ini%rho_d_ed_0(f) stt%rho_dip_edg(from:upto,e) = ini%rho_d_ed_0(f)
stt%rho_dip_scr(from:upto,e) = ini%rho_d_sc_0(f) stt%rho_dip_scr(from:upto,e) = ini%rho_d_sc_0(f)
enddo end do
enddo end do
endif end if
end associate end associate
@ -1688,8 +1688,8 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, T,
dv_dtau(s) = v(s)**2.0_pReal * (dtSolidSolution_dtau / lambda_S + prm%B / (prm%b_sl(s) * tauEff**2.0_pReal)) dv_dtau(s) = v(s)**2.0_pReal * (dtSolidSolution_dtau / lambda_S + prm%B / (prm%b_sl(s) * tauEff**2.0_pReal))
dv_dtauNS(s) = v(s)**2.0_pReal * dtPeierls_dtau / lambda_P dv_dtauNS(s) = v(s)**2.0_pReal * dtPeierls_dtau / lambda_P
endif end if
enddo end do
end associate end associate
@ -1760,8 +1760,8 @@ subroutine storeGeometry(ph)
do ce = 1, size(material_homogenizationEntry,1) do ce = 1, size(material_homogenizationEntry,1)
do co = 1, homogenization_maxNconstituents do co = 1, homogenization_maxNconstituents
if (material_phaseID(co,ce) == ph) geom(ph)%V_0(material_phaseEntry(co,ce)) = V(ce) if (material_phaseID(co,ce) == ph) geom(ph)%V_0(material_phaseEntry(co,ce)) = V(ce)
enddo end do
enddo end do
end subroutine end subroutine

30
src/phase_mechanical_plastic_phenopowerlaw.f90
View File

@ -95,8 +95,8 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
myPlasticity = plastic_active('phenopowerlaw') myPlasticity = plastic_active('phenopowerlaw')
if(count(myPlasticity) == 0) return if(count(myPlasticity) == 0) return
print'(/,a)', ' <<<+- phase:mechanical:plastic:phenopowerlaw init -+>>>' print'(/,1x,a)', '<<<+- phase:mechanical:plastic:phenopowerlaw init -+>>>'
print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
phases => config_material%get('phase') phases => config_material%get('phase')
@ -105,7 +105,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
allocate(dotState(phases%length)) allocate(dotState(phases%length))
do ph = 1, phases%length do ph = 1, phases%length
if(.not. myPlasticity(ph)) cycle if (.not. myPlasticity(ph)) cycle
associate(prm => param(ph), dot => dotState(ph), stt => state(ph)) associate(prm => param(ph), dot => dotState(ph), stt => state(ph))
@ -129,7 +129,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
else else
prm%P_nS_pos = prm%P_sl prm%P_nS_pos = prm%P_sl
prm%P_nS_neg = prm%P_sl prm%P_nS_neg = prm%P_sl
endif end if
prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dFloat('h_sl-sl'),phase_lattice(ph)) prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dFloat('h_sl-sl'),phase_lattice(ph))
xi_0_sl = pl%get_as1dFloat('xi_0_sl', requiredSize=size(N_sl)) xi_0_sl = pl%get_as1dFloat('xi_0_sl', requiredSize=size(N_sl))
@ -158,7 +158,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
xi_0_sl = emptyRealArray xi_0_sl = emptyRealArray
allocate(prm%xi_inf_sl,prm%h_int,source=emptyRealArray) allocate(prm%xi_inf_sl,prm%h_int,source=emptyRealArray)
allocate(prm%h_sl_sl(0,0)) allocate(prm%h_sl_sl(0,0))
endif slipActive end if slipActive
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! twin related parameters ! twin related parameters
@ -192,7 +192,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
xi_0_tw = emptyRealArray xi_0_tw = emptyRealArray
allocate(prm%gamma_char,source=emptyRealArray) allocate(prm%gamma_char,source=emptyRealArray)
allocate(prm%h_tw_tw(0,0)) allocate(prm%h_tw_tw(0,0))
endif twinActive end if twinActive
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! slip-twin related parameters ! slip-twin related parameters
@ -206,7 +206,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
allocate(prm%h_sl_tw(prm%sum_N_sl,prm%sum_N_tw)) ! at least one dimension is 0 allocate(prm%h_sl_tw(prm%sum_N_sl,prm%sum_N_tw)) ! at least one dimension is 0
allocate(prm%h_tw_sl(prm%sum_N_tw,prm%sum_N_sl)) ! at least one dimension is 0 allocate(prm%h_tw_sl(prm%sum_N_tw,prm%sum_N_sl)) ! at least one dimension is 0
prm%h_0_tw_sl = 0.0_pReal prm%h_0_tw_sl = 0.0_pReal
endif slipAndTwinActive end if slipAndTwinActive
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! output pararameters ! output pararameters
@ -263,7 +263,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
! exit if any parameter is out of range ! exit if any parameter is out of range
if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(phenopowerlaw)') if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(phenopowerlaw)')
enddo end do
end function plastic_phenopowerlaw_init end function plastic_phenopowerlaw_init
@ -306,7 +306,7 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ ddot_gamma_dtau_sl_pos(i) * prm%P_sl(k,l,i) * prm%P_nS_pos(m,n,i) & + ddot_gamma_dtau_sl_pos(i) * prm%P_sl(k,l,i) * prm%P_nS_pos(m,n,i) &
+ ddot_gamma_dtau_sl_neg(i) * prm%P_sl(k,l,i) * prm%P_nS_neg(m,n,i) + ddot_gamma_dtau_sl_neg(i) * prm%P_sl(k,l,i) * prm%P_nS_neg(m,n,i)
enddo slipSystems end do slipSystems
call kinetics_tw(Mp,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw) call kinetics_tw(Mp,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw)
twinSystems: do i = 1, prm%sum_N_tw twinSystems: do i = 1, prm%sum_N_tw
@ -314,7 +314,7 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ ddot_gamma_dtau_tw(i)*prm%P_tw(k,l,i)*prm%P_tw(m,n,i) + ddot_gamma_dtau_tw(i)*prm%P_tw(k,l,i)*prm%P_tw(m,n,i)
enddo twinSystems end do twinSystems
end associate end associate
@ -397,7 +397,7 @@ module subroutine plastic_phenopowerlaw_results(ph,group)
end select end select
enddo end do
end associate end associate
@ -438,7 +438,7 @@ pure subroutine kinetics_sl(Mp,ph,en, &
tau_sl_pos(i) = math_tensordot(Mp,prm%P_nS_pos(1:3,1:3,i)) tau_sl_pos(i) = math_tensordot(Mp,prm%P_nS_pos(1:3,1:3,i))
tau_sl_neg(i) = merge(math_tensordot(Mp,prm%P_nS_neg(1:3,1:3,i)), & tau_sl_neg(i) = merge(math_tensordot(Mp,prm%P_nS_neg(1:3,1:3,i)), &
0.0_pReal, prm%nonSchmidActive) 0.0_pReal, prm%nonSchmidActive)
enddo end do
where(dNeq0(tau_sl_pos)) where(dNeq0(tau_sl_pos))
dot_gamma_sl_pos = prm%dot_gamma_0_sl * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active dot_gamma_sl_pos = prm%dot_gamma_0_sl * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active
@ -460,14 +460,14 @@ pure subroutine kinetics_sl(Mp,ph,en, &
else where else where
ddot_gamma_dtau_sl_pos = 0.0_pReal ddot_gamma_dtau_sl_pos = 0.0_pReal
end where end where
endif end if
if (present(ddot_gamma_dtau_sl_neg)) then if (present(ddot_gamma_dtau_sl_neg)) then
where(dNeq0(dot_gamma_sl_neg)) where(dNeq0(dot_gamma_sl_neg))
ddot_gamma_dtau_sl_neg = dot_gamma_sl_neg*prm%n_sl/tau_sl_neg ddot_gamma_dtau_sl_neg = dot_gamma_sl_neg*prm%n_sl/tau_sl_neg
else where else where
ddot_gamma_dtau_sl_neg = 0.0_pReal ddot_gamma_dtau_sl_neg = 0.0_pReal
end where end where
endif end if
end associate end associate
@ -517,7 +517,7 @@ pure subroutine kinetics_tw(Mp,ph,en,&
else where else where
ddot_gamma_dtau_tw = 0.0_pReal ddot_gamma_dtau_tw = 0.0_pReal
end where end where
endif end if
end associate end associate

2
src/phase_thermal.f90
View File

@ -86,7 +86,7 @@ module subroutine thermal_init(phases)
Nmembers Nmembers
print'(/,a)', ' <<<+- phase:thermal init -+>>>' print'(/,1x,a)', '<<<+- phase:thermal init -+>>>'
allocate(current(phases%length)) allocate(current(phases%length))

14
src/phase_thermal_dissipation.f90
View File

@ -36,8 +36,8 @@ module function dissipation_init(source_length) result(mySources)
mySources = thermal_active('dissipation',source_length) mySources = thermal_active('dissipation',source_length)
if(count(mySources) == 0) return if(count(mySources) == 0) return
print'(/,a)', ' <<<+- phase:thermal:dissipation init -+>>>' print'(/,1x,a)', '<<<+- phase:thermal:dissipation init -+>>>'
print'(a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT) print'(/,a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT)
phases => config_material%get('phase') phases => config_material%get('phase')
@ -45,11 +45,11 @@ module function dissipation_init(source_length) result(mySources)
do ph = 1, phases%length do ph = 1, phases%length
phase => phases%get(ph) phase => phases%get(ph)
if(count(mySources(:,ph)) == 0) cycle !ToDo: error if > 1 if (count(mySources(:,ph)) == 0) cycle !ToDo: error if > 1
thermal => phase%get('thermal') thermal => phase%get('thermal')
sources => thermal%get('source') sources => thermal%get('source')
do so = 1, sources%length do so = 1, sources%length
if(mySources(so,ph)) then if (mySources(so,ph)) then
associate(prm => param(ph)) associate(prm => param(ph))
src => sources%get(so) src => sources%get(so)
@ -58,9 +58,9 @@ module function dissipation_init(source_length) result(mySources)
call phase_allocateState(thermalState(ph)%p(so),Nmembers,0,0,0) call phase_allocateState(thermalState(ph)%p(so),Nmembers,0,0,0)
end associate end associate
endif end if
enddo end do
enddo end do
end function dissipation_init end function dissipation_init

16
src/phase_thermal_externalheat.f90
View File

@ -43,8 +43,8 @@ module function externalheat_init(source_length) result(mySources)
mySources = thermal_active('externalheat',source_length) mySources = thermal_active('externalheat',source_length)
if(count(mySources) == 0) return if(count(mySources) == 0) return
print'(/,a)', ' <<<+- phase:thermal:externalheat init -+>>>' print'(/,1x,a)', '<<<+- phase:thermal:externalheat init -+>>>'
print'(a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT) print'(/,a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT)
phases => config_material%get('phase') phases => config_material%get('phase')
@ -53,11 +53,11 @@ module function externalheat_init(source_length) result(mySources)
do ph = 1, phases%length do ph = 1, phases%length
phase => phases%get(ph) phase => phases%get(ph)
if(count(mySources(:,ph)) == 0) cycle if (count(mySources(:,ph)) == 0) cycle
thermal => phase%get('thermal') thermal => phase%get('thermal')
sources => thermal%get('source') sources => thermal%get('source')
do so = 1, sources%length do so = 1, sources%length
if(mySources(so,ph)) then if (mySources(so,ph)) then
source_thermal_externalheat_offset(ph) = so source_thermal_externalheat_offset(ph) = so
associate(prm => param(ph)) associate(prm => param(ph))
src => sources%get(so) src => sources%get(so)
@ -70,9 +70,9 @@ module function externalheat_init(source_length) result(mySources)
Nmembers = count(material_phaseID == ph) Nmembers = count(material_phaseID == ph)
call phase_allocateState(thermalState(ph)%p(so),Nmembers,1,1,0) call phase_allocateState(thermalState(ph)%p(so),Nmembers,1,1,0)
end associate end associate
endif end if
enddo end do
enddo end do
end function externalheat_init end function externalheat_init
@ -125,7 +125,7 @@ module function externalheat_f_T(ph,en) result(f_T)
f_T = prm%f_T(interval ) * (1.0_pReal - frac_time) + & f_T = prm%f_T(interval ) * (1.0_pReal - frac_time) + &
prm%f_T(interval+1) * frac_time ! interpolate heat rate between segment boundaries... prm%f_T(interval+1) * frac_time ! interpolate heat rate between segment boundaries...
! ...or extrapolate if outside of bounds ! ...or extrapolate if outside of bounds
enddo end do
end associate end associate
end function externalheat_f_T end function externalheat_f_T

16
src/prec.f90
View File

@ -69,15 +69,15 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine prec_init subroutine prec_init
print'(/,a)', ' <<<+- prec init -+>>>' print'(/,1x,a)', '<<<+- prec init -+>>>'
print'(a,i3)', ' Size of integer in bit: ',bit_size(0) print'(/,a,i3)', ' integer size / bit: ',bit_size(0)
print'(a,i19)', ' Maximum value: ',huge(0) print'( a,i19)', ' maximum value: ',huge(0)
print'(/,a,i3)', ' Size of float in bit: ',storage_size(0.0_pReal) print'(/,a,i3)', ' float size / bit: ',storage_size(0.0_pReal)
print'(a,e10.3)', ' Maximum value: ',huge(0.0_pReal) print'( a,e10.3)', ' maximum value: ',huge(0.0_pReal)
print'(a,e10.3)', ' Minimum value: ',PREAL_MIN print'( a,e10.3)', ' minimum value: ',PREAL_MIN
print'(a,e10.3)', ' Epsilon value: ',PREAL_EPSILON print'( a,e10.3)', ' epsilon value: ',PREAL_EPSILON
print'(a,i3)', ' Decimal precision: ',precision(0.0_pReal) print'( a,i3)', ' decimal precision: ',precision(0.0_pReal)
call selfTest call selfTest

42
src/results.f90
View File

@ -70,10 +70,10 @@ subroutine results_init(restart)
character(len=:), allocatable :: date character(len=:), allocatable :: date
print'(/,a)', ' <<<+- results init -+>>>'; flush(IO_STDOUT) print'(/,1x,a)', '<<<+- results init -+>>>'; flush(IO_STDOUT)
print*, 'M. Diehl et al., Integrating Materials and Manufacturing Innovation 6(1):83–91, 2017' print'(/,1x,a)', 'M. Diehl et al., Integrating Materials and Manufacturing Innovation 6(1):83–91, 2017'
print*, 'https://doi.org/10.1007/s40192-017-0084-5'//IO_EOL print'( 1x,a)', 'https://doi.org/10.1007/s40192-017-0084-5'
if (.not. restart) then if (.not. restart) then
resultsFile = HDF5_openFile(getSolverJobName()//'.hdf5','w') resultsFile = HDF5_openFile(getSolverJobName()//'.hdf5','w')
@ -98,7 +98,7 @@ subroutine results_init(restart)
call results_closeGroup(results_addGroup('setup')) call results_closeGroup(results_addGroup('setup'))
call results_addAttribute('description','input data used to run the simulation','setup') call results_addAttribute('description','input data used to run the simulation','setup')
call h5gmove_f(resultsFile,'tmp','setup/previous',hdferr) call h5gmove_f(resultsFile,'tmp','setup/previous',hdferr)
endif end if
call results_closeJobFile call results_closeJobFile
@ -222,7 +222,7 @@ subroutine results_addAttribute_str(attrLabel,attrValue,path)
call HDF5_addAttribute(resultsFile,attrLabel, attrValue, path) call HDF5_addAttribute(resultsFile,attrLabel, attrValue, path)
else else
call HDF5_addAttribute(resultsFile,attrLabel, attrValue) call HDF5_addAttribute(resultsFile,attrLabel, attrValue)
endif end if
end subroutine results_addAttribute_str end subroutine results_addAttribute_str
@ -241,7 +241,7 @@ subroutine results_addAttribute_int(attrLabel,attrValue,path)
call HDF5_addAttribute(resultsFile,attrLabel, attrValue, path) call HDF5_addAttribute(resultsFile,attrLabel, attrValue, path)
else else
call HDF5_addAttribute(resultsFile,attrLabel, attrValue) call HDF5_addAttribute(resultsFile,attrLabel, attrValue)
endif end if
end subroutine results_addAttribute_int end subroutine results_addAttribute_int
@ -260,7 +260,7 @@ subroutine results_addAttribute_real(attrLabel,attrValue,path)
call HDF5_addAttribute(resultsFile,attrLabel, attrValue, path) call HDF5_addAttribute(resultsFile,attrLabel, attrValue, path)
else else
call HDF5_addAttribute(resultsFile,attrLabel, attrValue) call HDF5_addAttribute(resultsFile,attrLabel, attrValue)
endif end if
end subroutine results_addAttribute_real end subroutine results_addAttribute_real
@ -279,7 +279,7 @@ subroutine results_addAttribute_str_array(attrLabel,attrValue,path)
call HDF5_addAttribute(resultsFile,attrLabel, attrValue, path) call HDF5_addAttribute(resultsFile,attrLabel, attrValue, path)
else else
call HDF5_addAttribute(resultsFile,attrLabel, attrValue) call HDF5_addAttribute(resultsFile,attrLabel, attrValue)
endif end if
end subroutine results_addAttribute_str_array end subroutine results_addAttribute_str_array
@ -298,7 +298,7 @@ subroutine results_addAttribute_int_array(attrLabel,attrValue,path)
call HDF5_addAttribute(resultsFile,attrLabel, attrValue, path) call HDF5_addAttribute(resultsFile,attrLabel, attrValue, path)
else else
call HDF5_addAttribute(resultsFile,attrLabel, attrValue) call HDF5_addAttribute(resultsFile,attrLabel, attrValue)
endif end if
end subroutine results_addAttribute_int_array end subroutine results_addAttribute_int_array
@ -317,7 +317,7 @@ subroutine results_addAttribute_real_array(attrLabel,attrValue,path)
call HDF5_addAttribute(resultsFile,attrLabel, attrValue, path) call HDF5_addAttribute(resultsFile,attrLabel, attrValue, path)
else else
call HDF5_addAttribute(resultsFile,attrLabel, attrValue) call HDF5_addAttribute(resultsFile,attrLabel, attrValue)
endif end if
end subroutine results_addAttribute_real_array end subroutine results_addAttribute_real_array
@ -390,7 +390,7 @@ subroutine results_writeVectorDataset_real(dataset,group,label,description,SIuni
if (present(systems)) then if (present(systems)) then
if (size(systems)*size(dataset,2) == 0 ) return !ToDo: maybe also implement for other results_write (not sure about scalar) if (size(systems)*size(dataset,2) == 0 ) return !ToDo: maybe also implement for other results_write (not sure about scalar)
endif end if
groupHandle = results_openGroup(group) groupHandle = results_openGroup(group)
call HDF5_write(dataset,groupHandle,label) call HDF5_write(dataset,groupHandle,label)
@ -422,7 +422,7 @@ subroutine results_writeTensorDataset_real(dataset,group,label,description,SIuni
transposed_ = transposed transposed_ = transposed
else else
transposed_ = .true. transposed_ = .true.
endif end if
groupHandle = results_openGroup(group) groupHandle = results_openGroup(group)
if(transposed_) then if(transposed_) then
@ -430,11 +430,11 @@ subroutine results_writeTensorDataset_real(dataset,group,label,description,SIuni
allocate(dataset_transposed,mold=dataset) allocate(dataset_transposed,mold=dataset)
do i=1,size(dataset_transposed,3) do i=1,size(dataset_transposed,3)
dataset_transposed(:,:,i) = transpose(dataset(:,:,i)) dataset_transposed(:,:,i) = transpose(dataset(:,:,i))
enddo end do
call HDF5_write(dataset_transposed,groupHandle,label) call HDF5_write(dataset_transposed,groupHandle,label)
else else
call HDF5_write(dataset,groupHandle,label) call HDF5_write(dataset,groupHandle,label)
endif end if
call executionStamp(group//'/'//label,description,SIunit) call executionStamp(group//'/'//label,description,SIunit)
call HDF5_closeGroup(groupHandle) call HDF5_closeGroup(groupHandle)
@ -456,7 +456,7 @@ subroutine results_writeVectorDataset_int(dataset,group,label,description,SIunit
if (present(systems)) then if (present(systems)) then
if (size(systems)*size(dataset,2) == 0 ) return !ToDo: maybe also implement for other results_write (not sure about scalar) if (size(systems)*size(dataset,2) == 0 ) return !ToDo: maybe also implement for other results_write (not sure about scalar)
endif end if
groupHandle = results_openGroup(group) groupHandle = results_openGroup(group)
call HDF5_write(dataset,groupHandle,label) call HDF5_write(dataset,groupHandle,label)
@ -542,16 +542,16 @@ subroutine results_mapping_phase(ID,entry,label)
do co = 1, size(ID,1) do co = 1, size(ID,1)
do ce = 1, size(ID,2) do ce = 1, size(ID,2)
entryOffset(ID(co,ce),worldrank) = entryOffset(ID(co,ce),worldrank) +1 entryOffset(ID(co,ce),worldrank) = entryOffset(ID(co,ce),worldrank) +1
enddo end do
enddo end do
call MPI_Allreduce(MPI_IN_PLACE,entryOffset,size(entryOffset),MPI_INT,MPI_SUM,MPI_COMM_WORLD,ierr)! get offset at each process call MPI_Allreduce(MPI_IN_PLACE,entryOffset,size(entryOffset),MPI_INT,MPI_SUM,MPI_COMM_WORLD,ierr)! get offset at each process
if(ierr /= 0) error stop 'MPI error' if(ierr /= 0) error stop 'MPI error'
entryOffset(:,worldrank) = sum(entryOffset(:,0:worldrank-1),2) entryOffset(:,worldrank) = sum(entryOffset(:,0:worldrank-1),2)
do co = 1, size(ID,1) do co = 1, size(ID,1)
do ce = 1, size(ID,2) do ce = 1, size(ID,2)
entryGlobal(co,ce) = entry(co,ce) -1 + entryOffset(ID(co,ce),worldrank) entryGlobal(co,ce) = entry(co,ce) -1 + entryOffset(ID(co,ce),worldrank)
enddo end do
enddo end do
#endif #endif
myShape = int([size(ID,1),writeSize(worldrank)], HSIZE_T) myShape = int([size(ID,1),writeSize(worldrank)], HSIZE_T)
@ -694,13 +694,13 @@ subroutine results_mapping_homogenization(ID,entry,label)
entryOffset = 0 entryOffset = 0
do ce = 1, size(ID,1) do ce = 1, size(ID,1)
entryOffset(ID(ce),worldrank) = entryOffset(ID(ce),worldrank) +1 entryOffset(ID(ce),worldrank) = entryOffset(ID(ce),worldrank) +1
enddo end do
call MPI_Allreduce(MPI_IN_PLACE,entryOffset,size(entryOffset),MPI_INT,MPI_SUM,MPI_COMM_WORLD,ierr)! get offset at each process call MPI_Allreduce(MPI_IN_PLACE,entryOffset,size(entryOffset),MPI_INT,MPI_SUM,MPI_COMM_WORLD,ierr)! get offset at each process
if(ierr /= 0) error stop 'MPI error' if(ierr /= 0) error stop 'MPI error'
entryOffset(:,worldrank) = sum(entryOffset(:,0:worldrank-1),2) entryOffset(:,worldrank) = sum(entryOffset(:,0:worldrank-1),2)
do ce = 1, size(ID,1) do ce = 1, size(ID,1)
entryGlobal(ce) = entry(ce) -1 + entryOffset(ID(ce),worldrank) entryGlobal(ce) = entry(ce) -1 + entryOffset(ID(ce),worldrank)
enddo end do
#endif #endif
myShape = int([writeSize(worldrank)], HSIZE_T) myShape = int([writeSize(worldrank)], HSIZE_T)

6
src/rotations.f90
View File

@ -103,10 +103,10 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine rotations_init subroutine rotations_init
print'(/,a)', ' <<<+- rotations init -+>>>'; flush(IO_STDOUT) print'(/,1x,a)', '<<<+- rotations init -+>>>'; flush(IO_STDOUT)
print*, 'D. Rowenhorst et al., Modelling and Simulation in Materials Science and Engineering 23:083501, 2015' print'(/,1x,a)', 'D. Rowenhorst et al., Modelling and Simulation in Materials Science and Engineering 23:083501, 2015'
print*, 'https://doi.org/10.1088/0965-0393/23/8/083501' print'( 1x,a)', 'https://doi.org/10.1088/0965-0393/23/8/083501'
call selfTest call selfTest