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

improved IO functionality 

- naming for file opening follows python
- damage modules do not read from file any more
This commit is contained in:
Martin Diehl 2019-03-08 23:16:56 +01:00
parent 48cfc35996
commit af707c671c
9 changed files with 355 additions and 498 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
src/CPFEM2.f90
View File

@ -87,11 +87,6 @@ end subroutine CPFEM_initAll
!> @brief allocate the arrays defined in module CPFEM and initialize them
!--------------------------------------------------------------------------------------------------
subroutine CPFEM_init
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
use prec, only: &
pInt, pReal
use IO, only: &
@ -136,8 +131,6 @@ subroutine CPFEM_init
integer(HID_T) :: fileHandle, groupPlasticID, groupHomogID
write(6,'(/,a)') ' <<<+- CPFEM init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
flush(6)
! *** restore the last converged values of each essential variable from the binary file
@ -223,9 +216,6 @@ subroutine CPFEM_age()
crystallite_dPdF, &
crystallite_Tstar0_v, &
crystallite_Tstar_v
use IO, only: &
IO_write_jobRealFile, &
IO_warning
use HDF5_utilities, only: &
HDF5_openFile, &
HDF5_closeFile, &

221
src/IO.f90
View File

@ -25,10 +25,8 @@ module IO
IO_recursiveRead, &
IO_open_file_stat, &
IO_open_file, &
IO_open_jobFile_binary, &
IO_write_jobFile, &
IO_write_jobRealFile, &
IO_read_realFile, &
IO_read_intFile, &
IO_isBlank, &
IO_getTag, &
IO_stringPos, &
@ -229,7 +227,6 @@ recursive function IO_recursiveRead(fileName,cnt) result(fileContent)
end function IO_recursiveRead
!--------------------------------------------------------------------------------------------------
!> @brief opens existing file for reading to given unit. Path to file is relative to working
!! directory
@ -250,6 +247,61 @@ subroutine IO_open_file(fileUnit,path)
end subroutine IO_open_file
!--------------------------------------------------------------------------------------------------
!> @brief opens an existing file for reading or a new file for writing. Name is the job name
!> @details replaces an existing file when writing
!--------------------------------------------------------------------------------------------------
integer function IO_open_jobFile_binary(extension,mode)
use DAMASK_interface, only: &
getSolverJobName
implicit none
character(len=*), intent(in) :: extension
character, intent(in), optional :: mode
if (present(mode)) then
IO_open_jobFile_binary = IO_open_binary(trim(getSolverJobName())//'.'//trim(extension),mode)
else
IO_open_jobFile_binary = IO_open_binary(trim(getSolverJobName())//'.'//trim(extension))
endif
end function IO_open_jobFile_binary
!--------------------------------------------------------------------------------------------------
!> @brief opens an existing file for reading or a new file for writing.
!> @details replaces an existing file when writing
!--------------------------------------------------------------------------------------------------
integer function IO_open_binary(fileName,mode)
implicit none
character(len=*), intent(in) :: fileName
character, intent(in), optional :: mode
character :: m
integer :: ierr
if (present(mode)) then
m = mode
else
m = 'r'
endif
if (m == 'w') then
open(newunit=IO_open_binary, file=trim(fileName),&
status='replace',access='stream',action='write',iostat=ierr)
if (ierr /= 0) call IO_error(100,ext_msg='could not open file (w): '//trim(fileName))
elseif(m == 'r') then
open(newunit=IO_open_binary, file=trim(fileName),&
status='old', access='stream',action='read', iostat=ierr)
if (ierr /= 0) call IO_error(100,ext_msg='could not open file (r): '//trim(fileName))
else
call IO_error(100,ext_msg='unknown access mode: '//m)
endif
end function IO_open_binary
!--------------------------------------------------------------------------------------------------
!> @brief opens existing file for reading to given unit. Path to file is relative to working
!! directory
@ -277,7 +329,6 @@ end function IO_open_file_stat
!--------------------------------------------------------------------------------------------------
subroutine IO_open_inputFile(fileUnit,modelName)
use DAMASK_interface, only: &
getSolverJobName, &
inputFileExtension
implicit none
@ -411,80 +462,6 @@ subroutine IO_write_jobFile(fileUnit,ext)
end subroutine IO_write_jobFile
!--------------------------------------------------------------------------------------------------
!> @brief opens binary file containing array of pReal numbers to given unit for writing. File is
!! named after solver job name plus given extension and located in current working directory
!--------------------------------------------------------------------------------------------------
subroutine IO_write_jobRealFile(fileUnit,ext)
use DAMASK_interface, only: &
getSolverJobName
implicit none
integer(pInt), intent(in) :: fileUnit !< file unit
character(len=*), intent(in) :: ext !< extension of file
integer(pInt) :: myStat
character(len=1024) :: path
path = trim(getSolverJobName())//'.'//ext
open(fileUnit,status='replace',form='unformatted',access='direct', &
recl=pReal,iostat=myStat,file=path)
if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=path)
end subroutine IO_write_jobRealFile
!--------------------------------------------------------------------------------------------------
!> @brief opens binary file containing array of pReal numbers to given unit for reading. File is
!! located in current working directory
!--------------------------------------------------------------------------------------------------
subroutine IO_read_realFile(fileUnit,ext,modelName)
implicit none
integer(pInt), intent(in) :: fileUnit !< file unit
character(len=*), intent(in) :: ext, & !< extension of file
modelName !< model name, in case of restart not solver job name
integer(pInt) :: myStat
character(len=1024) :: path
path = trim(modelName)//'.'//ext
open(fileUnit,status='old',form='unformatted',access='direct', &
recl=pReal,iostat=myStat,file=path)
if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=path)
end subroutine IO_read_realFile
!--------------------------------------------------------------------------------------------------
!> @brief opens binary file containing array of pInt numbers to given unit for reading. File is
!! located in current working directory
!--------------------------------------------------------------------------------------------------
subroutine IO_read_intFile(fileUnit,ext,modelName)
implicit none
integer(pInt), intent(in) :: fileUnit !< file unit
character(len=*), intent(in) :: ext, & !< extension of file
modelName !< model name, in case of restart not solver job name
integer(pInt) :: myStat
character(len=1024) :: path
path = trim(modelName)//'.'//ext
open(fileUnit,status='old',form='unformatted',access='direct', &
recl=pInt,iostat=myStat,file=path)
if (myStat /= 0) call IO_error(100_pInt,ext_msg=path)
end subroutine IO_read_intFile
!--------------------------------------------------------------------------------------------------
!> @brief identifies strings without content
!--------------------------------------------------------------------------------------------------
@ -1401,27 +1378,27 @@ function IO_continuousIntValues(fileUnit,maxN,lookupName,lookupMap,lookupMaxN)
!> @brief returns verified integer value in given string
!--------------------------------------------------------------------------------------------------
integer(pInt) function IO_verifyIntValue (string,validChars,myName)
implicit none
character(len=*), intent(in) :: string, & !< string for conversion to int value. Must not contain spaces!
validChars, & !< valid characters in string
myName !< name of caller function (for debugging)
integer(pInt) :: readStatus, invalidWhere
IO_verifyIntValue = 0_pInt
invalidWhere = verify(string,validChars)
if (invalidWhere == 0_pInt) then
read(UNIT=string,iostat=readStatus,FMT=*) IO_verifyIntValue ! no offending chars found
if (readStatus /= 0_pInt) & ! error during string to integer conversion
call IO_warning(203_pInt,ext_msg=myName//'"'//string//'"')
else
call IO_warning(202_pInt,ext_msg=myName//'"'//string//'"') ! complain about offending characters
read(UNIT=string(1_pInt:invalidWhere-1_pInt),iostat=readStatus,FMT=*) IO_verifyIntValue ! interpret remaining string
if (readStatus /= 0_pInt) & ! error during string to integer conversion
call IO_warning(203_pInt,ext_msg=myName//'"'//string(1_pInt:invalidWhere-1_pInt)//'"')
endif
implicit none
character(len=*), intent(in) :: string, & !< string for conversion to int value. Must not contain spaces!
validChars, & !< valid characters in string
myName !< name of caller function (for debugging)
integer :: readStatus, invalidWhere
IO_verifyIntValue = 0
invalidWhere = verify(string,validChars)
if (invalidWhere == 0) then
read(UNIT=string,iostat=readStatus,FMT=*) IO_verifyIntValue ! no offending chars found
if (readStatus /= 0) & ! error during string to integer conversion
call IO_warning(203,ext_msg=myName//'"'//string//'"')
else
call IO_warning(202,ext_msg=myName//'"'//string//'"') ! complain about offending characters
read(UNIT=string(1:invalidWhere-1),iostat=readStatus,FMT=*) IO_verifyIntValue ! interpret remaining string
if (readStatus /= 0) & ! error during string to integer conversion
call IO_warning(203,ext_msg=myName//'"'//string(1:invalidWhere-1)//'"')
endif
end function IO_verifyIntValue
@ -1429,28 +1406,28 @@ end function IO_verifyIntValue
!> @brief returns verified float value in given string
!--------------------------------------------------------------------------------------------------
real(pReal) function IO_verifyFloatValue (string,validChars,myName)
implicit none
character(len=*), intent(in) :: string, & !< string for conversion to int value. Must not contain spaces!
validChars, & !< valid characters in string
myName !< name of caller function (for debugging)
integer(pInt) :: readStatus, invalidWhere
IO_verifyFloatValue = 0.0_pReal
invalidWhere = verify(string,validChars)
if (invalidWhere == 0_pInt) then
read(UNIT=string,iostat=readStatus,FMT=*) IO_verifyFloatValue ! no offending chars found
if (readStatus /= 0_pInt) & ! error during string to float conversion
call IO_warning(203_pInt,ext_msg=myName//'"'//string//'"')
else
call IO_warning(202_pInt,ext_msg=myName//'"'//string//'"') ! complain about offending characters
read(UNIT=string(1_pInt:invalidWhere-1_pInt),iostat=readStatus,FMT=*) IO_verifyFloatValue ! interpret remaining string
if (readStatus /= 0_pInt) & ! error during string to float conversion
call IO_warning(203_pInt,ext_msg=myName//'"'//string(1_pInt:invalidWhere-1_pInt)//'"')
endif
implicit none
character(len=*), intent(in) :: string, & !< string for conversion to int value. Must not contain spaces!
validChars, & !< valid characters in string
myName !< name of caller function (for debugging)
integer :: readStatus, invalidWhere
IO_verifyFloatValue = 0.0_pReal
invalidWhere = verify(string,validChars)
if (invalidWhere == 0) then
read(UNIT=string,iostat=readStatus,FMT=*) IO_verifyFloatValue ! no offending chars found
if (readStatus /= 0) & ! error during string to float conversion
call IO_warning(203,ext_msg=myName//'"'//string//'"')
else
call IO_warning(202,ext_msg=myName//'"'//string//'"') ! complain about offending characters
read(UNIT=string(1:invalidWhere-1),iostat=readStatus,FMT=*) IO_verifyFloatValue ! interpret remaining string
if (readStatus /= 0) & ! error during string to float conversion
call IO_warning(203,ext_msg=myName//'"'//string(1:invalidWhere-1)//'"')
endif
end function IO_verifyFloatValue
end module IO

152
src/damage_local.f90
View File

@ -9,9 +9,6 @@ module damage_local
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
damage_local_sizePostResults !< cumulative size of post results
integer(pInt), dimension(:,:), allocatable, target, public :: &
damage_local_sizePostResult !< size of each post result output
@ -27,7 +24,15 @@ module damage_local
end enum
integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: &
damage_local_outputID !< ID of each post result output
type, private :: tParameters
integer(kind(undefined_ID)), dimension(:), allocatable :: &
outputID
end type tParameters
type(tparameters), dimension(:), allocatable, private :: &
param
public :: &
damage_local_init, &
damage_local_updateState, &
@ -38,27 +43,10 @@ module damage_local
contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields, reads information from material configuration file
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine damage_local_init(fileUnit)
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
use IO, only: &
IO_read, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_intValue, &
IO_warning, &
IO_error, &
IO_timeStamp, &
IO_EOF
subroutine damage_local_init
use material, only: &
damage_type, &
damage_typeInstance, &
@ -72,94 +60,65 @@ subroutine damage_local_init(fileUnit)
damage, &
damage_initialPhi
use config, only: &
material_partHomogenization
config_homogenization
implicit none
integer(pInt), intent(in) :: fileUnit
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: maxNinstance,mySize=0_pInt,homog,instance,o
integer(pInt) :: maxNinstance,homog,instance,o,i
integer(pInt) :: sizeState
integer(pInt) :: NofMyHomog
character(len=65536) :: &
tag = '', &
line = ''
integer(pInt) :: NofMyHomog, h
integer(kind(undefined_ID)) :: &
outputID
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
character(len=65536), dimension(:), allocatable :: &
outputs
write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_local_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
maxNinstance = int(count(damage_type == DAMAGE_local_ID),pInt)
if (maxNinstance == 0_pInt) return
allocate(damage_local_sizePostResults(maxNinstance), source=0_pInt)
allocate(damage_local_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0_pInt)
allocate(damage_local_output (maxval(homogenization_Noutput),maxNinstance))
damage_local_output = ''
allocate(damage_local_outputID (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID)
allocate(damage_local_Noutput (maxNinstance), source=0_pInt)
rewind(fileUnit)
homog = 0_pInt
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization)! wind forward to <homogenization>
line = IO_read(fileUnit)
enddo
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of homog part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next homog section
homog = homog + 1_pInt ! advance homog section counter
cycle ! skip to next line
endif
allocate(param(maxNinstance))
do h = 1, size(damage_type)
if (damage_type(h) /= DAMAGE_LOCAL_ID) cycle
associate(prm => param(damage_typeInstance(h)), &
config => config_homogenization(h))
if (homog > 0_pInt ) then; if (damage_type(homog) == DAMAGE_local_ID) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
instance = damage_typeInstance(homog) ! which instance of my damage is present homog
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag)
case ('(output)')
select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))
case ('damage')
damage_local_Noutput(instance) = damage_local_Noutput(instance) + 1_pInt
damage_local_outputID(damage_local_Noutput(instance),instance) = damage_ID
damage_local_output(damage_local_Noutput(instance),instance) = &
IO_lc(IO_stringValue(line,chunkPos,2_pInt))
end select
end select
endif; endif
enddo parsingFile
initializeInstances: do homog = 1_pInt, size(damage_type)
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
myhomog: if (damage_type(homog) == DAMAGE_local_ID) then
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
case ('damage')
damage_local_output(i,damage_typeInstance(h)) = outputs(i)
damage_local_Noutput(instance) = damage_local_Noutput(instance) + 1
damage_local_sizePostResult(i,damage_typeInstance(h)) = 1
prm%outputID = [prm%outputID , damage_ID]
end select
enddo
homog = h
NofMyHomog = count(material_homog == homog)
instance = damage_typeInstance(homog)
!--------------------------------------------------------------------------------------------------
! Determine size of postResults array
outputsLoop: do o = 1_pInt,damage_local_Noutput(instance)
select case(damage_local_outputID(o,instance))
case(damage_ID)
mySize = 1_pInt
end select
if (mySize > 0_pInt) then ! any meaningful output found
damage_local_sizePostResult(o,instance) = mySize
damage_local_sizePostResults(instance) = damage_local_sizePostResults(instance) + mySize
endif
enddo outputsLoop
! allocate state arrays
sizeState = 1_pInt
damageState(homog)%sizeState = sizeState
damageState(homog)%sizePostResults = damage_local_sizePostResults(instance)
damageState(homog)%sizePostResults = sum(damage_local_sizePostResult(:,instance))
allocate(damageState(homog)%state0 (sizeState,NofMyHomog), source=damage_initialPhi(homog))
allocate(damageState(homog)%subState0(sizeState,NofMyHomog), source=damage_initialPhi(homog))
allocate(damageState(homog)%state (sizeState,NofMyHomog), source=damage_initialPhi(homog))
@ -169,8 +128,8 @@ subroutine damage_local_init(fileUnit)
deallocate(damage(homog)%p)
damage(homog)%p => damageState(homog)%state(1,:)
endif myhomog
enddo initializeInstances
end associate
enddo
end subroutine damage_local_init
@ -193,7 +152,7 @@ function damage_local_updateState(subdt, ip, el)
el !< element number
real(pReal), intent(in) :: &
subdt
logical, dimension(2) :: &
logical, dimension(2) :: &
damage_local_updateState
integer(pInt) :: &
homog, &
@ -303,7 +262,7 @@ function damage_local_postResults(ip,el)
integer(pInt), intent(in) :: &
ip, & !< integration point
el !< element
real(pReal), dimension(damage_local_sizePostResults(damage_typeInstance(mappingHomogenization(2,ip,el)))) :: &
real(pReal), dimension(sum(damage_local_sizePostResult(:,damage_typeInstance(mappingHomogenization(2,ip,el))))) :: &
damage_local_postResults
integer(pInt) :: &
@ -312,18 +271,19 @@ function damage_local_postResults(ip,el)
homog = mappingHomogenization(2,ip,el)
offset = damageMapping(homog)%p(ip,el)
instance = damage_typeInstance(homog)
associate(prm => param(instance))
c = 0_pInt
damage_local_postResults = 0.0_pReal
do o = 1_pInt,damage_local_Noutput(instance)
select case(damage_local_outputID(o,instance))
outputsLoop: do o = 1_pInt,size(prm%outputID)
select case(prm%outputID(o))
case (damage_ID)
damage_local_postResults(c+1_pInt) = damage(homog)%p(offset)
c = c + 1
end select
enddo
enddo outputsLoop
end associate
end function damage_local_postResults
end module damage_local

171
src/damage_nonlocal.f90
View File

@ -10,9 +10,6 @@ module damage_nonlocal
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
damage_nonlocal_sizePostResults !< cumulative size of post results
integer(pInt), dimension(:,:), allocatable, target, public :: &
damage_nonlocal_sizePostResult !< size of each post result output
@ -26,9 +23,14 @@ module damage_nonlocal
enumerator :: undefined_ID, &
damage_ID
end enum
integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: &
damage_nonlocal_outputID !< ID of each post result output
type, private :: tParameters
integer(kind(undefined_ID)), dimension(:), allocatable :: &
outputID
end type tParameters
type(tparameters), dimension(:), allocatable, private :: &
param
public :: &
damage_nonlocal_init, &
@ -40,30 +42,11 @@ module damage_nonlocal
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine damage_nonlocal_init(fileUnit)
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
use IO, only: &
IO_read, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_intValue, &
IO_warning, &
IO_error, &
IO_timeStamp, &
IO_EOF
subroutine damage_nonlocal_init
use material, only: &
damage_type, &
damage_typeInstance, &
@ -77,105 +60,75 @@ subroutine damage_nonlocal_init(fileUnit)
damage, &
damage_initialPhi
use config, only: &
material_partHomogenization
config_homogenization
implicit none
integer(pInt), intent(in) :: fileUnit
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: maxNinstance,mySize=0_pInt,section,instance,o
integer(pInt) :: maxNinstance,homog,instance,o,i
integer(pInt) :: sizeState
integer(pInt) :: NofMyHomog
character(len=65536) :: &
tag = '', &
line = ''
integer(pInt) :: NofMyHomog, h
integer(kind(undefined_ID)) :: &
outputID
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
character(len=65536), dimension(:), allocatable :: &
outputs
write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_nonlocal_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
maxNinstance = int(count(damage_type == DAMAGE_nonlocal_ID),pInt)
if (maxNinstance == 0_pInt) return
allocate(damage_nonlocal_sizePostResults(maxNinstance), source=0_pInt)
allocate(damage_nonlocal_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0_pInt)
allocate(damage_nonlocal_output (maxval(homogenization_Noutput),maxNinstance))
damage_nonlocal_output = ''
allocate(damage_nonlocal_outputID (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID)
allocate(damage_nonlocal_Noutput (maxNinstance), source=0_pInt)
rewind(fileUnit)
section = 0_pInt
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization)! wind forward to <homogenization>
line = IO_read(fileUnit)
enddo
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of homog part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next homog section
section = section + 1_pInt ! advance homog section counter
cycle ! skip to next line
endif
allocate(param(maxNinstance))
do h = 1, size(damage_type)
if (damage_type(h) /= DAMAGE_NONLOCAL_ID) cycle
associate(prm => param(damage_typeInstance(h)), &
config => config_homogenization(h))
if (section > 0_pInt ) then; if (damage_type(section) == DAMAGE_nonlocal_ID) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
instance = damage_typeInstance(section) ! which instance of my damage is present homog
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag)
case ('(output)')
select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
case ('damage')
damage_nonlocal_Noutput(instance) = damage_nonlocal_Noutput(instance) + 1_pInt
damage_nonlocal_outputID(damage_nonlocal_Noutput(instance),instance) = damage_ID
damage_nonlocal_output(damage_nonlocal_Noutput(instance),instance) = &
IO_lc(IO_stringValue(line,chunkPos,2_pInt))
damage_nonlocal_output(i,damage_typeInstance(h)) = outputs(i)
damage_nonlocal_Noutput(instance) = damage_nonlocal_Noutput(instance) + 1
damage_nonlocal_sizePostResult(i,damage_typeInstance(h)) = 1
prm%outputID = [prm%outputID , damage_ID]
end select
enddo
end select
endif; endif
enddo parsingFile
initializeInstances: do section = 1_pInt, size(damage_type)
if (damage_type(section) == DAMAGE_nonlocal_ID) then
NofMyHomog=count(material_homog==section)
instance = damage_typeInstance(section)
!--------------------------------------------------------------------------------------------------
! Determine size of postResults array
outputsLoop: do o = 1_pInt,damage_nonlocal_Noutput(instance)
select case(damage_nonlocal_outputID(o,instance))
case(damage_ID)
mySize = 1_pInt
end select
if (mySize > 0_pInt) then ! any meaningful output found
damage_nonlocal_sizePostResult(o,instance) = mySize
damage_nonlocal_sizePostResults(instance) = damage_nonlocal_sizePostResults(instance) + mySize
endif
enddo outputsLoop
homog = h
NofMyHomog = count(material_homog == homog)
instance = damage_typeInstance(homog)
! allocate state arrays
sizeState = 0_pInt
damageState(section)%sizeState = sizeState
damageState(section)%sizePostResults = damage_nonlocal_sizePostResults(instance)
allocate(damageState(section)%state0 (sizeState,NofMyHomog))
allocate(damageState(section)%subState0(sizeState,NofMyHomog))
allocate(damageState(section)%state (sizeState,NofMyHomog))
sizeState = 1_pInt
damageState(homog)%sizeState = sizeState
damageState(homog)%sizePostResults = sum(damage_nonlocal_sizePostResult(:,instance))
allocate(damageState(homog)%state0 (sizeState,NofMyHomog), source=damage_initialPhi(homog))
allocate(damageState(homog)%subState0(sizeState,NofMyHomog), source=damage_initialPhi(homog))
allocate(damageState(homog)%state (sizeState,NofMyHomog), source=damage_initialPhi(homog))
nullify(damageMapping(section)%p)
damageMapping(section)%p => mappingHomogenization(1,:,:)
deallocate(damage(section)%p)
allocate(damage(section)%p(NofMyHomog), source=damage_initialPhi(section))
nullify(damageMapping(homog)%p)
damageMapping(homog)%p => mappingHomogenization(1,:,:)
deallocate(damage(homog)%p)
damage(homog)%p => damageState(homog)%state(1,:)
endif
enddo initializeInstances
end associate
enddo
end subroutine damage_nonlocal_init
!--------------------------------------------------------------------------------------------------
@ -221,7 +174,7 @@ subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip,
do grain = 1, homogenization_Ngrains(mappingHomogenization(2,ip,el))
phase = phaseAt(grain,ip,el)
constituent = phasememberAt(grain,ip,el)
do source = 1_pInt, phase_Nsources(phase)
do source = 1, phase_Nsources(phase)
select case(phase_source(source,phase))
case (SOURCE_damage_isoBrittle_ID)
call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
@ -349,33 +302,35 @@ function damage_nonlocal_postResults(ip,el)
use material, only: &
mappingHomogenization, &
damage_typeInstance, &
damageMapping, &
damage
implicit none
integer(pInt), intent(in) :: &
ip, & !< integration point
el !< element
real(pReal), dimension(damage_nonlocal_sizePostResults(damage_typeInstance(mappingHomogenization(2,ip,el)))) :: &
real(pReal), dimension(sum(damage_nonlocal_sizePostResult(:,damage_typeInstance(mappingHomogenization(2,ip,el))))) :: &
damage_nonlocal_postResults
integer(pInt) :: &
instance, homog, offset, o, c
homog = mappingHomogenization(2,ip,el)
offset = mappingHomogenization(1,ip,el)
offset = damageMapping(homog)%p(ip,el)
instance = damage_typeInstance(homog)
associate(prm => param(instance))
c = 0_pInt
damage_nonlocal_postResults = 0.0_pReal
do o = 1_pInt,damage_nonlocal_Noutput(instance)
select case(damage_nonlocal_outputID(o,instance))
outputsLoop: do o = 1_pInt,size(prm%outputID)
select case(prm%outputID(o))
case (damage_ID)
damage_nonlocal_postResults(c+1_pInt) = damage(homog)%p(offset)
c = c + 1
end select
enddo
enddo outputsLoop
end associate
end function damage_nonlocal_postResults
end module damage_nonlocal

22
src/homogenization.f90
View File

@ -57,11 +57,6 @@ contains
!> @brief module initialization
!--------------------------------------------------------------------------------------------------
subroutine homogenization_init
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
use math, only: &
math_I3
use debug, only: &
@ -79,8 +74,6 @@ subroutine homogenization_init
use crystallite, only: &
crystallite_maxSizePostResults
use config, only: &
material_configFile, &
material_localFileExt, &
config_deallocate, &
config_homogenization, &
homogenization_name
@ -116,16 +109,9 @@ subroutine homogenization_init
if (any(thermal_type == THERMAL_adiabatic_ID)) call thermal_adiabatic_init
if (any(thermal_type == THERMAL_conduction_ID)) call thermal_conduction_init
!--------------------------------------------------------------------------------------------------
! open material.config
if (.not. IO_open_jobFile_stat(FILEUNIT,material_localFileExt)) & ! no local material configuration present...
call IO_open_file(FILEUNIT,material_configFile) ! ... open material.config file
if (any(damage_type == DAMAGE_none_ID)) &
call damage_none_init()
if (any(damage_type == DAMAGE_local_ID)) &
call damage_local_init(FILEUNIT)
if (any(damage_type == DAMAGE_nonlocal_ID)) &
call damage_nonlocal_init(FILEUNIT)
if (any(damage_type == DAMAGE_none_ID)) call damage_none_init
if (any(damage_type == DAMAGE_local_ID)) call damage_local_init
if (any(damage_type == DAMAGE_nonlocal_ID)) call damage_nonlocal_init
!--------------------------------------------------------------------------------------------------
! write description file for homogenization output
@ -265,8 +251,6 @@ subroutine homogenization_init
allocate(materialpoint_results(materialpoint_sizeResults,theMesh%elem%nIPs,theMesh%nElems))
write(6,'(/,a)') ' <<<+- homogenization init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt) then
#ifdef TODO

11
src/spectral_damage.f90
View File

@ -57,15 +57,8 @@ contains
!> @brief allocates all neccessary fields and fills them with data, potentially from restart info
!--------------------------------------------------------------------------------------------------
subroutine spectral_damage_init()
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
use IO, only: &
IO_intOut, &
IO_read_realFile, &
IO_timeStamp
IO_intOut
use spectral_utilities, only: &
wgt
use mesh, only: &
@ -87,8 +80,6 @@ subroutine spectral_damage_init()
write(6,'(/,a)') ' <<<+- spectral_damage init -+>>>'
write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, volume in press, '
write(6,'(a,/)') ' chapter Spectral Solvers for Crystal Plasticity and Multi-Physics Simulations. Springer, 2018 '
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc

69
src/spectral_mech_Basic.f90
View File

@ -73,16 +73,10 @@ contains
!> @brief allocates all necessary fields and fills them with data, potentially from restart info
!--------------------------------------------------------------------------------------------------
subroutine basic_init
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
use IO, only: &
IO_intOut, &
IO_error, &
IO_read_realFile, &
IO_timeStamp
IO_open_jobFile_binary
use debug, only: &
debug_level, &
debug_spectral, &
@ -115,14 +109,12 @@ subroutine basic_init
PetscErrorCode :: ierr
PetscScalar, pointer, dimension(:,:,:,:) :: F
PetscInt, dimension(:), allocatable :: localK
integer(pInt) :: proc
integer :: proc, fileUnit
character(len=1024) :: rankStr
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasic init -+>>>'
write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity, 66:31–45, 2015'
write(6,'(a,/)') ' https://doi.org/10.1016/j.ijplas.2014.02.006'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
!--------------------------------------------------------------------------------------------------
! allocate global fields
@ -134,7 +126,7 @@ subroutine basic_init
call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr)
allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3
do proc = 1, worldsize
do proc = 1, worldsize !ToDo: there are smarter options in MPI
call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr)
enddo
call DMDACreate3d(PETSC_COMM_WORLD, &
@ -166,13 +158,17 @@ subroutine basic_init
'reading values of increment ', restartInc, ' from file'
flush(6)
endif
fileUnit = IO_open_jobFile_binary('F_aimDot')
read(fileUnit) F_aimDot; close(fileUnit)
write(rankStr,'(a1,i0)')'_',worldrank
call IO_read_realFile(777,'F'//trim(rankStr),trim(getSolverJobName()),size(F))
read (777,rec=1) F; close (777)
call IO_read_realFile(777,'F_lastInc'//trim(rankStr),trim(getSolverJobName()),size(F_lastInc))
read (777,rec=1) F_lastInc; close (777)
call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(F_aimDot))
read (777,rec=1) F_aimDot; close (777)
fileUnit = IO_open_jobFile_binary('F'//trim(rankStr))
read(fileUnit) F; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr))
read(fileUnit) F_lastInc; close (fileUnit)
F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F
call MPI_Allreduce(MPI_IN_PLACE,F_aim,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='F_aim')
@ -198,12 +194,12 @@ subroutine basic_init
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') &
'reading more values of increment ', restartInc, ' from file'
flush(6)
call IO_read_realFile(777,'C_volAvg',trim(getSolverJobName()),size(C_volAvg))
read (777,rec=1) C_volAvg; close (777)
call IO_read_realFile(777,'C_volAvgLastInc',trim(getSolverJobName()),size(C_volAvgLastInc))
read (777,rec=1) C_volAvgLastInc; close (777)
call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg))
read (777,rec=1) C_minMaxAvg; close (777)
fileUnit = IO_open_jobFile_binary('C_volAvg')
read(fileUnit) C_volAvg; close(fileUnit)
fileUnit = IO_open_jobFile_binary('C_volAvgLastInv')
read(fileUnit) C_volAvgLastInc; close(fileUnit)
fileUnit = IO_open_jobFile_binary('C_ref')
read(fileUnit) C_minMaxAvg; close(fileUnit)
endif restartRead
call Utilities_updateGamma(C_minMaxAvg,.true.)
@ -450,7 +446,7 @@ subroutine Basic_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,s
tBoundaryCondition, &
cutBack
use IO, only: &
IO_write_JobRealFile
IO_open_jobFile_binary
use FEsolving, only: &
restartWrite
@ -468,7 +464,8 @@ subroutine Basic_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,s
rotation_BC
PetscErrorCode :: ierr
PetscScalar, dimension(:,:,:,:), pointer :: F
integer :: fileUnit
character(len=32) :: rankStr
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
@ -483,20 +480,20 @@ subroutine Basic_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,s
write(6,'(/,a)') ' writing converged results for restart'
flush(6)
if (worldrank == 0_pInt) then
call IO_write_jobRealFile(777,'C_volAvg',size(C_volAvg))
write (777,rec=1) C_volAvg; close(777)
call IO_write_jobRealFile(777,'C_volAvgLastInc',size(C_volAvgLastInc))
write (777,rec=1) C_volAvgLastInc; close(777)
call IO_write_jobRealFile(777,'F_aimDot',size(F_aimDot))
write (777,rec=1) F_aimDot; close(777)
if (worldrank == 0) then
fileUnit = IO_open_jobFile_binary('C_volAvg','w')
write(fileUnit) C_volAvg; close(fileUnit)
fileUnit = IO_open_jobFile_binary('C_volAvgLastInv','w')
write(fileUnit) C_volAvgLastInc; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aimDot','w')
write(fileUnit) F_aimDot; close(fileUnit)
endif
write(rankStr,'(a1,i0)')'_',worldrank
call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file
write (777,rec=1) F; close (777)
call IO_write_jobRealFile(777,'F_lastInc'//trim(rankStr),size(F_lastInc)) ! writing F_lastInc field to file
write (777,rec=1) F_lastInc; close (777)
fileUnit = IO_open_jobFile_binary('F'//trim(rankStr),'w')
write(fileUnit) F; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr),'w')
write(fileUnit) F_lastInc; close (fileUnit)
endif
call CPFEM_age() ! age state and kinematics

89
src/spectral_mech_Polarisation.f90
View File

@ -80,16 +80,10 @@ contains
!> @brief allocates all necessary fields and fills them with data, potentially from restart info
!--------------------------------------------------------------------------------------------------
subroutine Polarisation_init
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
use IO, only: &
IO_intOut, &
IO_error, &
IO_read_realFile, &
IO_timeStamp
IO_open_jobFile_binary
use debug, only: &
debug_level, &
debug_spectral, &
@ -125,14 +119,12 @@ subroutine Polarisation_init
F, & ! specific (sub)pointer
F_tau ! specific (sub)pointer
PetscInt, dimension(:), allocatable :: localK
integer(pInt) :: proc
integer :: proc, fileUnit
character(len=1024) :: rankStr
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverPolarisation init -+>>>'
write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity, 66:31–45, 2015'
write(6,'(a,/)') ' https://doi.org/10.1016/j.ijplas.2014.02.006'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
!--------------------------------------------------------------------------------------------------
! allocate global fields
@ -146,7 +138,7 @@ subroutine Polarisation_init
call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr)
allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3
do proc = 1, worldsize
do proc = 1, worldsize !ToDo: there are smarter options in MPI
call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr)
enddo
call DMDACreate3d(PETSC_COMM_WORLD, &
@ -173,23 +165,28 @@ subroutine Polarisation_init
call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! places pointer on PETSc data
F => FandF_tau( 0: 8,:,:,:)
F_tau => FandF_tau( 9:17,:,:,:)
restart: if (restartInc > 0_pInt) then
restart: if (restartInc > 0) then
if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0) then
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') &
'reading values of increment ', restartInc, ' from file'
flush(6)
endif
fileUnit = IO_open_jobFile_binary('F_aimDot')
read(fileUnit) F_aimDot; close(fileUnit)
write(rankStr,'(a1,i0)')'_',worldrank
call IO_read_realFile(777,'F'//trim(rankStr),trim(getSolverJobName()),size(F))
read (777,rec=1) F; close (777)
call IO_read_realFile(777,'F_lastInc'//trim(rankStr),trim(getSolverJobName()),size(F_lastInc))
read (777,rec=1) F_lastInc; close (777)
call IO_read_realFile(777,'F_tau'//trim(rankStr),trim(getSolverJobName()),size(F_tau))
read (777,rec=1) F_tau; close (777)
call IO_read_realFile(777,'F_tau_lastInc'//trim(rankStr),trim(getSolverJobName()),size(F_tau_lastInc))
read (777,rec=1) F_tau_lastInc; close (777)
call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(F_aimDot))
read (777,rec=1) F_aimDot; close (777)
fileUnit = IO_open_jobFile_binary('F'//trim(rankStr))
read(fileUnit) F; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr))
read(fileUnit) F_lastInc; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_tau'//trim(rankStr))
read(fileUnit) F_tau; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_tau_lastInc'//trim(rankStr))
read(fileUnit) F_tau_lastInc; close (fileUnit)
F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F
call MPI_Allreduce(MPI_IN_PLACE,F_aim,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='F_aim')
@ -218,12 +215,12 @@ subroutine Polarisation_init
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') &
'reading more values of increment ', restartInc, ' from file'
flush(6)
call IO_read_realFile(777,'C_volAvg',trim(getSolverJobName()),size(C_volAvg))
read (777,rec=1) C_volAvg; close (777)
call IO_read_realFile(777,'C_volAvgLastInc',trim(getSolverJobName()),size(C_volAvgLastInc))
read (777,rec=1) C_volAvgLastInc; close (777)
call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg))
read (777,rec=1) C_minMaxAvg; close (777)
fileUnit = IO_open_jobFile_binary('C_volAvg')
read(fileUnit) C_volAvg; close(fileUnit)
fileUnit = IO_open_jobFile_binary('C_volAvgLastInv')
read(fileUnit) C_volAvgLastInc; close(fileUnit)
fileUnit = IO_open_jobFile_binary('C_ref')
read(fileUnit) C_minMaxAvg; close(fileUnit)
endif restartRead
call Utilities_updateGamma(C_minMaxAvg,.true.)
@ -552,7 +549,7 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,deformati
tBoundaryCondition, &
cutBack
use IO, only: &
IO_write_JobRealFile
IO_open_jobFile_binary
use FEsolving, only: &
restartWrite
@ -572,6 +569,8 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,deformati
PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau, F, F_tau
integer(pInt) :: i, j, k
real(pReal), dimension(3,3) :: F_lambda33
integer :: fileUnit
character(len=32) :: rankStr
!--------------------------------------------------------------------------------------------------
@ -590,24 +589,25 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,deformati
write(6,'(/,a)') ' writing converged results for restart'
flush(6)
if (worldrank == 0_pInt) then
call IO_write_jobRealFile(777,'C_volAvg',size(C_volAvg))
write (777,rec=1) C_volAvg; close(777)
call IO_write_jobRealFile(777,'C_volAvgLastInc',size(C_volAvgLastInc))
write (777,rec=1) C_volAvgLastInc; close(777)
call IO_write_jobRealFile(777,'F_aimDot',size(F_aimDot))
write (777,rec=1) F_aimDot; close(777)
if (worldrank == 0) then
fileUnit = IO_open_jobFile_binary('C_volAvg','w')
write(fileUnit) C_volAvg; close(fileUnit)
fileUnit = IO_open_jobFile_binary('C_volAvgLastInv','w')
write(fileUnit) C_volAvgLastInc; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aimDot','w')
write(fileUnit) F_aimDot; close(fileUnit)
endif
write(rankStr,'(a1,i0)')'_',worldrank
call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file
write (777,rec=1) F; close (777)
call IO_write_jobRealFile(777,'F_lastInc'//trim(rankStr),size(F_lastInc)) ! writing F_lastInc field to file
write (777,rec=1) F_lastInc; close (777)
call IO_write_jobRealFile(777,'F_tau'//trim(rankStr),size(F_tau)) ! writing deformation gradient field to file
write (777,rec=1) F_tau; close (777)
call IO_write_jobRealFile(777,'F_tau_lastInc'//trim(rankStr),size(F_tau_lastInc)) ! writing F_tau_lastInc field to file
write (777,rec=1) F_tau_lastInc; close (777)
fileUnit = IO_open_jobFile_binary('F'//trim(rankStr),'w')
write(fileUnit) F; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr),'w')
write(fileUnit) F_lastInc; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_tau'//trim(rankStr),'w')
write(fileUnit) F_tau; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_tau_lastInc'//trim(rankStr),'w')
write(fileUnit) F_tau_lastInc; close (fileUnit)
endif
call CPFEM_age() ! age state and kinematics
@ -618,6 +618,7 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,deformati
F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess)
F_aim_lastInc = F_aim
!--------------------------------------------------------------------------------------------------
! calculate rate for aim
if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F

108
src/spectral_utilities.f90
View File

@ -362,62 +362,64 @@ end subroutine utilities_init
!> Also writes out the current reference stiffness for restart.
!--------------------------------------------------------------------------------------------------
subroutine utilities_updateGamma(C,saveReference)
use IO, only: &
IO_write_jobRealFile
use numerics, only: &
memory_efficient, &
worldrank
use mesh, only: &
grid3Offset, &
grid3,&
grid
use math, only: &
math_det33, &
math_invert
implicit none
real(pReal), intent(in), dimension(3,3,3,3) :: C !< input stiffness to store as reference stiffness
logical , intent(in) :: saveReference !< save reference stiffness to file for restart
complex(pReal), dimension(3,3) :: temp33_complex, xiDyad_cmplx
real(pReal), dimension(6,6) :: matA, matInvA
integer(pInt) :: &
i, j, k, &
l, m, n, o
logical :: err
C_ref = C
if (saveReference) then
if (worldrank == 0_pInt) then
write(6,'(/,a)') ' writing reference stiffness to file'
flush(6)
call IO_write_jobRealFile(777,'C_ref',size(C_ref))
write (777,rec=1) C_ref; close(777)
endif
endif
if(.not. memory_efficient) then
gamma_hat = cmplx(0.0_pReal,0.0_pReal,pReal) ! for the singular point and any non invertible A
do k = grid3Offset+1_pInt, grid3Offset+grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid1Red
if (any([i,j,k] /= 1_pInt)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) &
xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-grid3Offset))*xi1st(m,i,j,k-grid3Offset)
forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) &
temp33_complex(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx)
matA(1:3,1:3) = real(temp33_complex); matA(4:6,4:6) = real(temp33_complex)
matA(1:3,4:6) = aimag(temp33_complex); matA(4:6,1:3) = -aimag(temp33_complex)
if (abs(math_det33(matA(1:3,1:3))) > 1e-16) then
call math_invert(6_pInt, matA, matInvA, err)
temp33_complex = cmplx(matInvA(1:3,1:3),matInvA(1:3,4:6),pReal)
forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, o=1_pInt:3_pInt) &
gamma_hat(l,m,n,o,i,j,k-grid3Offset) = temp33_complex(l,n)* &
use IO, only: &
IO_open_jobFile_binary
use numerics, only: &
memory_efficient, &
worldrank
use mesh, only: &
grid3Offset, &
grid3,&
grid
use math, only: &
math_det33, &
math_invert2
implicit none
real(pReal), intent(in), dimension(3,3,3,3) :: C !< input stiffness to store as reference stiffness
logical , intent(in) :: saveReference !< save reference stiffness to file for restart
complex(pReal), dimension(3,3) :: temp33_complex, xiDyad_cmplx
real(pReal), dimension(6,6) :: A, A_inv
integer :: &
i, j, k, &
l, m, n, o, &
fileUnit
logical :: err
C_ref = C
if (saveReference) then
if (worldrank == 0_pInt) then
write(6,'(/,a)') ' writing reference stiffness to file'
flush(6)
fileUnit = IO_open_jobFile_binary('C_ref','w')
write(fileUnit) C_ref; close(fileUnit)
endif
endif
if(.not. memory_efficient) then
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
if (any([i,j,k] /= 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) &
xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-grid3Offset))*xi1st(m,i,j,k-grid3Offset)
forall(l = 1:3, m = 1:3) &
temp33_complex(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx)
A(1:3,1:3) = real(temp33_complex); A(4:6,4:6) = real(temp33_complex)
A(1:3,4:6) = aimag(temp33_complex); A(4:6,1:3) = -aimag(temp33_complex)
if (abs(math_det33(A(1:3,1:3))) > 1e-16) then
call math_invert2(A_inv, err, A)
temp33_complex = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal)
forall(l=1:3, m=1:3, n=1:3, o=1:3) &
gamma_hat(l,m,n,o,i,j,k-grid3Offset) = temp33_complex(l,n)* &
conjg(-xi1st(o,i,j,k-grid3Offset))*xi1st(m,i,j,k-grid3Offset)
endif
endif
enddo; enddo; enddo
endif
endif
endif
enddo; enddo; enddo
endif
end subroutine utilities_updateGamma
!--------------------------------------------------------------------------------------------------
!> @brief forward FFT of data in field_real to field_fourier
!> @details Does an unweighted filtered FFT transform from real to complex