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

reordered (and partly redistributed) error message identifiers, deleted those which are not in use anymore. 

all constitutive as well numerics now raises an error, if an unknown keyword is found in the respective config file
This commit is contained in:
Christoph Kords 2012-02-13 17:41:27 +00:00
parent dd51e1da81
commit c786336af3
19 changed files with 796 additions and 763 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

193
code/CPFEM.f90
View File

@ -149,54 +149,54 @@ subroutine CPFEM_init()
write(6,'(a)') '<< CPFEM >> Restored state variables of last converged step from binary files'
!$OMP END CRITICAL (write2out)
endif
if (IO_read_jobBinaryFile(777,'recordedPhase',FEmodelGeometry,size(material_phase))) then
read (777,rec=1) material_phase
close (777)
endif
if (IO_read_jobBinaryFile(777,'convergedF',FEmodelGeometry,size(crystallite_F0))) then
read (777,rec=1) crystallite_F0
close (777)
endif
if (IO_read_jobBinaryFile(777,'convergedFp',FEmodelGeometry,size(crystallite_Fp0))) then
read (777,rec=1) crystallite_Fp0
close (777)
endif
if (IO_read_jobBinaryFile(777,'convergedLp',FEmodelGeometry,size(crystallite_Lp0))) then
read (777,rec=1) crystallite_Lp0
close (777)
endif
if (IO_read_jobBinaryFile(777,'convergeddPdF',FEmodelGeometry,size(crystallite_dPdF0))) then
read (777,rec=1) crystallite_dPdF0
close (777)
endif
if (IO_read_jobBinaryFile(777,'convergedTstar',FEmodelGeometry,size(crystallite_Tstar0_v))) then
read (777,rec=1) crystallite_Tstar0_v
close (777)
endif
if (IO_read_jobBinaryFile(777,'convergedStateConst',FEmodelGeometry)) then
m = 0_pInt
do i = 1,homogenization_maxNgrains; do j = 1,mesh_maxNips; do k = 1,mesh_NcpElems
do l = 1,size(constitutive_state0(i,j,k)%p)
m = m+1_pInt
read(777,rec=m) constitutive_state0(i,j,k)%p(l)
enddo
enddo; enddo; enddo
close (777)
endif
if (IO_read_jobBinaryFile(777,'convergedStateHomog',FEmodelGeometry)) then
m = 0_pInt
do k = 1,mesh_NcpElems; do j = 1,mesh_maxNips
do l = 1,homogenization_sizeState(j,k)
m = m+1_pInt
read(777,rec=m) homogenization_state0(j,k)%p(l)
enddo
enddo; enddo
close (777)
endif
if (IO_read_jobBinaryFile(777,'convergeddcsdE',FEmodelGeometry,size(CPFEM_dcsdE))) then
read (777,rec=1) CPFEM_dcsdE
close (777)
endif
call IO_read_jobBinaryFile(777,'recordedPhase',FEmodelGeometry,size(material_phase))
read (777,rec=1) material_phase
close (777)
call IO_read_jobBinaryFile(777,'convergedF',FEmodelGeometry,size(crystallite_F0))
read (777,rec=1) crystallite_F0
close (777)
call IO_read_jobBinaryFile(777,'convergedFp',FEmodelGeometry,size(crystallite_Fp0))
read (777,rec=1) crystallite_Fp0
close (777)
call IO_read_jobBinaryFile(777,'convergedLp',FEmodelGeometry,size(crystallite_Lp0))
read (777,rec=1) crystallite_Lp0
close (777)
call IO_read_jobBinaryFile(777,'convergeddPdF',FEmodelGeometry,size(crystallite_dPdF0))
read (777,rec=1) crystallite_dPdF0
close (777)
call IO_read_jobBinaryFile(777,'convergedTstar',FEmodelGeometry,size(crystallite_Tstar0_v))
read (777,rec=1) crystallite_Tstar0_v
close (777)
call IO_read_jobBinaryFile(777,'convergedStateConst',FEmodelGeometry)
m = 0_pInt
do i = 1,homogenization_maxNgrains; do j = 1,mesh_maxNips; do k = 1,mesh_NcpElems
do l = 1,size(constitutive_state0(i,j,k)%p)
m = m+1_pInt
read(777,rec=m) constitutive_state0(i,j,k)%p(l)
enddo
enddo; enddo; enddo
close (777)
call IO_read_jobBinaryFile(777,'convergedStateHomog',FEmodelGeometry)
m = 0_pInt
do k = 1,mesh_NcpElems; do j = 1,mesh_maxNips
do l = 1,homogenization_sizeState(j,k)
m = m+1_pInt
read(777,rec=m) homogenization_state0(j,k)%p(l)
enddo
enddo; enddo
close (777)
call IO_read_jobBinaryFile(777,'convergeddcsdE',FEmodelGeometry,size(CPFEM_dcsdE))
read (777,rec=1) CPFEM_dcsdE
close (777)
restartRead = .false.
endif
! *** end of restoring
@ -427,54 +427,55 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
write(6,'(a)') '<< CPFEM >> Writing state variables of last converged step to binary files'
!$OMP END CRITICAL (write2out)
endif
if (IO_write_jobBinaryFile(777,'recordedPhase',size(material_phase))) then
write (777,rec=1) material_phase
close (777)
endif
if (IO_write_jobBinaryFile(777,'convergedF',size(crystallite_F0))) then
write (777,rec=1) crystallite_F0
close (777)
endif
if (IO_write_jobBinaryFile(777,'convergedFp',size(crystallite_Fp0))) then
write (777,rec=1) crystallite_Fp0
close (777)
endif
if (IO_write_jobBinaryFile(777,'convergedLp',size(crystallite_Lp0))) then
write (777,rec=1) crystallite_Lp0
close (777)
endif
if (IO_write_jobBinaryFile(777,'convergeddPdF',size(crystallite_dPdF0))) then
write (777,rec=1) crystallite_dPdF0
close (777)
endif
if (IO_write_jobBinaryFile(777,'convergedTstar',size(crystallite_Tstar0_v))) then
write (777,rec=1) crystallite_Tstar0_v
close (777)
endif
if (IO_write_jobBinaryFile(777,'convergedStateConst')) then
m = 0_pInt
do i = 1,homogenization_maxNgrains; do j = 1,mesh_maxNips; do k = 1,mesh_NcpElems
do l = 1,size(constitutive_state0(i,j,k)%p)
m = m+1_pInt
write(777,rec=m) constitutive_state0(i,j,k)%p(l)
enddo
enddo; enddo; enddo
close (777)
endif
if (IO_write_jobBinaryFile(777,'convergedStateHomog')) then
m = 0_pInt
do k = 1,mesh_NcpElems; do j = 1,mesh_maxNips
do l = 1,homogenization_sizeState(j,k)
m = m+1_pInt
write(777,rec=m) homogenization_state0(j,k)%p(l)
enddo
enddo; enddo
close (777)
endif
if (IO_write_jobBinaryFile(777,'convergeddcsdE',size(CPFEM_dcsdE))) then
write (777,rec=1) CPFEM_dcsdE
close (777)
endif
call IO_write_jobBinaryFile(777,'recordedPhase',size(material_phase))
write (777,rec=1) material_phase
close (777)
call IO_write_jobBinaryFile(777,'convergedF',size(crystallite_F0))
write (777,rec=1) crystallite_F0
close (777)
call IO_write_jobBinaryFile(777,'convergedFp',size(crystallite_Fp0))
write (777,rec=1) crystallite_Fp0
close (777)
call IO_write_jobBinaryFile(777,'convergedLp',size(crystallite_Lp0))
write (777,rec=1) crystallite_Lp0
close (777)
call IO_write_jobBinaryFile(777,'convergeddPdF',size(crystallite_dPdF0))
write (777,rec=1) crystallite_dPdF0
close (777)
call IO_write_jobBinaryFile(777,'convergedTstar',size(crystallite_Tstar0_v))
write (777,rec=1) crystallite_Tstar0_v
close (777)
call IO_write_jobBinaryFile(777,'convergedStateConst')
m = 0_pInt
do i = 1,homogenization_maxNgrains; do j = 1,mesh_maxNips; do k = 1,mesh_NcpElems
do l = 1,size(constitutive_state0(i,j,k)%p)
m = m+1_pInt
write(777,rec=m) constitutive_state0(i,j,k)%p(l)
enddo
enddo; enddo; enddo
close (777)
call IO_write_jobBinaryFile(777,'convergedStateHomog')
m = 0_pInt
do k = 1,mesh_NcpElems; do j = 1,mesh_maxNips
do l = 1,homogenization_sizeState(j,k)
m = m+1_pInt
write(777,rec=m) homogenization_state0(j,k)%p(l)
enddo
enddo; enddo
close (777)
call IO_write_jobBinaryFile(777,'convergeddcsdE',size(CPFEM_dcsdE))
write (777,rec=1) CPFEM_dcsdE
close (777)
endif
! * end of dumping
endif

2
code/DAMASK_marc.f90
View File

@ -400,7 +400,7 @@ subroutine plotv(&
real(pReal) v, t(*)
integer(pInt) m, nn, layer, ndi, nshear, jpltcd
if (jpltcd > materialpoint_sizeResults) call IO_error(667,jpltcd) ! complain about out of bounds error
if (jpltcd > materialpoint_sizeResults) call IO_error(700_pInt,jpltcd) ! complain about out of bounds error
v = materialpoint_results(jpltcd,nn,mesh_FEasCP('elem', m))
return

57
code/DAMASK_spectral.f90
View File

@ -210,7 +210,7 @@ program DAMASK_spectral
!--------------------------------------------------------------------------------------------------
! reading the load case file and allocate data structure containing load cases
path = getLoadcaseName()
if (.not. IO_open_file(myUnit,path)) call IO_error(error_ID = 30_pInt,ext_msg = trim(path))
call IO_open_file(myUnit,path)
rewind(myUnit)
do
read(myUnit,'(a1024)',END = 100) line
@ -232,7 +232,7 @@ program DAMASK_spectral
100 N_Loadcases = N_n
if ((N_l + N_Fdot /= N_n) .or. (N_n /= N_t)) & ! sanity check
call IO_error(error_ID=37_pInt,ext_msg = trim(path)) ! error message for incomplete loadcase
call IO_error(error_ID=837_pInt,ext_msg = trim(path)) ! error message for incomplete loadcase
allocate (bc(N_Loadcases))
@ -310,14 +310,13 @@ program DAMASK_spectral
!-------------------------------------------------------------------------------------------------- ToDo: if temperature at CPFEM is treated properly, move this up immediately after interface init
! initialization of all related DAMASK modules (e.g. mesh.f90 reads in geometry)
call CPFEM_initAll(bc(1)%temperature,1_pInt,1_pInt)
if (update_gamma .and. .not. memory_efficient) call IO_error(error_ID = 47_pInt)
if (update_gamma .and. .not. memory_efficient) call IO_error(error_ID = 847_pInt)
!--------------------------------------------------------------------------------------------------
! read header of geom file to get size information. complete geom file is intepretated by mesh.f90
path = getModelName()
if (.not. IO_open_file(myUnit,trim(path)//InputFileExtension))&
call IO_error(error_ID=101_pInt,ext_msg = trim(path)//InputFileExtension)
call IO_open_file(myUnit,trim(path)//InputFileExtension)
rewind(myUnit)
read(myUnit,'(a1024)') line
positions = IO_stringPos(line,2_pInt)
@ -325,7 +324,7 @@ program DAMASK_spectral
if (keyword(1:4) == 'head') then
headerLength = IO_intValue(line,positions,1_pInt) + 1_pInt
else
call IO_error(error_ID=42_pInt)
call IO_error(error_ID=842_pInt)
endif
rewind(myUnit)
@ -367,12 +366,12 @@ program DAMASK_spectral
!--------------------------------------------------------------------------------------------------
! sanity checks of geometry parameters
if (.not.(gotDimension .and. gotHomogenization .and. gotResolution))&
call IO_error(error_ID = 45_pInt)
if (any(geomdim<=0.0_pReal)) call IO_error(error_ID = 102_pInt)
call IO_error(error_ID = 845_pInt)
if (any(geomdim<=0.0_pReal)) call IO_error(error_ID = 802_pInt)
if(mod(res(1),2_pInt)/=0_pInt .or.&
mod(res(2),2_pInt)/=0_pInt .or.&
(mod(res(3),2_pInt)/=0_pInt .and. res(3)/= 1_pInt))&
call IO_error(error_ID = 103_pInt)
call IO_error(error_ID = 803_pInt)
!--------------------------------------------------------------------------------------------------
! variables derived from resolution
@ -413,7 +412,7 @@ program DAMASK_spectral
if (bc(loadcase)%velGradApplied) then
do j = 1_pInt, 3_pInt
if (any(bc(loadcase)%maskDeformation(j,1:3) .eqv. .true.) .and. &
any(bc(loadcase)%maskDeformation(j,1:3) .eqv. .false.)) errorID = 32_pInt ! each row should be either fully or not at all defined
any(bc(loadcase)%maskDeformation(j,1:3) .eqv. .false.)) errorID = 832_pInt ! each row should be either fully or not at all defined
enddo
print '(a)','velocity gradient:'
else
@ -433,17 +432,17 @@ program DAMASK_spectral
print '(a,i5)','output frequency: ',bc(loadcase)%outputfrequency
print '(a,i5)','restart frequency: ',bc(loadcase)%restartfrequency
if (any(bc(loadcase)%maskStress .eqv. bc(loadcase)%maskDeformation)) errorID = 31 ! exclusive or masking only
if (any(bc(loadcase)%maskStress .eqv. bc(loadcase)%maskDeformation)) errorID = 831_pInt ! exclusive or masking only
if (any(bc(loadcase)%maskStress .and. transpose(bc(loadcase)%maskStress) .and. &
reshape([ .false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false.],[ 3,3]))) &
errorID = 38_pInt ! no rotation is allowed by stress BC
errorID = 838_pInt ! no rotation is allowed by stress BC
if (any(abs(math_mul33x33(bc(loadcase)%rotation,math_transpose33(bc(loadcase)%rotation))&
-math_I3) > reshape(spread(rotation_tol,1,9),[ 3,3]))&
.or. abs(math_det33(bc(loadcase)%rotation)) > 1.0_pReal + rotation_tol)&
errorID = 46_pInt ! given rotation matrix contains strain
if (bc(loadcase)%time < 0.0_pReal) errorID = 34_pInt ! negative time increment
if (bc(loadcase)%incs < 1_pInt) errorID = 35_pInt ! non-positive incs count
if (bc(loadcase)%outputfrequency < 1_pInt) errorID = 36_pInt ! non-positive result frequency
errorID = 846_pInt ! given rotation matrix contains strain
if (bc(loadcase)%time < 0.0_pReal) errorID = 834_pInt ! negative time increment
if (bc(loadcase)%incs < 1_pInt) errorID = 835_pInt ! non-positive incs count
if (bc(loadcase)%outputfrequency < 1_pInt) errorID = 836_pInt ! non-positive result frequency
if (errorID > 0_pInt) call IO_error(error_ID = errorID, ext_msg = loadcase_string)
enddo
@ -483,7 +482,7 @@ program DAMASK_spectral
c0_reference = c_current * wgt ! linear reference material stiffness
c0_66 = math_Mandel3333to66(c0_reference)
call math_invert(6_pInt, c0_66, s0_66, i, errmatinv) ! invert in mandel notation
if(errmatinv) call IO_error(error_ID=800_pInt)
if(errmatinv) call IO_error(error_ID=400_pInt)
s0_reference = math_Mandel66to3333(s0_66)
!--------------------------------------------------------------------------------------------------
@ -491,11 +490,10 @@ program DAMASK_spectral
if (restartInc > 1_pInt) then ! using old values from file
if (debugRestart) print '(a,i6,a)' , 'Reading values of increment ',&
restartInc - 1_pInt,' from file'
if (IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',&
trim(getSolverJobName()),size(defgrad))) then
read (777,rec=1) defgrad
close (777)
endif
call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',&
trim(getSolverJobName()),size(defgrad))
read (777,rec=1) defgrad
close (777)
defgradold = defgrad
defgradAim = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
@ -544,11 +542,11 @@ program DAMASK_spectral
!--------------------------------------------------------------------------------------------------
! general initialization of fftw (see manual on fftw.org for more details)
if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=108_pInt) ! check for correct precision in C
if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=808_pInt) ! check for correct precision in C
#ifdef _OPENMP
if(DAMASK_NumThreadsInt > 0_pInt) then
ierr = fftw_init_threads()
if (ierr == 0_pInt) call IO_error(error_ID = 109_pInt)
if (ierr == 0_pInt) call IO_error(error_ID = 809_pInt)
call fftw_plan_with_nthreads(DAMASK_NumThreadsInt)
endif
#endif
@ -723,7 +721,7 @@ program DAMASK_spectral
c_reduced(k,j) = c_prev99(n,m)
endif; enddo; endif; enddo
call math_invert(size_reduced, c_reduced, s_reduced, i, errmatinv) ! invert reduced stiffness
if(errmatinv) call IO_error(error_ID=800_pInt)
if(errmatinv) call IO_error(error_ID=400_pInt)
s_prev99 = 0.0_pReal ! build full compliance
k = 0_pInt
do n = 1_pInt,9_pInt
@ -1097,10 +1095,9 @@ program DAMASK_spectral
mod(inc - 1_pInt,bc(loadcase)%restartFrequency) == 0_pInt) then ! at frequency of writing restart information set restart parameter for FEsolving (first call to CPFEM_general will write ToDo: true?)
restartWrite = .true.
print '(A)', 'writing converged results for restart'
if(IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(defgrad))) then ! writing deformation gradient field to file
write (777,rec=1) defgrad
close (777)
endif
call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(defgrad)) ! writing deformation gradient field to file
write (777,rec=1) defgrad
close (777)
restartInc=totalIncsCounter
endif
@ -1111,7 +1108,7 @@ program DAMASK_spectral
!c0_99 = math_Plain3333to99(c0_reference)
! call math_invert(9_pInt, s0_99, c0_99, i, errmatinv) ! invert reduced stiffness
! if(errmatinv) call IO_error(error_ID=800_pInt)
! if(errmatinv) call IO_error(error_ID=400_pInt)
! print*, (c0_reference - math_Plain99to3333(c0_99))/c0_reference
! pause
endif

103
code/FEsolving.f90
View File

@ -60,63 +60,61 @@ integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions
character(len=1024) line, commandLine
FEmodelGeometry = getModelName()
if (IO_open_inputFile(fileunit,FEmodelGeometry)) then
if (trim(FEsolver) == 'Spectral') then
call get_command(commandLine) ! may contain uppercase
do i=1,len(commandLine)
if(64 < iachar(commandLine(i:i)) .and. iachar(commandLine(i:i)) < 91) &
commandLine(i:i) = achar(iachar(commandLine(i:i))+32) ! make lowercase
enddo
if (index(commandLine,'-r ',.true.)>0) & ! look for -r
start = index(commandLine,'-r ',.true.) + 3_pInt ! set to position after trailing space
if (index(commandLine,'--restart ',.true.)>0) & ! look for --restart
start = index(commandLine,'--restart ',.true.) + 10_pInt ! set to position after trailing space
call IO_open_inputFile(fileunit,FEmodelGeometry)
if (trim(FEsolver) == 'Spectral') then
call get_command(commandLine) ! may contain uppercase
do i=1,len(commandLine)
if(64 < iachar(commandLine(i:i)) .and. iachar(commandLine(i:i)) < 91) &
commandLine(i:i) = achar(iachar(commandLine(i:i))+32) ! make lowercase
enddo
if (index(commandLine,'-r ',.true.)>0) & ! look for -r
start = index(commandLine,'-r ',.true.) + 3_pInt ! set to position after trailing space
if (index(commandLine,'--restart ',.true.)>0) & ! look for --restart
start = index(commandLine,'--restart ',.true.) + 10_pInt ! set to position after trailing space
if(start /= 0_pInt) then ! found something
length = verify(commandLine(start:len(commandLine)),'0123456789',.false.) ! where is first non number after argument?
read(commandLine(start:start+length),'(I12)') restartInc ! read argument
restartRead = restartInc > 0_pInt
if(restartInc <= 0_pInt) then
call IO_warning(warning_ID=34_pInt)
restartInc = 1_pInt
endif
if(start /= 0_pInt) then ! found something
length = verify(commandLine(start:len(commandLine)),'0123456789',.false.) ! where is first non number after argument?
read(commandLine(start:start+length),'(I12)') restartInc ! read argument
restartRead = restartInc > 0_pInt
if(restartInc <= 0_pInt) then
call IO_warning(warning_ID=34_pInt)
restartInc = 1_pInt
endif
else
rewind(fileunit)
do
read (fileunit,'(a1024)',END=100) line
positions = IO_stringPos(line,maxNchunks)
tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key
select case(tag)
case ('solver')
read (fileunit,'(a1024)',END=100) line ! next line
positions = IO_stringPos(line,maxNchunks)
symmetricSolver = (IO_intValue(line,positions,2_pInt) /= 1_pInt)
case ('restart')
read (fileunit,'(a1024)',END=100) line ! next line
positions = IO_stringPos(line,maxNchunks)
restartWrite = iand(IO_intValue(line,positions,1_pInt),1_pInt) > 0_pInt
restartRead = iand(IO_intValue(line,positions,1_pInt),2_pInt) > 0_pInt
case ('*restart')
do i=2,positions(1)
restartWrite = (IO_lc(IO_StringValue(line,positions,i)) == 'write') .or. restartWrite
restartRead = (IO_lc(IO_StringValue(line,positions,i)) == 'read') .or. restartRead
enddo
if(restartWrite) then
do i=2,positions(1)
restartWrite = (IO_lc(IO_StringValue(line,positions,i)) /= 'frequency=0') .and. restartWrite
enddo
endif
end select
enddo
endif
else
call IO_error(101_pInt, ext_msg=FEmodelGeometry) ! cannot open input file
rewind(fileunit)
do
read (fileunit,'(a1024)',END=100) line
positions = IO_stringPos(line,maxNchunks)
tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key
select case(tag)
case ('solver')
read (fileunit,'(a1024)',END=100) line ! next line
positions = IO_stringPos(line,maxNchunks)
symmetricSolver = (IO_intValue(line,positions,2_pInt) /= 1_pInt)
case ('restart')
read (fileunit,'(a1024)',END=100) line ! next line
positions = IO_stringPos(line,maxNchunks)
restartWrite = iand(IO_intValue(line,positions,1_pInt),1_pInt) > 0_pInt
restartRead = iand(IO_intValue(line,positions,1_pInt),2_pInt) > 0_pInt
case ('*restart')
do i=2,positions(1)
restartWrite = (IO_lc(IO_StringValue(line,positions,i)) == 'write') .or. restartWrite
restartRead = (IO_lc(IO_StringValue(line,positions,i)) == 'read') .or. restartRead
enddo
if(restartWrite) then
do i=2,positions(1)
restartWrite = (IO_lc(IO_StringValue(line,positions,i)) /= 'frequency=0') .and. restartWrite
enddo
endif
end select
enddo
endif
100 close(fileunit)
if (restartRead) then
if(FEsolver == 'Marc' .and. IO_open_logFile(fileunit)) then
if(FEsolver == 'Marc') then
call IO_open_logFile(fileunit)
rewind(fileunit)
do
read (fileunit,'(a1024)',END=200) line
@ -127,7 +125,8 @@ integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions
IO_lc(IO_stringValue(line,positions,4_pInt)) == 'id' ) &
FEmodelGeometry = IO_StringValue(line,positions,6_pInt)
enddo
elseif (FEsolver == 'Abaqus' .and. IO_open_inputFile(fileunit,FEmodelGeometry)) then
elseif (FEsolver == 'Abaqus') then
call IO_open_inputFile(fileunit,FEmodelGeometry)
rewind(fileunit)
do
read (fileunit,'(a1024)',END=200) line
@ -138,10 +137,6 @@ integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions
FEmodelGeometry = IO_StringValue(line,positions,1_pInt)
endif
enddo
elseif (FEsolver == 'Spectral') then
!do nothing
else
call IO_error(106_pInt) ! cannot open file for old job info
endif
endif

518
code/IO.f90
View File

@ -25,9 +25,9 @@
CONTAINS
!---------------------------
! function IO_abaqus_assembleInputFile
! function IO_open_file(unit,relPath)
! function IO_open_inputFile(unit, model)
! function IO_open_logFile(unit)
! subroutine IO_open_file(unit,relPath)
! subroutine IO_open_inputFile(unit, model)
! subroutine IO_open_logFile(unit)
! function IO_hybridIA(Nast,ODFfileName)
! private function hybridIA_reps(dV_V,steps,C)
! function IO_stringPos(line,maxN)
@ -153,21 +153,49 @@ end function
! open existing file to given unit
! path to file is relative to working directory
!********************************************************************
logical function IO_open_file(unit,relPath)
logical function IO_open_file_stat(unit,relPath)
use prec, only: pInt
use DAMASK_interface
implicit none
character(len=*) relPath
integer(pInt) unit
integer(pInt), intent(in) :: unit
character(len=*), intent(in) :: relPath
integer stat
character path
IO_open_file_stat = .false.
path = trim(getSolverWorkingDirectoryName())//relPath
open(unit,status='old',iostat=stat,file=path)
if (stat == 0) then
IO_open_file_stat = .true.
endif
endfunction
IO_open_file = .false.
!********************************************************************
! open existing file to given unit
! path to file is relative to working directory
!********************************************************************
subroutine IO_open_file(unit,relPath)
use prec, only: pInt
use DAMASK_interface
open(unit,status='old',err=100,file=trim(getSolverWorkingDirectoryName())//relPath)
IO_open_file = .true.
implicit none
integer(pInt), intent(in) :: unit
character(len=*), intent(in) :: relPath
character path
integer(pInt) stat
100 endfunction
path = trim(getSolverWorkingDirectoryName())//relPath
open(unit,status='old',iostat=stat,file=path)
if (stat /= 0) then
call IO_error(100_pInt,ext_msg=path)
endif
endsubroutine
!********************************************************************
@ -176,59 +204,74 @@ end function
! : changed the function to open *.inp_assembly, which is basically
! the input file without comment lines and possibly assembled includes
!********************************************************************
logical function IO_open_inputFile(unit,model)
subroutine IO_open_inputFile(unit,model)
use prec, only: pReal, pInt
use DAMASK_interface
implicit none
integer(pInt), intent(in) :: unit
character(len=*) model
character(len=*), intent(in) :: model
integer(pInt) stat
character path
IO_open_inputFile = .false.
if (FEsolver == 'Abaqus') then
open(unit+1,status='old',err=100,&
file=trim(getSolverWorkingDirectoryName())//&
trim(model)//InputFileExtension)
open(unit,err=100,file=trim(getSolverWorkingDirectoryName())//&
trim(model)//InputFileExtension//'_assembly')
IO_open_inputFile = IO_abaqus_assembleInputFile(unit,unit+1_pInt) ! strip comments and concatenate any "include"s
path = trim(getSolverWorkingDirectoryName())//trim(model)//InputFileExtension
open(unit+1,status='old',iostat=stat,file=path)
if (stat /= 0) then
call IO_error(100_pInt,ext_msg=path)
endif
path = trim(getSolverWorkingDirectoryName())//trim(model)//InputFileExtension//'_assembly'
open(unit,iostat=stat,file=path)
if (stat /= 0) then
call IO_error(100_pInt,ext_msg=path)
endif
if (IO_abaqus_assembleInputFile(unit,unit+1_pInt)) then ! strip comments and concatenate any "include"s
call IO_error(103_pInt)
endif
close(unit+1_pInt)
else
open(unit,status='old',err=100,file=trim(getSolverWorkingDirectoryName())//&
trim(model)//InputFileExtension)
IO_open_inputFile = .true.
path = trim(getSolverWorkingDirectoryName())//trim(model)//InputFileExtension
open(unit,status='old',iostat=stat,file=path)
if (stat /= 0) then
call IO_error(100_pInt,ext_msg=path)
endif
endif
100 endfunction
endsubroutine
!********************************************************************
! open FEM logfile to given unit
!********************************************************************
logical function IO_open_logFile(unit)
subroutine IO_open_logFile(unit)
use prec, only: pReal, pInt
use DAMASK_interface
implicit none
character path
integer(pInt) stat
integer(pInt), intent(in) :: unit
IO_open_logFile = .false.
open(unit,status='old',err=100,file=trim(getSolverWorkingDirectoryName())//&
trim(getSolverJobName())//LogFileExtension)
IO_open_logFile = .true.
path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//LogFileExtension
open(unit,status='old',iostat=stat,file=path)
if (stat /= 0) then
call IO_error(100_pInt,ext_msg=path)
endif
100 endfunction
endsubroutine
!********************************************************************
! open (write) file related to current job
! but with different extension to given unit
!********************************************************************
logical function IO_open_jobFile(unit,newExt)
logical function IO_open_jobFile_stat(unit,newExt)
use prec, only: pReal, pInt
use DAMASK_interface
@ -236,21 +279,24 @@ end function
integer(pInt), intent(in) :: unit
character(*), intent(in) :: newExt
character path
integer stat
IO_open_jobFile = .false.
IO_open_jobFile_stat = .false.
path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//newExt
open(unit,status='old',iostat=stat,file=path)
if (stat == 0) then
IO_open_jobFile_stat = .true.
endif
open(unit,status='old',err=100,file=trim(getSolverWorkingDirectoryName())//&
trim(getSolverJobName())//'.'//newExt)
IO_open_jobFile = .true.
100 endfunction
endfunction
!********************************************************************
! open (write) file related to current job
! but with different extension to given unit
!********************************************************************
logical function IO_write_jobFile(unit,newExt)
subroutine IO_open_jobFile(unit,newExt)
use prec, only: pReal, pInt
use DAMASK_interface
@ -258,21 +304,47 @@ end function
integer(pInt), intent(in) :: unit
character(*), intent(in) :: newExt
character path
integer(pInt) stat
IO_write_jobFile = .false.
path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//newExt
open(unit,status='old',iostat=stat,file=path)
if (stat /= 0) then
call IO_error(100_pInt,ext_msg=path)
endif
open(unit,status='replace',err=100,file=trim(getSolverWorkingDirectoryName())//&
trim(getSolverJobName())//'.'//newExt)
IO_write_jobFile = .true.
endsubroutine
!********************************************************************
! open (write) file related to current job
! but with different extension to given unit
!********************************************************************
subroutine IO_write_jobFile(unit,newExt)
use prec, only: pReal, pInt
use DAMASK_interface
implicit none
integer(pInt), intent(in) :: unit
character(*), intent(in) :: newExt
character path
integer(pInt) stat
path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//newExt
open(unit,status='replace',iostat=stat,file=path)
if (stat /= 0) then
call IO_error(100_pInt,ext_msg=path)
endif
100 endfunction
endsubroutine
!********************************************************************
! open (write) binary file related to current job
! but with different extension to given unit
!********************************************************************
logical function IO_write_jobBinaryFile(unit,newExt,recMultiplier)
subroutine IO_write_jobBinaryFile(unit,newExt,recMultiplier)
use prec, only: pReal, pInt
use DAMASK_interface
@ -281,27 +353,27 @@ end function
integer(pInt), intent(in) :: unit
integer(pInt), intent(in), optional :: recMultiplier
character(*), intent(in) :: newExt
character path
integer(pInt) stat
IO_write_jobBinaryFile = .false.
path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//newExt
if (present(recMultiplier)) then
open(unit,status='replace',form='unformatted',access='direct',recl=pReal*recMultiplier, &
err=100,file=trim(getSolverWorkingDirectoryName())//&
trim(getSolverJobName())//'.'//newExt)
else
open(unit,status='replace',form='unformatted',access='direct',recl=pReal, &
err=100,file=trim(getSolverWorkingDirectoryName())//&
trim(getSolverJobName())//'.'//newExt)
open(unit,status='replace',form='unformatted',access='direct',recl=pReal*recMultiplier,iostat=stat,file=path)
else
open(unit,status='replace',form='unformatted',access='direct',recl=pReal,iostat=stat,file=path)
endif
if (stat /= 0) then
call IO_error(100_pInt,ext_msg=path)
endif
IO_write_jobBinaryFile = .true.
100 endfunction
endsubroutine
!********************************************************************
! open (read) binary file related to restored job
! and with different extension to given unit
!********************************************************************
logical function IO_read_jobBinaryFile(unit,newExt,jobName,recMultiplier)
subroutine IO_read_jobBinaryFile(unit,newExt,jobName,recMultiplier)
use prec, only: pReal, pInt
use DAMASK_interface
@ -310,20 +382,20 @@ end function
integer(pInt), intent(in) :: unit
integer(pInt), intent(in), optional :: recMultiplier
character(*), intent(in) :: newExt, jobName
character path
integer(pInt) stat
IO_read_jobBinaryFile = .false.
path = trim(getSolverWorkingDirectoryName())//trim(jobName)//'.'//newExt
if (present(recMultiplier)) then
open(unit,status='old',form='unformatted',access='direct',recl=pReal*recMultiplier, &
err=100,file=trim(getSolverWorkingDirectoryName())//&
trim(jobName)//'.'//newExt)
else
open(unit,status='old',form='unformatted',access='direct',recl=pReal, &
err=100,file=trim(getSolverWorkingDirectoryName())//&
trim(jobName)//'.'//newExt)
open(unit,status='old',form='unformatted',access='direct',recl=pReal*recMultiplier,iostat=stat,file=path)
else
open(unit,status='old',form='unformatted',access='direct',recl=pReal,iostat=stat,file=path)
endif
if (stat /= 0) then
call IO_error(100_pInt,ext_msg=path)
endif
IO_read_jobBinaryFile = .true.
100 endfunction
endsubroutine
!********************************************************************
@ -374,7 +446,7 @@ end function
real(pReal), dimension(:,:,:), allocatable :: dV_V
real(pReal), dimension(3,Nast) :: IO_hybridIA
if (.not. IO_open_file(999_pInt,ODFfileName)) goto 100
call IO_open_file(999_pInt,ODFfileName)
!--- parse header of ODF file ---
!--- limits in phi1, Phi, phi2 ---
@ -1135,225 +1207,172 @@ endfunction
character(len=1024) msg
select case (error_ID)
case (30_pInt)
msg = 'could not open spectral loadcase'
case (31_pInt)
msg = 'mask consistency violated in spectral loadcase'
case (32_pInt)
msg = 'ill-defined L (each line should be either fully or not at all defined) in spectral loadcase'
case (34_pInt)
msg = 'negative time increment in spectral loadcase'
case (35_pInt)
msg = 'non-positive increments in spectral loadcase'
case (36_pInt)
msg = 'non-positive result frequency in spectral loadcase'
case (37_pInt)
msg = 'incomplete loadcase'
case (38_pInt)
msg = 'mixed boundary conditions allow rotation'
case (39_pInt)
msg = 'non-positive restart frequency in spectral loadcase'
case (40_pInt)
msg = 'path rectification error'
case (41_pInt)
msg = 'path too long'
case (42_pInt)
msg = 'missing header info in spectral mesh'
case (43_pInt)
msg = 'resolution in spectral mesh'
case (44_pInt)
msg = 'dimension in spectral mesh'
case (45_pInt)
msg = 'incomplete information in spectral mesh header'
case (46_pInt)
msg = 'not a rotation defined for loadcase rotation'
case (47_pInt)
msg = 'updating of gamma operator not possible if it is pre calculated'
case (50_pInt)
msg = 'writing constitutive output description'
!* file handling errors
case (100_pInt)
msg = 'opening material configuration'
msg = 'could not open file:'
case (101_pInt)
msg = 'opening input file'
msg = 'write error for file:'
case (102_pInt)
msg = 'non-positive dimension'
msg = 'could not read file:'
case (103_pInt)
msg = 'odd resolution given'
case (105_pInt)
msg = 'reading from ODF file'
case (106_pInt)
msg = 'reading info on old job'
case (107_pInt)
msg = 'writing spectralOut file'
case (108_pInt)
msg = 'precistion not suitable for FFTW'
case (109_pInt)
msg = 'initializing FFTW'
case (110_pInt)
msg = 'no homogenization specified via State Variable 2'
case (120_pInt)
msg = 'no microstructure specified via State Variable 3'
case (125_pInt)
msg = 'no entries in config part'
case (130_pInt)
msg = 'homogenization index out of bounds'
case (140_pInt)
msg = 'microstructure index out of bounds'
msg = 'could not assemble input files'
!* material error messages and related messages in mesh
case (150_pInt)
msg = 'crystallite index out of bounds'
case (155_pInt)
case (151_pInt)
msg = 'phase index out of bounds'
case (160_pInt)
case (152_pInt)
msg = 'texture index out of bounds'
case (170_pInt)
case (153_pInt)
msg = 'sum of phase fractions differs from 1'
case (154_pInt)
msg = 'homogenization index out of bounds'
case (155_pInt)
msg = 'microstructure index out of bounds'
case (156_pInt)
msg = 'reading from ODF file'
case (160_pInt)
msg = 'no entries in config part'
case (170_pInt)
msg = 'no homogenization specified via State Variable 2'
case (180_pInt)
msg = 'mismatch of microstructure count and a*b*c in geom file'
msg = 'no microstructure specified via State Variable 3'
!* constitutive error messages
case (200_pInt)
msg = 'unknown constitution specified'
case (201_pInt)
msg = 'unknown homogenization specified'
case (205_pInt)
msg = 'unknown lattice structure encountered'
case (210_pInt)
msg = 'negative initial resistance'
msg = 'unknown material parameter for j2 constitutive phase:'
case (211_pInt)
msg = 'non-positive reference shear rate'
msg = 'material parameter for j2 constitutive phase out of bounds:'
case (212_pInt)
msg = 'non-positive stress exponent'
case (213_pInt)
msg = 'non-positive saturation stress'
case (214_pInt)
msg = 'zero hardening exponent'
msg = 'unknown constitutive output for j2 constitution:'
case (220_pInt)
msg = 'negative initial dislocation density'
msg = 'unknown material parameter for phenopowerlaw constitutive phase:'
case (221_pInt)
msg = 'negative Burgers vector'
msg = 'material parameter for phenopowerlaw constitutive phase out of bounds:'
case (222_pInt)
msg = 'negative activation energy for edge dislocation glide'
case (223_pInt)
msg = 'zero stackin fault energy'
case (225_pInt)
msg = 'no slip systems specified'
case (226_pInt)
msg = 'non-positive prefactor for dislocation velocity'
case (227_pInt)
msg = 'non-positive prefactor for mean free path'
case (228_pInt)
msg = 'non-positive minimum stable dipole distance'
case (229_pInt)
msg = 'non-positive hardening interaction coefficients'
msg = 'unknown constitutive output for phenopowerlaw constitution:'
case (230_pInt)
msg = 'non-positive atomic volume'
msg = 'unknown material parameter for titanmod constitutive phase:'
case (231_pInt)
msg = 'non-positive prefactor for self-diffusion coefficient'
msg = 'material parameter for titanmod constitutive phase out of bounds:'
case (232_pInt)
msg = 'non-positive activation energy for self-diffusion'
case (233_pInt)
msg = 'non-positive relevant dislocation density'
case (234_pInt)
msg = 'error in shear banding input'
case (235_pInt)
msg = 'material parameter for nonlocal constitutive phase out of bounds:'
case (236_pInt)
msg = 'unknown material parameter for nonlocal constitutive phase:'
case (237_pInt)
msg = 'unknown constitutive output for nonlocal constitution:'
msg = 'unknown constitutive output for titanmod constitution:'
case (240_pInt)
msg = 'non-positive Taylor factor'
msg = 'unknown material parameter for dislotwin constitutive phase:'
case (241_pInt)
msg = 'non-positive hardening exponent'
msg = 'material parameter for dislotwin constitutive phase out of bounds:'
case (242_pInt)
msg = 'non-positive relevant slip resistance'
case (260_pInt)
msg = 'non-positive relevant strain'
case (261_pInt)
msg = 'frequency of stiffness update < 0'
case (262_pInt)
msg = 'frequency of Jacobian update of Lp residuum < 0'
case (263_pInt)
msg = 'non-positive perturbation value'
case (264_pInt)
msg = 'limit for homogenization loop too small'
case (265_pInt)
msg = 'limit for crystallite loop too small'
case (266_pInt)
msg = 'limit for state loop too small'
case (267_pInt)
msg = 'limit for stress loop too small'
case (268_pInt)
msg = 'non-positive minimum substep size'
case (269_pInt)
msg = 'non-positive relative state tolerance'
case (270_pInt)
msg = 'Non-positive relative stress tolerance'
case (271_pInt)
msg = 'Non-positive absolute stress tolerance'
msg = 'unknown constitutive output for dislotwin constitution:'
case (243_pInt)
msg = 'zero stackin fault energy'
!* Error messages related to RGC numerical parameters <<<updated 31.07.2009>>>
case (272_pInt)
msg = 'non-positive relative tolerance of residual in RGC'
case (273_pInt)
msg = 'non-positive absolute tolerance of residual in RGC'
case (274_pInt)
msg = 'non-positive relative maximum of residual in RGC'
case (275_pInt)
msg = 'non-positive absolute maximum of residual in RGC'
case (276_pInt)
msg = 'non-positive penalty perturbation in RGC'
case (277_pInt)
msg = 'non-positive relevant mismatch in RGC'
case (278_pInt)
msg = 'non-positive definite viscosity model in RGC'
case (288_pInt)
msg = 'non-positive maximum threshold of relaxation change in RGC'
case (289_pInt)
msg = 'non-positive definite volume discrepancy penalty in RGC'
case (250_pInt)
msg = 'unknown material parameter for nonlocal constitutive phase:'
case (251_pInt)
msg = 'material parameter for nonlocal constitutive phase out of bounds:'
case (252_pInt)
msg = 'unknown constitutive output for nonlocal constitution:'
case (253_pInt)
msg = 'element type not supported for nonlocal constitution'
case (294_pInt)
msg = 'non-positive tolerance for deformation gradient'
case (298_pInt)
msg = 'chosen integration method does not exist'
case (299_pInt)
msg = 'chosen perturbation method does not exist'
!* numerics error messages
case (300_pInt)
msg = 'this material can only be used with elements with three direct stress components'
case (500_pInt)
msg = 'unknown lattice type specified'
case (550_pInt)
msg = 'unknown numerics parameter:'
case (301_pInt)
msg = 'numerics parameter out of bounds:'
!* math errors
case (400_pInt)
msg = 'matrix inversion error'
case (401_pInt)
msg = 'math_check: quat -> axisAngle -> quat failed'
case (402_pInt)
msg = 'math_check: quat -> R -> quat failed'
case (403_pInt)
msg = 'math_check: quat -> euler -> quat failed'
case (404_pInt)
msg = 'math_check: R -> euler -> R failed'
case (405_pInt)
msg = 'I_TO_HALTON-error: An input base BASE is <= 1'
case (406_pInt)
msg = 'Prime-error: N must be between 0 and PRIME_MAX'
case (450_pInt)
msg = 'unknown symmetry type specified'
case (600_pInt)
msg = 'convergence not reached'
case (610_pInt)
msg = 'stress loop not converged'
case (666_pInt)
msg = 'memory leak detected'
case (667_pInt)
!* homogenization errors
case (500_pInt)
msg = 'unknown homogenization specified'
!* DAMASK_marc errors
case (700_pInt)
msg = 'invalid materialpoint result requested'
case (670_pInt)
msg = 'math_check: quat -> axisAngle -> quat failed'
case (671_pInt)
msg = 'math_check: quat -> R -> quat failed'
case (672_pInt)
msg = 'math_check: quat -> euler -> quat failed'
case (673_pInt)
msg = 'math_check: R -> euler -> R failed'
case (700_pInt)
msg = 'singular matrix in stress iteration'
!* errors related to spectral solver
case (800_pInt)
msg = 'matrix inversion error'
case (801_pInt)
msg = 'I_TO_HALTON-error: An input base BASE is <= 1'
case (802_pInt)
msg = 'Prime-error: N must be between 0 and PRIME_MAX'
msg = 'non-positive dimension'
case (803_pInt)
msg = 'odd resolution given'
case (808_pInt)
msg = 'precision not suitable for FFTW'
case (809_pInt)
msg = 'initializing FFTW'
case (831_pInt)
msg = 'mask consistency violated in spectral loadcase'
case (832_pInt)
msg = 'ill-defined L (each line should be either fully or not at all defined) in spectral loadcase'
case (834_pInt)
msg = 'negative time increment in spectral loadcase'
case (835_pInt)
msg = 'non-positive increments in spectral loadcase'
case (836_pInt)
msg = 'non-positive result frequency in spectral loadcase'
case (837_pInt)
msg = 'incomplete loadcase'
case (838_pInt)
msg = 'mixed boundary conditions allow rotation'
case (842_pInt)
msg = 'missing header info in spectral mesh'
case (843_pInt)
msg = 'resolution in spectral mesh'
case (844_pInt)
msg = 'dimension in spectral mesh'
case (845_pInt)
msg = 'incomplete information in spectral mesh header'
case (846_pInt)
msg = 'not a rotation defined for loadcase rotation'
case (847_pInt)
msg = 'updating of gamma operator not possible if it is pre calculated'
case (880_pInt)
msg = 'mismatch of microstructure count and a*b*c in geom file'
!* Error messages related to parsing of Abaqus input file
! Error messages related to parsing of Abaqus input file
case (900_pInt)
msg = 'PARSE ERROR: Improper definition of nodes in input file (Nnodes < 2)'
case (901_pInt)
@ -1378,8 +1397,13 @@ endfunction
msg = 'PARSE ERROR: Incorrect element type mapping in '
!* general error messages
case (666_pInt)
msg = 'memory leak detected'
case default
msg = 'Unknown error number...'
end select
!$OMP CRITICAL (write2out)

10
code/constitutive.f90
View File

@ -71,7 +71,7 @@ subroutine constitutive_init()
use prec, only: pReal,pInt
use debug, only: debug_verbosity, debug_selectiveDebugger, debug_e, debug_i, debug_g
use numerics, only: numerics_integrator
use IO, only: IO_error, IO_open_file, IO_open_jobFile, IO_write_jobFile
use IO, only: IO_error, IO_open_file, IO_open_jobFile_stat, IO_write_jobFile
use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips
use material
use constitutive_j2
@ -100,8 +100,8 @@ logical knownConstitution
! --- PARSE CONSTITUTIONS FROM CONFIG FILE ---
if (.not. IO_open_jobFile(fileunit,material_localFileExt)) then ! no local material configuration present...
if (.not. IO_open_file(fileunit,material_configFile)) call IO_error(100_pInt) ! ...and cannot open material.config file
if (.not. IO_open_jobFile_stat(fileunit,material_localFileExt)) then ! no local material configuration present...
call IO_open_file(fileunit,material_configFile) ! ... open material.config file
endif
call constitutive_j2_init(fileunit)
call constitutive_phenopowerlaw_init(fileunit)
@ -113,9 +113,7 @@ close(fileunit)
! --- WRITE DESCRIPTION FILE FOR CONSTITUTIVE PHASE OUTPUT ---
if(.not. IO_write_jobFile(fileunit,'outputConstitutive')) then ! problems in writing file
call IO_error (50_pInt)
endif
call IO_write_jobFile(fileunit,'outputConstitutive')
do p = 1,material_Nphase
i = phase_constitutionInstance(p) ! which instance of a constitution is present phase
knownConstitution = .true. ! assume valid

39
code/constitutive_dislotwin.f90
View File

@ -393,6 +393,8 @@ do ! read thru sections of
constitutive_dislotwin_sbResistance(i) = IO_floatValue(line,positions,2_pInt)
case ('shearbandvelocity')
constitutive_dislotwin_sbVelocity(i) = IO_floatValue(line,positions,2_pInt)
case default
call IO_error(240_pInt,ext_msg=tag)
end select
endif
enddo
@ -403,32 +405,31 @@ enddo
myStructure = constitutive_dislotwin_structure(i)
!* Sanity checks
if (myStructure < 1_pInt .or. myStructure > 3_pInt) call IO_error(205_pInt)
if (sum(constitutive_dislotwin_Nslip(:,i)) <= 0_pInt) call IO_error(225_pInt)
if (sum(constitutive_dislotwin_Ntwin(:,i)) < 0_pInt) call IO_error(225_pInt) !***
if (myStructure < 1_pInt .or. myStructure > 3_pInt) call IO_error(205_pInt,e=i)
if (sum(constitutive_dislotwin_Nslip(:,i)) <= 0_pInt) call IO_error(241_pInt,e=i,ext_msg='nslip')
if (sum(constitutive_dislotwin_Ntwin(:,i)) < 0_pInt) call IO_error(241_pInt,e=i,ext_msg='ntwin')
do f = 1,lattice_maxNslipFamily
if (constitutive_dislotwin_Nslip(f,i) > 0_pInt) then
if (constitutive_dislotwin_rhoEdge0(f,i) < 0.0_pReal) call IO_error(220_pInt)
if (constitutive_dislotwin_rhoEdgeDip0(f,i) < 0.0_pReal) call IO_error(220_pInt)
if (constitutive_dislotwin_burgersPerSlipFamily(f,i) <= 0.0_pReal) call IO_error(221_pInt)
if (constitutive_dislotwin_v0PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(226_pInt)
if (constitutive_dislotwin_rhoEdge0(f,i) < 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='rhoEdge0')
if (constitutive_dislotwin_rhoEdgeDip0(f,i) < 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='rhoEdgeDip0')
if (constitutive_dislotwin_burgersPerSlipFamily(f,i) <= 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='slipburgers')
if (constitutive_dislotwin_v0PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='v0')
endif
enddo
do f = 1,lattice_maxNtwinFamily
if (constitutive_dislotwin_Ntwin(f,i) > 0_pInt) then
if (constitutive_dislotwin_burgersPerTwinFamily(f,i) <= 0.0_pReal) call IO_error(221_pInt) !***
if (constitutive_dislotwin_Ndot0PerTwinFamily(f,i) < 0.0_pReal) call IO_error(226_pInt) !***
if (constitutive_dislotwin_burgersPerTwinFamily(f,i) <= 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='twinburgers')
if (constitutive_dislotwin_Ndot0PerTwinFamily(f,i) < 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='ndot0')
endif
enddo
! if (any(constitutive_dislotwin_interactionSlipSlip(1:maxval(lattice_interactionSlipSlip(:,:,myStructure)),i) < 1.0_pReal)) call IO_error(229)
if (constitutive_dislotwin_CAtomicVolume(i) <= 0.0_pReal) call IO_error(230_pInt)
if (constitutive_dislotwin_D0(i) <= 0.0_pReal) call IO_error(231_pInt)
if (constitutive_dislotwin_Qsd(i) <= 0.0_pReal) call IO_error(232_pInt)
if (constitutive_dislotwin_aTolRho(i) <= 0.0_pReal) call IO_error(233_pInt)
if (constitutive_dislotwin_sbResistance(i) <= 0.0_pReal) call IO_error(234_pInt)
if (constitutive_dislotwin_sbVelocity(i) < 0.0_pReal) call IO_error(235_pInt)
if (constitutive_dislotwin_SFE_0K(i) == 0.0_pReal .and. &
constitutive_dislotwin_dSFE_dT(i) == 0.0_pReal) call IO_error(223_pInt)
if (constitutive_dislotwin_CAtomicVolume(i) <= 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='cAtomicVolume')
if (constitutive_dislotwin_D0(i) <= 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='D0')
if (constitutive_dislotwin_Qsd(i) <= 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='Qsd')
if (constitutive_dislotwin_aTolRho(i) <= 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='aTolRho')
if (constitutive_dislotwin_sbResistance(i) <= 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='sbResistance')
if (constitutive_dislotwin_sbVelocity(i) < 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='sbVelocity')
if (constitutive_dislotwin_SFE_0K(i) == 0.0_pReal .AND. &
constitutive_dislotwin_dSFE_dT(i) == 0.0_pReal) call IO_error(243_pInt,e=i)
!* Determine total number of active slip or twin systems
constitutive_dislotwin_Nslip(:,i) = min(lattice_NslipSystem(:,myStructure),constitutive_dislotwin_Nslip(:,i))
@ -534,7 +535,7 @@ do i = 1_pInt,maxNinstance
case('sb_eigenvectors')
mySize = 9_pInt
case default
mySize = 0_pInt
call IO_error(242_pInt,ext_msg=constitutive_dislotwin_output(o,i))
end select
if (mySize > 0_pInt) then ! any meaningful output found

89
code/constitutive_j2.f90
View File

@ -83,15 +83,14 @@ subroutine constitutive_j2_init(file)
!**************************************
!* Module initialization *
!**************************************
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: pInt, pReal
use math, only: math_Mandel3333to66, math_Voigt66to3333
use IO
use material
use debug, only: debug_verbosity
integer(pInt), intent(in) :: file
integer(pInt), parameter :: maxNchunks = 7_pInt
integer(pInt), dimension(1+2_pInt*maxNchunks) :: positions
integer(pInt), parameter :: maxNchunks = 7
integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt) section, maxNinstance, i,j,k, output, mySize
character(len=64) tag
character(len=1024) line
@ -104,9 +103,9 @@ subroutine constitutive_j2_init(file)
!$OMP END CRITICAL (write2out)
maxNinstance = count(phase_constitution == constitutive_j2_label)
if (maxNinstance == 0_pInt) return
if (maxNinstance == 0) return
if (debug_verbosity > 0_pInt) then
if (debug_verbosity > 0) then
!$OMP CRITICAL (write2out)
write(6,'(a16,1x,i5)') '# instances:',maxNinstance
write(6,*)
@ -127,12 +126,12 @@ subroutine constitutive_j2_init(file)
allocate(constitutive_j2_n(maxNinstance)) ; constitutive_j2_n = 0.0_pReal
allocate(constitutive_j2_h0(maxNinstance)) ; constitutive_j2_h0 = 0.0_pReal
allocate(constitutive_j2_tausat(maxNinstance)) ; constitutive_j2_tausat = 0.0_pReal
allocate(constitutive_j2_a(maxNinstance)) ; constitutive_j2_a = 0.0_pReal
allocate(constitutive_j2_a(maxNinstance)) ; constitutive_j2_a = 0.0_pReal
allocate(constitutive_j2_aTolResistance(maxNinstance)) ; constitutive_j2_aTolResistance = 0.0_pReal
rewind(file)
line = ''
section = 0_pInt
section = 0
do while (IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to <phase>
read(file,'(a1024)',END=100) line
@ -143,60 +142,62 @@ subroutine constitutive_j2_init(file)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1_pInt
output = 0_pInt ! reset output counter
section = section + 1
output = 0 ! reset output counter
endif
if (section > 0_pInt .and. phase_constitution(section) == constitutive_j2_label) then ! one of my sections
if (section > 0 .and. phase_constitution(section) == constitutive_j2_label) then ! one of my sections
i = phase_constitutionInstance(section) ! which instance of my constitution is present phase
positions = IO_stringPos(line,maxNchunks)
tag = IO_lc(IO_stringValue(line,positions,1)) ! extract key
select case(tag)
case ('(output)')
output = output + 1_pInt
constitutive_j2_output(output,i) = IO_lc(IO_stringValue(line,positions,2_pInt))
output = output + 1
constitutive_j2_output(output,i) = IO_lc(IO_stringValue(line,positions,2))
case ('c11')
constitutive_j2_C11(i) = IO_floatValue(line,positions,2_pInt)
constitutive_j2_C11(i) = IO_floatValue(line,positions,2)
case ('c12')
constitutive_j2_C12(i) = IO_floatValue(line,positions,2_pInt)
constitutive_j2_C12(i) = IO_floatValue(line,positions,2)
case ('tau0')
constitutive_j2_tau0(i) = IO_floatValue(line,positions,2_pInt)
constitutive_j2_tau0(i) = IO_floatValue(line,positions,2)
case ('gdot0')
constitutive_j2_gdot0(i) = IO_floatValue(line,positions,2_pInt)
constitutive_j2_gdot0(i) = IO_floatValue(line,positions,2)
case ('n')
constitutive_j2_n(i) = IO_floatValue(line,positions,2_pInt)
constitutive_j2_n(i) = IO_floatValue(line,positions,2)
case ('h0')
constitutive_j2_h0(i) = IO_floatValue(line,positions,2_pInt)
constitutive_j2_h0(i) = IO_floatValue(line,positions,2)
case ('tausat')
constitutive_j2_tausat(i) = IO_floatValue(line,positions,2_pInt)
constitutive_j2_tausat(i) = IO_floatValue(line,positions,2)
case ('a', 'w0')
constitutive_j2_a(i) = IO_floatValue(line,positions,2_pInt)
constitutive_j2_a(i) = IO_floatValue(line,positions,2)
case ('taylorfactor')
constitutive_j2_fTaylor(i) = IO_floatValue(line,positions,2_pInt)
constitutive_j2_fTaylor(i) = IO_floatValue(line,positions,2)
case ('atol_resistance')
constitutive_j2_aTolResistance(i) = IO_floatValue(line,positions,2_pInt)
constitutive_j2_aTolResistance(i) = IO_floatValue(line,positions,2)
case default
call IO_error(210_pInt,ext_msg=tag)
end select
endif
enddo
100 do i = 1_pInt,maxNinstance ! sanity checks
if (constitutive_j2_tau0(i) < 0.0_pReal) call IO_error(210_pInt)
if (constitutive_j2_gdot0(i) <= 0.0_pReal) call IO_error(211_pInt)
if (constitutive_j2_n(i) <= 0.0_pReal) call IO_error(212_pInt)
if (constitutive_j2_tausat(i) <= 0.0_pReal) call IO_error(213_pInt)
if (constitutive_j2_a(i) <= 0.0_pReal) call IO_error(241_pInt)
if (constitutive_j2_fTaylor(i) <= 0.0_pReal) call IO_error(240_pInt)
if (constitutive_j2_aTolResistance(i) <= 0.0_pReal) call IO_error(242_pInt)
100 do i = 1,maxNinstance ! sanity checks
if (constitutive_j2_tau0(i) < 0.0_pReal) call IO_error(211_pInt,ext_msg='tau0')
if (constitutive_j2_gdot0(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='gdot0')
if (constitutive_j2_n(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='n')
if (constitutive_j2_tausat(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='tausat')
if (constitutive_j2_a(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='a')
if (constitutive_j2_fTaylor(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='taylorfactor')
if (constitutive_j2_aTolResistance(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='aTol_resistance')
enddo
do i = 1_pInt,maxNinstance
do j = 1_pInt,maxval(phase_Noutput)
do i = 1,maxNinstance
do j = 1,maxval(phase_Noutput)
select case(constitutive_j2_output(j,i))
case('flowstress')
mySize = 1_pInt
mySize = 1_pInt
case('strainrate')
mySize = 1_pInt
mySize = 1_pInt
case default
mySize = 0_pInt
call IO_error(212_pInt,ext_msg=constitutive_j2_output(j,i))
end select
if (mySize > 0_pInt) then ! any meaningful output found
@ -206,24 +207,22 @@ subroutine constitutive_j2_init(file)
endif
enddo
constitutive_j2_sizeDotState(i) = 1_pInt
constitutive_j2_sizeState(i) = 1_pInt
constitutive_j2_sizeDotState(i) = 1
constitutive_j2_sizeState(i) = 1
forall(k=1_pInt:3_pInt)
forall(j=1_pInt:3_pInt) &
forall(k=1:3)
forall(j=1:3) &
constitutive_j2_Cslip_66(k,j,i) = constitutive_j2_C12(i)
constitutive_j2_Cslip_66(k,k,i) = constitutive_j2_C11(i)
constitutive_j2_Cslip_66(k+3_pInt,k+3_pInt,i) = 0.5_pReal*(constitutive_j2_C11(i)-constitutive_j2_C12(i))
constitutive_j2_Cslip_66(k+3,k+3,i) = 0.5_pReal*(constitutive_j2_C11(i)-constitutive_j2_C12(i))
end forall
constitutive_j2_Cslip_66(1:6,1:6,i) = &
math_Mandel3333to66(math_Voigt66to3333(constitutive_j2_Cslip_66(1:6,1:6,i)))
enddo
return
endsubroutine
@ -244,7 +243,6 @@ pure function constitutive_j2_stateInit(myInstance)
endfunction
!*********************************************************************
!* relevant microstructural state *
!*********************************************************************
@ -265,7 +263,6 @@ real(pReal), dimension(constitutive_j2_sizeState(myInstance)) :: &
constitutive_j2_aTolState = constitutive_j2_aTolResistance(myInstance)
endfunction
@ -293,7 +290,6 @@ function constitutive_j2_homogenizedC(state,ipc,ip,el)
return
endfunction
@ -316,7 +312,6 @@ subroutine constitutive_j2_microstructure(Temperature,state,ipc,ip,el)
real(pReal) Temperature
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
matID = phase_constitutionInstance(material_phase(ipc,ip,el))
endsubroutine
@ -342,7 +337,6 @@ pure subroutine constitutive_j2_LpAndItsTangent(Lp, dLp_dTstar_99, Tstar_dev_v,
material_phase, &
phase_constitutionInstance
implicit none
!*** input variables ***!
@ -466,7 +460,6 @@ pure function constitutive_j2_dotState(Tstar_v, Temperature, state, g, ip, el)
endfunction
!****************************************************************
!* calculates the rate of change of temperature *
!****************************************************************

138
code/constitutive_nonlocal.f90
View File

@ -341,13 +341,13 @@ do
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1_pInt
output = 0_pInt ! reset output counter
output = 0_pInt ! reset output counter
cycle
endif
if (section > 0_pInt .and. phase_constitution(section) == constitutive_nonlocal_label) then ! one of my sections
if (section > 0_pInt .and. phase_constitution(section) == constitutive_nonlocal_label) then ! one of my sections
i = phase_constitutionInstance(section) ! which instance of my constitution is present phase
positions = IO_stringPos(line,maxNchunks)
tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key
tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key
select case(tag)
case('constitution','/nonlocal/')
cycle
@ -369,40 +369,40 @@ do
case ('c44')
constitutive_nonlocal_C44(i) = IO_floatValue(line,positions,2_pInt)
case ('nslip')
forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_Nslip(f,i)&
= IO_intValue(line,positions,1_pInt+f)
forall (f = 1_pInt:lattice_maxNslipFamily) &
constitutive_nonlocal_Nslip(f,i) = IO_intValue(line,positions,1_pInt+f)
case ('rhosgledgepos0')
forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_rhoSglEdgePos0(f,i)&
= IO_floatValue(line,positions,1_pInt+f)
forall (f = 1_pInt:lattice_maxNslipFamily) &
constitutive_nonlocal_rhoSglEdgePos0(f,i) = IO_floatValue(line,positions,1_pInt+f)
case ('rhosgledgeneg0')
forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_rhoSglEdgeNeg0(f,i)&
= IO_floatValue(line,positions,1_pInt+f)
forall (f = 1_pInt:lattice_maxNslipFamily) &
constitutive_nonlocal_rhoSglEdgeNeg0(f,i) = IO_floatValue(line,positions,1_pInt+f)
case ('rhosglscrewpos0')
forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_rhoSglScrewPos0(f,i)&
= IO_floatValue(line,positions,1_pInt+f)
forall (f = 1_pInt:lattice_maxNslipFamily) &
constitutive_nonlocal_rhoSglScrewPos0(f,i) = IO_floatValue(line,positions,1_pInt+f)
case ('rhosglscrewneg0')
forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_rhoSglScrewNeg0(f,i)&
= IO_floatValue(line,positions,1_pInt+f)
forall (f = 1_pInt:lattice_maxNslipFamily) &
constitutive_nonlocal_rhoSglScrewNeg0(f,i) = IO_floatValue(line,positions,1_pInt+f)
case ('rhodipedge0')
forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_rhoDipEdge0(f,i)&
= IO_floatValue(line,positions,1_pInt+f)
forall (f = 1_pInt:lattice_maxNslipFamily) &
constitutive_nonlocal_rhoDipEdge0(f,i) = IO_floatValue(line,positions,1_pInt+f)
case ('rhodipscrew0')
forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_rhoDipScrew0(f,i)&
= IO_floatValue(line,positions,1_pInt+f)
forall (f = 1_pInt:lattice_maxNslipFamily) &
constitutive_nonlocal_rhoDipScrew0(f,i) = IO_floatValue(line,positions,1_pInt+f)
case ('lambda0')
forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_lambda0PerSlipFamily(f,i)&
= IO_floatValue(line,positions,1_pInt+f)
forall (f = 1_pInt:lattice_maxNslipFamily) &
constitutive_nonlocal_lambda0PerSlipFamily(f,i) = IO_floatValue(line,positions,1_pInt+f)
case ('burgers')
forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_burgersPerSlipFamily(f,i)&
= IO_floatValue(line,positions,1_pInt+f)
forall (f = 1_pInt:lattice_maxNslipFamily) &
constitutive_nonlocal_burgersPerSlipFamily(f,i) = IO_floatValue(line,positions,1_pInt+f)
case('cutoffradius','r')
constitutive_nonlocal_R(i) = IO_floatValue(line,positions,2_pInt)
case('minimumdipoleheightedge','ddipminedge')
forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_minimumDipoleHeightPerSlipFamily(f,1_pInt,i)&
= IO_floatValue(line,positions,1_pInt+f)
forall (f = 1_pInt:lattice_maxNslipFamily) &
constitutive_nonlocal_minimumDipoleHeightPerSlipFamily(f,1_pInt,i) = IO_floatValue(line,positions,1_pInt+f)
case('minimumdipoleheightscrew','ddipminscrew')
forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_minimumDipoleHeightPerSlipFamily(f,2_pInt,i)&
= IO_floatValue(line,positions,1_pInt+f)
forall (f = 1_pInt:lattice_maxNslipFamily) &
constitutive_nonlocal_minimumDipoleHeightPerSlipFamily(f,2_pInt,i) = IO_floatValue(line,positions,1_pInt+f)
case('atomicvolume')
constitutive_nonlocal_atomicVolume(i) = IO_floatValue(line,positions,2_pInt)
case('selfdiffusionprefactor','dsd0')
@ -412,14 +412,14 @@ do
case('atol_rho')
constitutive_nonlocal_aTolRho(i) = IO_floatValue(line,positions,2_pInt)
case ('interaction_slipslip')
forall (it = 1_pInt:lattice_maxNinteraction) constitutive_nonlocal_interactionSlipSlip(it,i)&
= IO_floatValue(line,positions,1_pInt+it)
forall (it = 1_pInt:lattice_maxNinteraction) &
constitutive_nonlocal_interactionSlipSlip(it,i) = IO_floatValue(line,positions,1_pInt+it)
case('peierlsstressedge')
forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_peierlsStressPerSlipFamily(f,1_pInt,i)&
= IO_floatValue(line,positions,1_pInt+f)
forall (f = 1_pInt:lattice_maxNslipFamily) &
constitutive_nonlocal_peierlsStressPerSlipFamily(f,1_pInt,i) = IO_floatValue(line,positions,1_pInt+f)
case('peierlsstressscrew')
forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_peierlsStressPerSlipFamily(f,2_pInt,i)&
= IO_floatValue(line,positions,1_pInt+f)
forall (f = 1_pInt:lattice_maxNslipFamily) &
constitutive_nonlocal_peierlsStressPerSlipFamily(f,2_pInt,i) = IO_floatValue(line,positions,1_pInt+f)
case('doublekinkwidth')
constitutive_nonlocal_doublekinkwidth(i) = IO_floatValue(line,positions,2_pInt)
case('solidsolutionenergy')
@ -441,7 +441,7 @@ do
case('surfacetransmissivity')
constitutive_nonlocal_surfaceTransmissivity(i) = IO_floatValue(line,positions,2_pInt)
case default
call IO_error(236_pInt,ext_msg=tag)
call IO_error(250_pInt,ext_msg=tag)
end select
endif
enddo
@ -457,46 +457,46 @@ enddo
!*** sanity checks
if (myStructure < 1 .or. myStructure > 3) call IO_error(205_pInt)
if (sum(constitutive_nonlocal_Nslip(:,i)) <= 0_pInt) call IO_error(235_pInt,ext_msg='Nslip')
if (sum(constitutive_nonlocal_Nslip(:,i)) <= 0_pInt) call IO_error(251_pInt,ext_msg='Nslip')
do o = 1,maxval(phase_Noutput)
if(len(constitutive_nonlocal_output(o,i)) > 64) call IO_error(666_pInt)
enddo
do f = 1,lattice_maxNslipFamily
if (constitutive_nonlocal_Nslip(f,i) > 0_pInt) then
if (constitutive_nonlocal_rhoSglEdgePos0(f,i) < 0.0_pReal) call IO_error(235_pInt,ext_msg='rhoSglEdgePos0')
if (constitutive_nonlocal_rhoSglEdgeNeg0(f,i) < 0.0_pReal) call IO_error(235_pInt,ext_msg='rhoSglEdgeNeg0')
if (constitutive_nonlocal_rhoSglScrewPos0(f,i) < 0.0_pReal) call IO_error(235_pInt,ext_msg='rhoSglScrewPos0')
if (constitutive_nonlocal_rhoSglScrewNeg0(f,i) < 0.0_pReal) call IO_error(235_pInt,ext_msg='rhoSglScrewNeg0')
if (constitutive_nonlocal_rhoDipEdge0(f,i) < 0.0_pReal) call IO_error(235_pInt,ext_msg='rhoDipEdge0')
if (constitutive_nonlocal_rhoDipScrew0(f,i) < 0.0_pReal) call IO_error(235_pInt,ext_msg='rhoDipScrew0')
if (constitutive_nonlocal_burgersPerSlipFamily(f,i) <= 0.0_pReal) call IO_error(235_pInt,ext_msg='burgers')
if (constitutive_nonlocal_lambda0PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(235_pInt,ext_msg='lambda0')
if (constitutive_nonlocal_rhoSglEdgePos0(f,i) < 0.0_pReal) call IO_error(251_pInt,ext_msg='rhoSglEdgePos0')
if (constitutive_nonlocal_rhoSglEdgeNeg0(f,i) < 0.0_pReal) call IO_error(251_pInt,ext_msg='rhoSglEdgeNeg0')
if (constitutive_nonlocal_rhoSglScrewPos0(f,i) < 0.0_pReal) call IO_error(251_pInt,ext_msg='rhoSglScrewPos0')
if (constitutive_nonlocal_rhoSglScrewNeg0(f,i) < 0.0_pReal) call IO_error(251_pInt,ext_msg='rhoSglScrewNeg0')
if (constitutive_nonlocal_rhoDipEdge0(f,i) < 0.0_pReal) call IO_error(251_pInt,ext_msg='rhoDipEdge0')
if (constitutive_nonlocal_rhoDipScrew0(f,i) < 0.0_pReal) call IO_error(251_pInt,ext_msg='rhoDipScrew0')
if (constitutive_nonlocal_burgersPerSlipFamily(f,i) <= 0.0_pReal) call IO_error(251_pInt,ext_msg='burgers')
if (constitutive_nonlocal_lambda0PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(251_pInt,ext_msg='lambda0')
if (constitutive_nonlocal_minimumDipoleHeightPerSlipFamily(f,1,i) <= 0.0_pReal) &
call IO_error(235_pInt,ext_msg='minimumDipoleHeightEdge')
call IO_error(251_pInt,ext_msg='minimumDipoleHeightEdge')
if (constitutive_nonlocal_minimumDipoleHeightPerSlipFamily(f,2,i) <= 0.0_pReal) &
call IO_error(235_pInt,ext_msg='minimumDipoleHeightScrew')
if (constitutive_nonlocal_peierlsStressPerSlipFamily(f,1,i) <= 0.0_pReal) call IO_error(235_pInt,ext_msg='peierlsStressEdge')
if (constitutive_nonlocal_peierlsStressPerSlipFamily(f,2,i) <= 0.0_pReal) call IO_error(235_pInt,ext_msg='peierlsStressScrew')
call IO_error(251_pInt,ext_msg='minimumDipoleHeightScrew')
if (constitutive_nonlocal_peierlsStressPerSlipFamily(f,1,i) <= 0.0_pReal) call IO_error(251_pInt,ext_msg='peierlsStressEdge')
if (constitutive_nonlocal_peierlsStressPerSlipFamily(f,2,i) <= 0.0_pReal) call IO_error(251_pInt,ext_msg='peierlsStressScrew')
endif
enddo
if (any(constitutive_nonlocal_interactionSlipSlip(1:maxval(lattice_interactionSlipSlip(:,:,myStructure)),i) < 0.0_pReal)) &
call IO_error(235_pInt,ext_msg='interaction_SlipSlip')
if (constitutive_nonlocal_R(i) < 0.0_pReal) call IO_error(235_pInt,ext_msg='r')
if (constitutive_nonlocal_atomicVolume(i) <= 0.0_pReal) call IO_error(235_pInt,ext_msg='atomicVolume')
if (constitutive_nonlocal_Dsd0(i) <= 0.0_pReal) call IO_error(235_pInt,ext_msg='selfDiffusionPrefactor')
if (constitutive_nonlocal_Qsd(i) <= 0.0_pReal) call IO_error(235_pInt,ext_msg='selfDiffusionEnergy')
if (constitutive_nonlocal_aTolRho(i) <= 0.0_pReal) call IO_error(235_pInt,ext_msg='aTol_rho')
if (constitutive_nonlocal_doublekinkwidth(i) <= 0.0_pReal) call IO_error(235_pInt,ext_msg='doublekinkwidth')
if (constitutive_nonlocal_solidSolutionEnergy(i) <= 0.0_pReal) call IO_error(235_pInt,ext_msg='solidSolutionEnergy')
if (constitutive_nonlocal_solidSolutionSize(i) <= 0.0_pReal) call IO_error(235_pInt,ext_msg='solidSolutionSize')
if (constitutive_nonlocal_solidSolutionConcentration(i) <= 0.0_pReal) call IO_error(235_pInt,ext_msg='solidSolutionConcentration')
if (constitutive_nonlocal_p(i) <= 0.0_pReal .or. constitutive_nonlocal_p(i) > 1.0_pReal) call IO_error(235_pInt,ext_msg='p')
if (constitutive_nonlocal_q(i) < 1.0_pReal .or. constitutive_nonlocal_q(i) > 2.0_pReal) call IO_error(235_pInt,ext_msg='q')
if (constitutive_nonlocal_viscosity(i) <= 0.0_pReal) call IO_error(235_pInt,ext_msg='viscosity')
if (constitutive_nonlocal_fattack(i) <= 0.0_pReal) call IO_error(235_pInt,ext_msg='attackFrequency')
if (constitutive_nonlocal_rhoSglScatter(i) < 0.0_pReal) call IO_error(235_pInt,ext_msg='rhoSglScatter')
call IO_error(251_pInt,ext_msg='interaction_SlipSlip')
if (constitutive_nonlocal_R(i) < 0.0_pReal) call IO_error(251_pInt,ext_msg='r')
if (constitutive_nonlocal_atomicVolume(i) <= 0.0_pReal) call IO_error(251_pInt,ext_msg='atomicVolume')
if (constitutive_nonlocal_Dsd0(i) <= 0.0_pReal) call IO_error(251_pInt,ext_msg='selfDiffusionPrefactor')
if (constitutive_nonlocal_Qsd(i) <= 0.0_pReal) call IO_error(251_pInt,ext_msg='selfDiffusionEnergy')
if (constitutive_nonlocal_aTolRho(i) <= 0.0_pReal) call IO_error(251_pInt,ext_msg='aTol_rho')
if (constitutive_nonlocal_doublekinkwidth(i) <= 0.0_pReal) call IO_error(251_pInt,ext_msg='doublekinkwidth')
if (constitutive_nonlocal_solidSolutionEnergy(i) <= 0.0_pReal) call IO_error(251_pInt,ext_msg='solidSolutionEnergy')
if (constitutive_nonlocal_solidSolutionSize(i) <= 0.0_pReal) call IO_error(251_pInt,ext_msg='solidSolutionSize')
if (constitutive_nonlocal_solidSolutionConcentration(i) <= 0.0_pReal) call IO_error(251_pInt,ext_msg='solidSolutionConcentration')
if (constitutive_nonlocal_p(i) <= 0.0_pReal .or. constitutive_nonlocal_p(i) > 1.0_pReal) call IO_error(251_pInt,ext_msg='p')
if (constitutive_nonlocal_q(i) < 1.0_pReal .or. constitutive_nonlocal_q(i) > 2.0_pReal) call IO_error(251_pInt,ext_msg='q')
if (constitutive_nonlocal_viscosity(i) <= 0.0_pReal) call IO_error(251_pInt,ext_msg='viscosity')
if (constitutive_nonlocal_fattack(i) <= 0.0_pReal) call IO_error(251_pInt,ext_msg='attackFrequency')
if (constitutive_nonlocal_rhoSglScatter(i) < 0.0_pReal) call IO_error(251_pInt,ext_msg='rhoSglScatter')
if (constitutive_nonlocal_surfaceTransmissivity(i) < 0.0_pReal &
.or. constitutive_nonlocal_surfaceTransmissivity(i) > 1.0_pReal) call IO_error(235_pInt,ext_msg='surfaceTransmissivity')
.or. constitutive_nonlocal_surfaceTransmissivity(i) > 1.0_pReal) call IO_error(251_pInt,ext_msg='surfaceTransmissivity')
!*** determine total number of active slip systems
@ -649,7 +649,7 @@ do i = 1,maxNinstance
case('dislocationstress')
mySize = 6_pInt
case default
call IO_error(237,ext_msg=constitutive_nonlocal_output(o,i))
call IO_error(252_pInt,ext_msg=constitutive_nonlocal_output(o,i))
end select
if (mySize > 0_pInt) then ! any meaningful output found
@ -662,13 +662,13 @@ do i = 1,maxNinstance
!*** elasticity matrix and shear modulus according to material.config
select case (myStructure)
case(1_pInt:2_pInt) ! cubic(s)
case(1_pInt:2_pInt) ! cubic(s)
forall(k=1_pInt:3_pInt)
forall(j=1_pInt:3_pInt) constitutive_nonlocal_Cslip_66(k,j,i) = constitutive_nonlocal_C12(i)
constitutive_nonlocal_Cslip_66(k,k,i) = constitutive_nonlocal_C11(i)
constitutive_nonlocal_Cslip_66(k+3_pInt,k+3_pInt,i) = constitutive_nonlocal_C44(i)
end forall
case(3_pInt:) ! all hex
case(3_pInt:) ! all hex
constitutive_nonlocal_Cslip_66(1,1,i) = constitutive_nonlocal_C11(i)
constitutive_nonlocal_Cslip_66(2,2,i) = constitutive_nonlocal_C11(i)
constitutive_nonlocal_Cslip_66(3,3,i) = constitutive_nonlocal_C33(i)
@ -1646,7 +1646,8 @@ logical considerEnteringFlux, &
considerLeavingFlux
#ifndef _OPENMP
if (debug_verbosity > 6_pInt .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then
if (debug_verbosity > 6_pInt .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) &
.or. .not. debug_selectiveDebugger)) then
write(6,*)
write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_dotState at el ip g ',el,ip,g
write(6,*)
@ -1657,7 +1658,7 @@ select case(mesh_element(2,el))
case (1_pInt,6_pInt,7_pInt,8_pInt,9_pInt)
! all fine
case default
call IO_error(-1_pInt,el,ip,g,'element type not supported for nonlocal constitution')
call IO_error(253_pInt,el,ip,g)
end select
myInstance = phase_constitutionInstance(material_phase(g,ip,el))
@ -1699,11 +1700,12 @@ endif
forall (t = 1_pInt:4_pInt) &
gdot(1_pInt:ns,t) = rhoSgl(1_pInt:ns,t) * constitutive_nonlocal_burgers(1:ns,myInstance) * v(1:ns,t)
forall (s = 1_pInt:ns, t = 1_pInt:4_pInt, rhoSgl(s,t+4_pInt) * v(s,t) < 0.0_pReal) & ! contribution of used rho for changing sign of v
forall (s = 1_pInt:ns, t = 1_pInt:4_pInt, rhoSgl(s,t+4_pInt) * v(s,t) < 0.0_pReal) & ! contribution of used rho for changing sign of v
gdot(s,t) = gdot(s,t) + abs(rhoSgl(s,t+4)) * constitutive_nonlocal_burgers(s,myInstance) * v(s,t)
#ifndef _OPENMP
if (debug_verbosity > 6_pInt .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then
if (debug_verbosity > 6_pInt .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) &
.or. .not. debug_selectiveDebugger)) then
write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> rho / 1/m^2', rhoSgl, rhoDip
write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> gdot / 1/s',gdot
endif

22
code/constitutive_phenopowerlaw.f90
View File

@ -345,6 +345,8 @@ subroutine constitutive_phenopowerlaw_init(file)
case ('interaction_twintwin')
forall (j = 1_pInt:lattice_maxNinteraction) &
constitutive_phenopowerlaw_interaction_twintwin(j,i) = IO_floatValue(line,positions,1_pInt+j)
case default
call IO_error(220_pInt,ext_msg=tag)
end select
endif
enddo
@ -362,21 +364,21 @@ subroutine constitutive_phenopowerlaw_init(file)
constitutive_phenopowerlaw_totalNslip(i) = sum(constitutive_phenopowerlaw_Nslip(:,i)) ! how many slip systems altogether
constitutive_phenopowerlaw_totalNtwin(i) = sum(constitutive_phenopowerlaw_Ntwin(:,i)) ! how many twin systems altogether
if (constitutive_phenopowerlaw_structure(i) < 1 ) call IO_error(205_pInt,i)
if (constitutive_phenopowerlaw_structure(i) < 1 ) call IO_error(205_pInt,e=i)
if (any(constitutive_phenopowerlaw_tau0_slip(:,i) < 0.0_pReal .and. &
constitutive_phenopowerlaw_Nslip(:,i) > 0)) call IO_error(210_pInt,i)
if (constitutive_phenopowerlaw_gdot0_slip(i) <= 0.0_pReal) call IO_error(211_pInt,i)
if (constitutive_phenopowerlaw_n_slip(i) <= 0.0_pReal) call IO_error(212_pInt,i)
constitutive_phenopowerlaw_Nslip(:,i) > 0)) call IO_error(221_pInt,e=i,ext_msg='tau0_slip')
if (constitutive_phenopowerlaw_gdot0_slip(i) <= 0.0_pReal) call IO_error(221_pInt,e=i,ext_msg='gdot0_slip')
if (constitutive_phenopowerlaw_n_slip(i) <= 0.0_pReal) call IO_error(221_pInt,e=i,ext_msg='n_slip')
if (any(constitutive_phenopowerlaw_tausat_slip(:,i) <= 0.0_pReal .and. &
constitutive_phenopowerlaw_Nslip(:,i) > 0)) call IO_error(213_pInt,i)
constitutive_phenopowerlaw_Nslip(:,i) > 0)) call IO_error(221_pInt,e=i,ext_msg='tausat_slip')
if (any(constitutive_phenopowerlaw_a_slip(i) == 0.0_pReal .and. &
constitutive_phenopowerlaw_Nslip(:,i) > 0)) call IO_error(214_pInt,i)
constitutive_phenopowerlaw_Nslip(:,i) > 0)) call IO_error(221_pInt,e=i,ext_msg='a_slip')
if (any(constitutive_phenopowerlaw_tau0_twin(:,i) < 0.0_pReal .and. &
constitutive_phenopowerlaw_Ntwin(:,i) > 0)) call IO_error(210_pInt,i)
constitutive_phenopowerlaw_Ntwin(:,i) > 0)) call IO_error(221_pInt,e=i,ext_msg='tau0_twin')
if ( constitutive_phenopowerlaw_gdot0_twin(i) <= 0.0_pReal .and. &
any(constitutive_phenopowerlaw_Ntwin(:,i) > 0)) call IO_error(211_pInt,i)
any(constitutive_phenopowerlaw_Ntwin(:,i) > 0)) call IO_error(221_pInt,e=i,ext_msg='gdot0_twin')
if ( constitutive_phenopowerlaw_n_twin(i) <= 0.0_pReal .and. &
any(constitutive_phenopowerlaw_Ntwin(:,i) > 0)) call IO_error(212_pInt,i)
any(constitutive_phenopowerlaw_Ntwin(:,i) > 0)) call IO_error(221_pInt,e=i,ext_msg='n_twin')
if (constitutive_phenopowerlaw_aTolResistance(i) <= 0.0_pReal) &
constitutive_phenopowerlaw_aTolResistance(i) = 1.0_pReal ! default absolute tolerance 1 Pa
@ -418,7 +420,7 @@ subroutine constitutive_phenopowerlaw_init(file)
)
mySize = 1_pInt
case default
mySize = 0_pInt
call IO_error(222_pInt,ext_msg=constitutive_phenopowerlaw_output(j,i))
end select
if (mySize > 0_pInt) then ! any meaningful output found

53
code/constitutive_titanmod.f90
View File

@ -583,6 +583,8 @@ do ! read thru sections of
forall (j = 1_pInt:lattice_maxNinteraction) &
constitutive_titanmod_interactionTwinTwin(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag
case default
call IO_error(230_pInt,ext_msg=tag)
end select
endif
enddo
@ -595,37 +597,38 @@ write(6,*) 'Material Property reading done'
myStructure = constitutive_titanmod_structure(i)
!* Sanity checks
if (myStructure < 1_pInt .or. myStructure > 3_pInt) call IO_error(205_pInt)
if (sum(constitutive_titanmod_Nslip(:,i)) <= 0_pInt) call IO_error(207_pInt)
if (sum(constitutive_titanmod_Ntwin(:,i)) < 0_pInt) call IO_error(208_pInt) !***
if (myStructure < 1_pInt .or. myStructure > 3_pInt) call IO_error(205_pInt,e=i)
if (sum(constitutive_titanmod_Nslip(:,i)) <= 0_pInt) call IO_error(231_pInt,e=i,ext_msg='nslip')
if (sum(constitutive_titanmod_Ntwin(:,i)) < 0_pInt) call IO_error(231_pInt,e=i,ext_msg='ntwin')
do f = 1,lattice_maxNslipFamily
if (constitutive_titanmod_Nslip(f,i) > 0_pInt) then
if (constitutive_titanmod_rho_edge0(f,i) < 0.0_pReal) call IO_error(209_pInt)
if (constitutive_titanmod_rho_screw0(f,i) < 0.0_pReal) call IO_error(210_pInt)
if (constitutive_titanmod_burgersPerSlipFamily(f,i) <= 0.0_pReal) call IO_error(211_pInt)
if (constitutive_titanmod_f0_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(212_pInt)
if (constitutive_titanmod_tau0e_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(229_pInt)
if (constitutive_titanmod_tau0s_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(233_pInt)
if (constitutive_titanmod_capre_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(234_pInt)
if (constitutive_titanmod_caprs_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(235_pInt)
if (constitutive_titanmod_v0e_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(226_pInt)
if (constitutive_titanmod_v0s_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(226_pInt)
if (constitutive_titanmod_kinkcriticallength_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(238_pInt)
if (constitutive_titanmod_rho_edge0(f,i) < 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='rho_edge0')
if (constitutive_titanmod_rho_screw0(f,i) < 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='rho_screw0')
if (constitutive_titanmod_burgersPerSlipFamily(f,i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='slipburgers')
if (constitutive_titanmod_f0_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='f0')
if (constitutive_titanmod_tau0e_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='tau0e')
if (constitutive_titanmod_tau0s_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='tau0s')
if (constitutive_titanmod_capre_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='capre')
if (constitutive_titanmod_caprs_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='carrs')
if (constitutive_titanmod_v0e_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='v0e')
if (constitutive_titanmod_v0s_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='v0s')
if (constitutive_titanmod_kinkcriticallength_PerSlipFamily(f,i) <= 0.0_pReal) &
call IO_error(231_pInt,e=i,ext_msg='kinkCriticalLength')
endif
enddo
do f = 1,lattice_maxNtwinFamily
if (constitutive_titanmod_Ntwin(f,i) > 0_pInt) then
if (constitutive_titanmod_burgersPerTwinFamily(f,i) <= 0.0_pReal) call IO_error(221_pInt) !***
if (constitutive_titanmod_twinf0_PerTwinFamily(f,i) <= 0.0_pReal) call IO_error(228_pInt)
if (constitutive_titanmod_twinshearconstant_PerTwinFamily(f,i) <= 0.0_pReal) call IO_error(228_pInt)
if (constitutive_titanmod_twintau0_PerTwinFamily(f,i) <= 0.0_pReal) call IO_error(229_pInt)
if (constitutive_titanmod_twingamma0_PerTwinFamily(f,i) <= 0.0_pReal) call IO_error(226_pInt)
if (constitutive_titanmod_burgersPerTwinFamily(f,i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='twinburgers')
if (constitutive_titanmod_twinf0_PerTwinFamily(f,i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='twinf0')
if (constitutive_titanmod_twinshearconstant_PerTwinFamily(f,i) <= 0.0_pReal) &
call IO_error(231_pInt,e=i,ext_msg='twinshearconstant')
if (constitutive_titanmod_twintau0_PerTwinFamily(f,i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='twintau0')
if (constitutive_titanmod_twingamma0_PerTwinFamily(f,i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='twingamma0')
endif
enddo
! if (any(constitutive_titanmod_interactionSlipSlip(1:maxval(lattice_interactionSlipSlip(:,:,myStructure)),i) < 1.0_pReal)) call IO_error(229)
if (constitutive_titanmod_dc(i) <= 0.0_pReal) call IO_error(231_pInt)
if (constitutive_titanmod_twinhpconstant(i) <= 0.0_pReal) call IO_error(232_pInt)
if (constitutive_titanmod_aTolRho(i) <= 0.0_pReal) call IO_error(233_pInt)
if (constitutive_titanmod_dc(i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='dc')
if (constitutive_titanmod_twinhpconstant(i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='twinhpconstant')
if (constitutive_titanmod_aTolRho(i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='aTolRho')
!* Determine total number of active slip or twin systems
constitutive_titanmod_Nslip(:,i) = min(lattice_NslipSystem(:,myStructure),constitutive_titanmod_Nslip(:,i))
@ -637,7 +640,7 @@ enddo
!* Allocation of variables whose size depends on the total number of active slip systems
maxTotalNslip = maxval(constitutive_titanmod_totalNslip)
maxTotalNtwin = maxval(constitutive_titanmod_totalNtwin)
maxTotalNtwin = maxval(constitutive_titanmod_totalNtwin)
write(6,*) 'maxTotalNslip',maxTotalNslip
write(6,*) 'maxTotalNtwin',maxTotalNtwin
allocate(constitutive_titanmod_burgersPerSlipSystem(maxTotalNslip, maxNinstance))
@ -796,7 +799,7 @@ do i = 1_pInt,maxNinstance
)
mySize = 1_pInt
case default
mySize = 0_pInt
call IO_error(232_pInt,ext_msg=constitutive_titanmod_output(o,i))
end select
if (mySize > 0_pInt) then ! any meaningful output found

6
code/crystallite.f90
View File

@ -234,8 +234,8 @@ allocate(crystallite_sizePostResult(maxval(crystallite_Noutput), &
material_Ncrystallite)) ; crystallite_sizePostResult = 0_pInt
if (.not. IO_open_jobFile(file,material_localFileExt)) then ! no local material configuration present...
if (.not. IO_open_file(file,material_configFile)) call IO_error(100_pInt) ! ...and cannot open material.config file
if (.not. IO_open_jobFile_stat(file,material_localFileExt)) then ! no local material configuration present...
call IO_open_file(file,material_configFile) ! ...open material.config file
endif
line = ''
section = 0_pInt
@ -299,7 +299,7 @@ enddo
! write description file for crystallite output
if(.not. IO_write_jobFile(file,'outputCrystallite')) call IO_error (50_pInt) ! problems in writing file
call IO_write_jobFile(file,'outputCrystallite')
do p = 1_pInt,material_Ncrystallite
write(file,*)

4
code/debug.f90
View File

@ -73,7 +73,7 @@ subroutine debug_init()
nMPstate, &
nHomog
use IO, only: IO_error, &
IO_open_file, &
IO_open_file_stat, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
@ -108,7 +108,7 @@ subroutine debug_init()
allocate(debug_MaterialpointLoopDistribution(nHomog+1)) ; debug_MaterialpointLoopDistribution = 0_pInt
! try to open the config file
if(IO_open_file(fileunit,debug_configFile)) then
if(IO_open_file_stat(fileunit,debug_configFile)) then
line = ''
! read variables from config file and overwrite parameters

12
code/homogenization.f90
View File

@ -75,7 +75,7 @@ subroutine homogenization_init(Temperature)
use prec, only: pReal,pInt
use math, only: math_I3
use debug, only: debug_verbosity
use IO, only: IO_error, IO_open_file, IO_open_jobFile, IO_write_jobFile
use IO, only: IO_error, IO_open_file, IO_open_jobFile_stat, IO_write_jobFile
use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips
use material
use constitutive, only: constitutive_maxSizePostResults
@ -95,8 +95,8 @@ logical knownHomogenization
! --- PARSE HOMOGENIZATIONS FROM CONFIG FILE ---
if (.not. IO_open_jobFile(fileunit,material_localFileExt)) then ! no local material configuration present...
if (.not. IO_open_file(fileunit,material_configFile)) call IO_error(100) ! ...and cannot open material.config file
if (.not. IO_open_jobFile_stat(fileunit,material_localFileExt)) then ! no local material configuration present...
call IO_open_file(fileunit,material_configFile) ! ... open material.config file
endif
call homogenization_isostrain_init(fileunit)
call homogenization_RGC_init(fileunit)
@ -105,9 +105,7 @@ close(fileunit)
! --- WRITE DESCRIPTION FILE FOR HOMOGENIZATION OUTPUT ---
if(.not. IO_write_jobFile(fileunit,'outputHomogenization')) then ! problems in writing file
call IO_error (50)
endif
call IO_write_jobFile(fileunit,'outputHomogenization')
do p = 1,material_Nhomogenization
i = homogenization_typeInstance(p) ! which instance of this homogenization type
knownHomogenization = .true. ! assume valid
@ -189,7 +187,7 @@ end forall
endif
homogenization_sizePostResults(i,e) = homogenization_RGC_sizePostResults(myInstance)
case default
call IO_error(201,ext_msg=homogenization_type(mesh_element(3,e))) ! unknown type 201 is homogenization!
call IO_error(500_pInt,ext_msg=homogenization_type(mesh_element(3,e))) ! unknown homogenization
end select
enddo
enddo

8
code/lattice.f90
View File

@ -728,7 +728,7 @@ subroutine lattice_init()
!* Module initialization *
!**************************************
use, intrinsic :: iso_fortran_env
use IO, only: IO_open_file,IO_open_jobFile,IO_countSections,IO_countTagInPart,IO_error
use IO, only: IO_open_file,IO_open_jobFile_stat,IO_countSections,IO_countTagInPart,IO_error
use material, only: material_configfile,material_localFileExt,material_partPhase
use debug, only: debug_verbosity
implicit none
@ -743,8 +743,8 @@ subroutine lattice_init()
#include "compilation_info.f90"
!$OMP END CRITICAL (write2out)
if (.not. IO_open_jobFile(fileunit,material_localFileExt)) then ! no local material configuration present...
if (.not. IO_open_file(fileunit,material_configFile)) call IO_error(100_pInt) ! ...and cannot open material.config file
if (.not. IO_open_jobFile_stat(fileunit,material_localFileExt)) then ! no local material configuration present...
call IO_open_file(fileunit,material_configFile) ! ... open material.config file
endif
Nsections = IO_countSections(fileunit,material_partPhase)
lattice_Nstructure = 2_pInt + sum(IO_countTagInPart(fileunit,material_partPhase,'covera_ratio',Nsections)) ! fcc + bcc + all hex
@ -922,7 +922,7 @@ function lattice_initializeStructure(struct,CoverA)
if (processMe) then
if (myStructure > lattice_Nstructure) &
call IO_error(666_pint,0_pInt,0_pInt,0_pInt,'structure index too large') ! check for memory leakage
call IO_error(666_pInt,0_pInt,0_pInt,0_pInt,'structure index too large') ! check for memory leakage
do i = 1,myNslip ! store slip system vectors and Schmid matrix for my structure
lattice_sd(1:3,i,myStructure) = sd(1:3,i)
lattice_st(1:3,i,myStructure) = st(1:3,i)

28
code/material.f90
View File

@ -109,7 +109,7 @@ subroutine material_init()
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: pReal,pInt
use IO, only: IO_error, IO_open_file, IO_open_jobFile
use IO, only: IO_error, IO_open_file, IO_open_jobFile_stat
use debug, only: debug_verbosity
implicit none
@ -124,8 +124,8 @@ subroutine material_init()
#include "compilation_info.f90"
!$OMP END CRITICAL (write2out)
if (.not. IO_open_jobFile(fileunit,material_localFileExt)) then ! no local material configuration present...
if (.not. IO_open_file(fileunit,material_configFile)) call IO_error(100_pInt) ! ...and cannot open material.config file
if (.not. IO_open_jobFile_stat(fileunit,material_localFileExt)) then ! no local material configuration present...
call IO_open_file(fileunit,material_configFile) ! ...open material.config file
endif
call material_parseHomogenization(fileunit,material_partHomogenization)
if (debug_verbosity > 0) then
@ -163,16 +163,16 @@ subroutine material_init()
if (microstructure_crystallite(i) < 1 .or. &
microstructure_crystallite(i) > material_Ncrystallite) call IO_error(150_pInt,i)
if (minval(microstructure_phase(1:microstructure_Nconstituents(i),i)) < 1 .or. &
maxval(microstructure_phase(1:microstructure_Nconstituents(i),i)) > material_Nphase) call IO_error(155_pInt,i)
maxval(microstructure_phase(1:microstructure_Nconstituents(i),i)) > material_Nphase) call IO_error(151_pInt,i)
if (minval(microstructure_texture(1:microstructure_Nconstituents(i),i)) < 1 .or. &
maxval(microstructure_texture(1:microstructure_Nconstituents(i),i)) > material_Ntexture) call IO_error(160_pInt,i)
maxval(microstructure_texture(1:microstructure_Nconstituents(i),i)) > material_Ntexture) call IO_error(152_pInt,i)
if (abs(sum(microstructure_fraction(:,i)) - 1.0_pReal) >= 1.0e-10_pReal) then
if (debug_verbosity > 0) then
!$OMP CRITICAL (write2out)
write(6,*)'sum of microstructure fraction = ',sum(microstructure_fraction(:,i))
!$OMP END CRITICAL (write2out)
endif
call IO_error(170_pInt,i)
call IO_error(153_pInt,i)
endif
enddo
if (debug_verbosity > 0) then
@ -227,7 +227,7 @@ subroutine material_parseHomogenization(file,myPart)
Nsections = IO_countSections(file,myPart)
material_Nhomogenization = Nsections
if (Nsections < 1_pInt) call IO_error(125_pInt,ext_msg=myPart)
if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart)
allocate(homogenization_name(Nsections)); homogenization_name = ''
allocate(homogenization_type(Nsections)); homogenization_type = ''
@ -295,7 +295,7 @@ subroutine material_parseMicrostructure(file,myPart)
Nsections = IO_countSections(file,myPart)
material_Nmicrostructure = Nsections
if (Nsections < 1_pInt) call IO_error(125_pInt,ext_msg=myPart)
if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart)
allocate(microstructure_name(Nsections)); microstructure_name = ''
allocate(microstructure_crystallite(Nsections)); microstructure_crystallite = 0_pInt
@ -372,7 +372,7 @@ subroutine material_parseCrystallite(file,myPart)
Nsections = IO_countSections(file,myPart)
material_Ncrystallite = Nsections
if (Nsections < 1_pInt) call IO_error(125_pInt,ext_msg=myPart)
if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart)
allocate(crystallite_name(Nsections)); crystallite_name = ''
allocate(crystallite_Noutput(Nsections)); crystallite_Noutput = 0_pInt
@ -418,7 +418,7 @@ subroutine material_parsePhase(file,myPart)
Nsections = IO_countSections(file,myPart)
material_Nphase = Nsections
if (Nsections < 1_pInt) call IO_error(125_pInt,ext_msg=myPart)
if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart)
allocate(phase_name(Nsections)); phase_name = ''
allocate(phase_constitution(Nsections)); phase_constitution = ''
@ -482,7 +482,7 @@ subroutine material_parseTexture(file,myPart)
Nsections = IO_countSections(file,myPart)
material_Ntexture = Nsections
if (Nsections < 1_pInt) call IO_error(125_pInt,ext_msg=myPart)
if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart)
allocate(texture_name(Nsections)); texture_name = ''
allocate(texture_ODFfile(Nsections)); texture_ODFfile = ''
@ -645,9 +645,9 @@ subroutine material_populateGrains()
homog = mesh_element(3,e)
micro = mesh_element(4,e)
if (homog < 1 .or. homog > material_Nhomogenization) & ! out of bounds
call IO_error(130_pInt,e,0_pInt,0_pInt)
call IO_error(154_pInt,e,0_pInt,0_pInt)
if (micro < 1 .or. micro > material_Nmicrostructure) & ! out of bounds
call IO_error(140_pInt,e,0_pInt,0_pInt)
call IO_error(155_pInt,e,0_pInt,0_pInt)
if (microstructure_elemhomo(micro)) then
dGrains = homogenization_Ngrains(homog)
else
@ -762,7 +762,7 @@ subroutine material_populateGrains()
else ! hybrid IA
! ---------
orientationOfGrain(:,grain+1:grain+myNorientations) = IO_hybridIA(myNorientations,texture_ODFfile(textureID))
if (all(orientationOfGrain(:,grain+1) == -1.0_pReal)) call IO_error(105_pInt)
if (all(orientationOfGrain(:,grain+1) == -1.0_pReal)) call IO_error(156_pInt)
constituentGrain = constituentGrain + myNorientations
endif

32
code/math.f90
View File

@ -201,7 +201,7 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = &
if ( any(abs( q-q2) > tol_math_check) .and. &
any(abs(-q-q2) > tol_math_check) ) then
write (error_msg, '(a,e14.6)' ) 'maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2)))
call IO_error(670_pInt,ext_msg=error_msg)
call IO_error(401_pInt,ext_msg=error_msg)
endif
! +++ q -> R -> q +++
@ -210,7 +210,7 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = &
if ( any(abs( q-q2) > tol_math_check) .and. &
any(abs(-q-q2) > tol_math_check) ) then
write (error_msg, '(a,e14.6)' ) 'maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2)))
call IO_error(671_pInt,ext_msg=error_msg)
call IO_error(402_pInt,ext_msg=error_msg)
endif
! +++ q -> euler -> q +++
@ -219,7 +219,7 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = &
if ( any(abs( q-q2) > tol_math_check) .and. &
any(abs(-q-q2) > tol_math_check) ) then
write (error_msg, '(a,e14.6)' ) 'maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2)))
call IO_error(672_pInt,ext_msg=error_msg)
call IO_error(403_pInt,ext_msg=error_msg)
endif
! +++ R -> euler -> R +++
@ -227,7 +227,7 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = &
R2 = math_EulerToR(Eulers)
if ( any(abs( R-R2) > tol_math_check) ) then
write (error_msg, '(a,e14.6)' ) 'maximum deviation ',maxval(abs( R-R2))
call IO_error(673_pInt,ext_msg=error_msg)
call IO_error(404_pInt,ext_msg=error_msg)
endif
ENDSUBROUTINE math_init
@ -1820,7 +1820,7 @@ function math_QuaternionDisorientation(Q1, Q2, symmetryType)
enddo; enddo; enddo
case default
call IO_error(550_pInt,symmetryType) ! complain about unknown symmetry
call IO_error(450_pInt,symmetryType) ! complain about unknown symmetry
end select
endfunction math_QuaternionDisorientation
@ -2553,7 +2553,7 @@ end subroutine
r(1:ndim) = 0.0_pReal
if (any (base(1:ndim) <= 1_pInt)) call IO_error(error_ID=801_pInt)
if (any (base(1:ndim) <= 1_pInt)) call IO_error(error_ID=405_pInt)
base_inv(1:ndim) = 1.0_pReal / real (base(1:ndim), pReal)
@ -2795,7 +2795,7 @@ end subroutine
else if (n <= prime_max) then
prime = npvec(n)
else
call IO_error(error_ID=802_pInt)
call IO_error(error_ID=406_pInt)
end if
endfunction prime
@ -3262,7 +3262,7 @@ subroutine deformed_fft(res,geomdim,defgrad_av,scaling,defgrad,coords)
res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c)
step = geomdim/real(res, pReal)
if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=102_pInt)
if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=808_pInt)
call fftw_set_timelimit(fftw_timelimit)
defgrad_fftw = fftw_alloc_complex(int(res1_red *res(2)*res(3)*9_pInt,C_SIZE_T)) !C_SIZE_T is of type integer(8)
call c_f_pointer(defgrad_fftw, defgrad_real, [res(1)+2_pInt,res(2),res(3),3_pInt,3_pInt])
@ -3385,7 +3385,7 @@ subroutine curl_fft(res,geomdim,vec_tens,field,curl)
wgt = 1.0_pReal/real(res(1)*res(2)*res(3),pReal)
res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c)
if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=102_pInt)
if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=808_pInt)
call fftw_set_timelimit(fftw_timelimit)
field_fftw = fftw_alloc_complex(int(res1_red *res(2)*res(3)*vec_tens*3_pInt,C_SIZE_T)) !C_SIZE_T is of type integer(8)
call c_f_pointer(field_fftw, field_real, [res(1)+2_pInt,res(2),res(3),vec_tens,3_pInt])
@ -3500,10 +3500,10 @@ subroutine divergence_fft(res,geomdim,vec_tens,field,divergence)
print '(a,3(i5))', ' Resolution:', res
endif
res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c)
res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c)
wgt = 1.0_pReal/real(res(1)*res(2)*res(3),pReal)
if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=102_pInt)
if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=808_pInt)
call fftw_set_timelimit(fftw_timelimit)
field_fftw = fftw_alloc_complex(int(res1_red*res(2)*res(3)*vec_tens*3_pInt,C_SIZE_T)) !C_SIZE_T is of type integer(8)
call c_f_pointer(field_fftw, field_real, [res(1)+2_pInt,res(2),res(3),vec_tens,3_pInt])
@ -3512,9 +3512,9 @@ if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=102_pInt)
call c_f_pointer(divergence_fftw, divergence_real, [res(1)+2_pInt,res(2),res(3),vec_tens])
call c_f_pointer(divergence_fftw, divergence_fourier,[res1_red ,res(2),res(3),vec_tens])
fftw_forth = fftw_plan_many_dft_r2c(3_pInt,(/res(3),res(2) ,res(1)/),vec_tens*3_pInt,& ! dimensions , length in each dimension in reversed order
fftw_forth = fftw_plan_many_dft_r2c(3_pInt,(/res(3),res(2) ,res(1)/),vec_tens*3_pInt,& ! dimensions , length in each dimension in reversed order
field_real,(/res(3),res(2) ,res(1)+2_pInt/),& ! input data , physical length in each dimension in reversed order
1_pInt, res(3)*res(2)*(res(1)+2_pInt),& ! striding , product of physical lenght in the 3 dimensions
1_pInt, res(3)*res(2)*(res(1)+2_pInt),& ! striding , product of physical lenght in the 3 dimensions
field_fourier,(/res(3),res(2) ,res1_red/),&
1_pInt, res(3)*res(2)* res1_red,fftw_planner_flag)
@ -3522,7 +3522,7 @@ if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=102_pInt)
divergence_fourier,(/res(3),res(2) ,res1_red/),&
1_pInt, res(3)*res(2)* res1_red,&
divergence_real,(/res(3),res(2) ,res(1)+2_pInt/),&
1_pInt, res(3)*res(2)*(res(1)+2_pInt),fftw_planner_flag) ! padding
1_pInt, res(3)*res(2)*(res(1)+2_pInt),fftw_planner_flag) ! padding
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
field_real(i,j,k,1:vec_tens,1:3) = field(i,j,k,1:vec_tens,1:3) ! ensure that data is aligned properly (fftw_alloc)
enddo; enddo; enddo
@ -3559,12 +3559,12 @@ if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=102_pInt)
call fftw_execute_dft_c2r(fftw_back, divergence_fourier, divergence_real)
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
divergence(i,j,k,1:vec_tens) = divergence_real(i,j,k,1:vec_tens) ! ensure that data is aligned properly (fftw_alloc)
divergence(i,j,k,1:vec_tens) = divergence_real(i,j,k,1:vec_tens) ! ensure that data is aligned properly (fftw_alloc)
enddo; enddo; enddo
divergence = divergence * wgt
call fftw_destroy_plan(fftw_forth); call fftw_destroy_plan(fftw_back)
call c_f_pointer(C_NULL_PTR, field_real, [res(1)+2_pInt,res(2),res(3),vec_tens,3_pInt]) ! let all pointers point on NULL-Type
call c_f_pointer(C_NULL_PTR, field_real, [res(1)+2_pInt,res(2),res(3),vec_tens,3_pInt]) ! let all pointers point on NULL-Type
call c_f_pointer(C_NULL_PTR, field_fourier, [res1_red ,res(2),res(3),vec_tens,3_pInt])
call c_f_pointer(C_NULL_PTR, divergence_real, [res(1)+2_pInt,res(2),res(3),vec_tens])
call c_f_pointer(C_NULL_PTR, divergence_fourier,[res1_red ,res(2),res(3),vec_tens])

121
code/mesh.f90
View File

@ -281,59 +281,56 @@
call mesh_build_FEdata() ! --- get properties of the different types of elements
if (IO_open_inputFile(fileUnit,FEmodelGeometry)) then ! --- parse info from input file...
call IO_open_inputFile(fileUnit,FEmodelGeometry) ! --- parse info from input file...
select case (FEsolver)
case ('Spectral')
call mesh_spectral_count_nodesAndElements(fileUnit)
call mesh_spectral_count_cpElements()
call mesh_spectral_map_elements()
call mesh_spectral_map_nodes()
call mesh_spectral_count_cpSizes()
call mesh_spectral_build_nodes(fileUnit)
call mesh_spectral_build_elements(fileUnit)
case ('Marc')
call mesh_marc_get_tableStyles(fileUnit)
call mesh_marc_count_nodesAndElements(fileUnit)
call mesh_marc_count_elementSets(fileUnit)
call mesh_marc_map_elementSets(fileUnit)
call mesh_marc_count_cpElements(fileUnit)
call mesh_marc_map_elements(fileUnit)
call mesh_marc_map_nodes(fileUnit)
call mesh_marc_build_nodes(fileUnit)
call mesh_marc_count_cpSizes(fileunit)
call mesh_marc_build_elements(fileUnit)
case ('Abaqus')
noPart = IO_abaqus_hasNoPart(fileUnit)
call mesh_abaqus_count_nodesAndElements(fileUnit)
call mesh_abaqus_count_elementSets(fileUnit)
call mesh_abaqus_count_materials(fileUnit)
call mesh_abaqus_map_elementSets(fileUnit)
call mesh_abaqus_map_materials(fileUnit)
call mesh_abaqus_count_cpElements(fileUnit)
call mesh_abaqus_map_elements(fileUnit)
call mesh_abaqus_map_nodes(fileUnit)
call mesh_abaqus_build_nodes(fileUnit)
call mesh_abaqus_count_cpSizes(fileunit)
call mesh_abaqus_build_elements(fileUnit)
end select
call mesh_get_damaskOptions(fileUnit)
close (fileUnit)
select case (FEsolver)
case ('Spectral')
call mesh_spectral_count_nodesAndElements(fileUnit)
call mesh_spectral_count_cpElements()
call mesh_spectral_map_elements()
call mesh_spectral_map_nodes()
call mesh_spectral_count_cpSizes()
call mesh_spectral_build_nodes(fileUnit)
call mesh_spectral_build_elements(fileUnit)
call mesh_build_subNodeCoords()
call mesh_build_ipCoordinates()
call mesh_build_ipVolumes()
call mesh_build_ipAreas()
call mesh_build_nodeTwins()
call mesh_build_sharedElems()
call mesh_build_ipNeighborhood()
call mesh_tell_statistics()
case ('Marc')
call mesh_marc_get_tableStyles(fileUnit)
call mesh_marc_count_nodesAndElements(fileUnit)
call mesh_marc_count_elementSets(fileUnit)
call mesh_marc_map_elementSets(fileUnit)
call mesh_marc_count_cpElements(fileUnit)
call mesh_marc_map_elements(fileUnit)
call mesh_marc_map_nodes(fileUnit)
call mesh_marc_build_nodes(fileUnit)
call mesh_marc_count_cpSizes(fileunit)
call mesh_marc_build_elements(fileUnit)
case ('Abaqus')
noPart = IO_abaqus_hasNoPart(fileUnit)
call mesh_abaqus_count_nodesAndElements(fileUnit)
call mesh_abaqus_count_elementSets(fileUnit)
call mesh_abaqus_count_materials(fileUnit)
call mesh_abaqus_map_elementSets(fileUnit)
call mesh_abaqus_map_materials(fileUnit)
call mesh_abaqus_count_cpElements(fileUnit)
call mesh_abaqus_map_elements(fileUnit)
call mesh_abaqus_map_nodes(fileUnit)
call mesh_abaqus_build_nodes(fileUnit)
call mesh_abaqus_count_cpSizes(fileunit)
call mesh_abaqus_build_elements(fileUnit)
end select
call mesh_get_damaskOptions(fileUnit)
close (fileUnit)
call mesh_build_subNodeCoords()
call mesh_build_ipCoordinates()
call mesh_build_ipVolumes()
call mesh_build_ipAreas()
call mesh_build_nodeTwins()
call mesh_build_sharedElems()
call mesh_build_ipNeighborhood()
call mesh_tell_statistics()
parallelExecution = (parallelExecution .and. (mesh_Nelems == mesh_NcpElems)) ! plus potential killer from non-local constitutive
else
call IO_error(error_ID=101_pInt) ! cannot open input file
endif
parallelExecution = (parallelExecution .and. (mesh_Nelems == mesh_NcpElems)) ! plus potential killer from non-local constitutive
FEsolving_execElem = [ 1_pInt,mesh_NcpElems]
allocate(FEsolving_execIP(2,mesh_NcpElems)); FEsolving_execIP = 1_pInt
@ -1487,7 +1484,7 @@ enddo
if (keyword(1:4) == 'head') then
headerLength = IO_intValue(line,myPos,1_pInt) + 1_pInt
else
call IO_error(error_ID=42_pInt)
call IO_error(error_ID=842_pInt)
endif
rewind(myUnit)
@ -2428,7 +2425,7 @@ subroutine mesh_marc_count_cpSizes (myUnit)
if (keyword(1:4) == 'head') then
headerLength = IO_intValue(line,myPos,1_pInt) + 1_pInt
else
call IO_error(error_ID=42_pInt)
call IO_error(error_ID=842_pInt)
endif
rewind(myUnit)
@ -2465,9 +2462,9 @@ subroutine mesh_marc_count_cpSizes (myUnit)
! --- sanity checks ---
if ((.not. gotDimension) .or. (.not. gotResolution)) call IO_error(error_ID=42_pInt)
if ((a < 1) .or. (b < 1) .or. (c < 0)) call IO_error(error_ID=43_pInt) ! 1_pInt is already added
if ((x <= 0.0_pReal) .or. (y <= 0.0_pReal) .or. (z <= 0.0_pReal)) call IO_error(error_ID=44_pInt)
if ((.not. gotDimension) .or. (.not. gotResolution)) call IO_error(error_ID=842_pInt)
if ((a < 1) .or. (b < 1) .or. (c < 0)) call IO_error(error_ID=843_pInt) ! 1_pInt is already added
if ((x <= 0.0_pReal) .or. (y <= 0.0_pReal) .or. (z <= 0.0_pReal)) call IO_error(error_ID=844_pInt)
forall (n = 0:mesh_Nnodes-1)
mesh_node0(1,n+1) = x * dble(mod(n,a)) / dble(a-1_pInt)
@ -2616,7 +2613,7 @@ subroutine mesh_marc_count_cpSizes (myUnit)
if (keyword(1:4) == 'head') then
headerLength = IO_intValue(line,myPos,1_pInt) + 1_pInt
else
call IO_error(error_ID=42_pInt)
call IO_error(error_ID=842_pInt)
endif
rewind(myUnit)
@ -2679,7 +2676,7 @@ subroutine mesh_marc_count_cpSizes (myUnit)
enddo
110 deallocate(microstructures)
if (e /= mesh_NcpElems) call IO_error(180_pInt,e)
if (e /= mesh_NcpElems) call IO_error(880_pInt,e)
endsubroutine
@ -3365,13 +3362,13 @@ character(len=64) fmt
integer(pInt) i,e,n,f,t
if (mesh_maxValStateVar(1) < 1_pInt) call IO_error(error_ID=110_pInt) ! no homogenization specified
if (mesh_maxValStateVar(2) < 1_pInt) call IO_error(error_ID=120_pInt) ! no microstructure specified
if (mesh_maxValStateVar(1) < 1_pInt) call IO_error(error_ID=170_pInt) ! no homogenization specified
if (mesh_maxValStateVar(2) < 1_pInt) call IO_error(error_ID=180_pInt) ! no microstructure specified
allocate (mesh_HomogMicro(mesh_maxValStateVar(1),mesh_maxValStateVar(2))); mesh_HomogMicro = 0_pInt
do e = 1,mesh_NcpElems
if (mesh_element(3,e) < 1_pInt) call IO_error(error_ID=110_pInt,e=e) ! no homogenization specified
if (mesh_element(4,e) < 1_pInt) call IO_error(error_ID=120_pInt,e=e) ! no microstructure specified
if (mesh_element(3,e) < 1_pInt) call IO_error(error_ID=170_pInt,e=e) ! no homogenization specified
if (mesh_element(4,e) < 1_pInt) call IO_error(error_ID=180_pInt,e=e) ! no microstructure specified
mesh_HomogMicro(mesh_element(3,e),mesh_element(4,e)) = &
mesh_HomogMicro(mesh_element(3,e),mesh_element(4,e)) + 1 ! count combinations of homogenization and microstructure
enddo
@ -3490,4 +3487,4 @@ deallocate(mesh_HomogMicro)
endsubroutine
END MODULE mesh
END MODULE mesh

124
code/numerics.f90
View File

@ -97,7 +97,7 @@ subroutine numerics_init()
use prec, only: pInt, &
pReal
use IO, only: IO_error, &
IO_open_file, &
IO_open_file_stat, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
@ -133,15 +133,17 @@ subroutine numerics_init()
!$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! ...and use it as number of threads for parallel execution
! try to open the config file
if(IO_open_file(fileunit,numerics_configFile)) then
if(IO_open_file_stat(fileunit,numerics_configFile)) then
!$OMP CRITICAL (write2out)
write(6,*) ' ... using values from config file'
write(6,*)
!$OMP END CRITICAL (write2out)
!* read variables from config file and overwrite parameters
line = ''
! read variables from config file and overwrite parameters
do
read(fileunit,'(a1024)',END=100) line
if (IO_isBlank(line)) cycle ! skip empty lines
@ -195,7 +197,8 @@ subroutine numerics_init()
case ('integratorstiffness')
numerics_integrator(2) = IO_intValue(line,positions,2_pInt)
!* RGC parameters:
!* RGC parameters:
case ('atol_rgc')
absTol_RGC = IO_floatValue(line,positions,2_pInt)
case ('rtol_rgc')
@ -223,7 +226,8 @@ subroutine numerics_init()
case ('discrepancypower_rgc')
volDiscrPow_RGC = IO_floatValue(line,positions,2_pInt)
!* spectral parameters
!* spectral parameters
case ('err_div_tol')
err_div_tol = IO_floatValue(line,positions,2_pInt)
case ('err_stress_tolrel')
@ -247,9 +251,13 @@ subroutine numerics_init()
case ('cut_off_value')
cut_off_value = IO_floatValue(line,positions,2_pInt)
!* Random seeding parameters
!* Random seeding parameters
case ('fixed_seed')
fixedSeed = IO_intValue(line,positions,2_pInt)
case default
call IO_error(300_pInt,ext_msg=tag)
endselect
enddo
100 close(fileunit)
@ -263,6 +271,7 @@ subroutine numerics_init()
!$OMP END CRITICAL (write2out)
endif
select case(IO_lc(fftw_planner_string)) ! setting parameters for the plan creation of FFTW. Basically a translation from fftw3.f
case('estimate','fftw_estimate') ! ordered from slow execution (but fast plan creation) to fast execution
fftw_planner_flag = 64_pInt
@ -277,8 +286,11 @@ subroutine numerics_init()
fftw_planner_flag = 32_pInt
end select
! writing parameters to output file
!* writing parameters to output file
!$OMP CRITICAL (write2out)
write(6,'(a24,1x,e8.1)') ' relevantStrain: ',relevantStrain
write(6,'(a24,1x,e8.1)') ' defgradTolerance: ',defgradTolerance
write(6,'(a24,1x,i8)') ' iJacoStiffness: ',iJacoStiffness
@ -303,7 +315,8 @@ subroutine numerics_init()
write(6,'(a24,1x,e8.1)') ' stepIncreaseHomog: ',stepIncreaseHomog
write(6,'(a24,1x,i8,/)') ' nMPstate: ',nMPstate
!* RGC parameters
!* RGC parameters
write(6,'(a24,1x,e8.1)') ' aTol_RGC: ',absTol_RGC
write(6,'(a24,1x,e8.1)') ' rTol_RGC: ',relTol_RGC
write(6,'(a24,1x,e8.1)') ' aMax_RGC: ',absMax_RGC
@ -317,7 +330,8 @@ subroutine numerics_init()
write(6,'(a24,1x,e8.1)') ' volDiscrepancyMod_RGC: ',volDiscrMod_RGC
write(6,'(a24,1x,e8.1,/)') ' discrepancyPower_RGC: ',volDiscrPow_RGC
!* spectral parameters
!* spectral parameters
write(6,'(a24,1x,e8.1)') ' err_div_tol: ',err_div_tol
write(6,'(a24,1x,e8.1)') ' err_stress_tolrel: ',err_stress_tolrel
write(6,'(a24,1x,i8)') ' itmax: ',itmax
@ -334,58 +348,66 @@ subroutine numerics_init()
write(6,'(a24,1x,L8,/)') ' update_gamma: ',update_gamma
write(6,'(a24,1x,L8,/)') ' simplified_algorithm: ',simplified_algorithm
write(6,'(a24,1x,e8.1)') ' cut_off_value: ',cut_off_value
!* Random seeding parameters
!* Random seeding parameters
write(6,'(a24,1x,i16,/)') ' fixed_seed: ',fixedSeed
!$OMP END CRITICAL (write2out)
!* openMP parameter
!* openMP parameter
!$ write(6,'(a24,1x,i8,/)') ' number of threads: ',DAMASK_NumThreadsInt
! sanity check
if (relevantStrain <= 0.0_pReal) call IO_error(260_pInt)
if (defgradTolerance <= 0.0_pReal) call IO_error(294_pInt)
if (iJacoStiffness < 1_pInt) call IO_error(261_pInt)
if (iJacoLpresiduum < 1_pInt) call IO_error(262_pInt)
if (pert_Fg <= 0.0_pReal) call IO_error(263_pInt)
!* sanity check
if (relevantStrain <= 0.0_pReal) call IO_error(301_pInt,ext_msg='relevantStrain')
if (defgradTolerance <= 0.0_pReal) call IO_error(301_pInt,ext_msg='defgradTolerance')
if (iJacoStiffness < 1_pInt) call IO_error(301_pInt,ext_msg='iJacoStiffness')
if (iJacoLpresiduum < 1_pInt) call IO_error(301_pInt,ext_msg='iJacoLpresiduum')
if (pert_Fg <= 0.0_pReal) call IO_error(301_pInt,ext_msg='pert_Fg')
if (pert_method <= 0_pInt .or. pert_method >= 4_pInt) &
call IO_error(299_pInt)
if (nHomog < 1_pInt) call IO_error(264_pInt)
if (nMPstate < 1_pInt) call IO_error(279_pInt) !! missing in IO !!
if (nCryst < 1_pInt) call IO_error(265_pInt)
if (nState < 1_pInt) call IO_error(266_pInt)
if (nStress < 1_pInt) call IO_error(267_pInt)
if (subStepMinCryst <= 0.0_pReal) call IO_error(268_pInt)
if (subStepSizeCryst <= 0.0_pReal) call IO_error(268_pInt)
if (stepIncreaseCryst <= 0.0_pReal) call IO_error(268_pInt)
if (subStepMinHomog <= 0.0_pReal) call IO_error(268_pInt)
if (subStepSizeHomog <= 0.0_pReal) call IO_error(268_pInt)
if (stepIncreaseHomog <= 0.0_pReal) call IO_error(268_pInt)
if (rTol_crystalliteState <= 0.0_pReal) call IO_error(269_pInt)
if (rTol_crystalliteTemperature <= 0.0_pReal) call IO_error(276_pInt) !! oops !!
if (rTol_crystalliteStress <= 0.0_pReal) call IO_error(270_pInt)
if (aTol_crystalliteStress <= 0.0_pReal) call IO_error(271_pInt)
call IO_error(301_pInt,ext_msg='pert_method')
if (nHomog < 1_pInt) call IO_error(301_pInt,ext_msg='nHomog')
if (nMPstate < 1_pInt) call IO_error(301_pInt,ext_msg='nMPstate')
if (nCryst < 1_pInt) call IO_error(301_pInt,ext_msg='nCryst')
if (nState < 1_pInt) call IO_error(301_pInt,ext_msg='nState')
if (nStress < 1_pInt) call IO_error(301_pInt,ext_msg='nStress')
if (subStepMinCryst <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepMinCryst')
if (subStepSizeCryst <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepSizeCryst')
if (stepIncreaseCryst <= 0.0_pReal) call IO_error(301_pInt,ext_msg='stepIncreaseCryst')
if (subStepMinHomog <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepMinHomog')
if (subStepSizeHomog <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepSizeHomog')
if (stepIncreaseHomog <= 0.0_pReal) call IO_error(301_pInt,ext_msg='stepIncreaseHomog')
if (rTol_crystalliteState <= 0.0_pReal) call IO_error(301_pInt,ext_msg='rTol_crystalliteState')
if (rTol_crystalliteTemperature <= 0.0_pReal) call IO_error(301_pInt,ext_msg='rTol_crystalliteTemperature')
if (rTol_crystalliteStress <= 0.0_pReal) call IO_error(301_pInt,ext_msg='rTol_crystalliteStress')
if (aTol_crystalliteStress <= 0.0_pReal) call IO_error(301_pInt,ext_msg='aTol_crystalliteStress')
if (any(numerics_integrator <= 0_pInt) .or. any(numerics_integrator >= 6_pInt)) &
call IO_error(298_pInt)
call IO_error(301_pInt,ext_msg='integrator')
!* RGC parameters
if (absTol_RGC <= 0.0_pReal) call IO_error(272_pInt)
if (relTol_RGC <= 0.0_pReal) call IO_error(273_pInt)
if (absMax_RGC <= 0.0_pReal) call IO_error(274_pInt)
if (relMax_RGC <= 0.0_pReal) call IO_error(275_pInt)
if (pPert_RGC <= 0.0_pReal) call IO_error(276_pInt) !! oops !!
if (xSmoo_RGC <= 0.0_pReal) call IO_error(277_pInt)
if (viscPower_RGC < 0.0_pReal) call IO_error(278_pInt)
if (viscModus_RGC < 0.0_pReal) call IO_error(278_pInt)
if (refRelaxRate_RGC <= 0.0_pReal) call IO_error(278_pInt)
if (maxdRelax_RGC <= 0.0_pReal) call IO_error(288_pInt)
if (maxVolDiscr_RGC <= 0.0_pReal) call IO_error(289_pInt)
if (volDiscrMod_RGC < 0.0_pReal) call IO_error(289_pInt)
if (volDiscrPow_RGC <= 0.0_pReal) call IO_error(289_pInt)
if (absTol_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='absTol_RGC')
if (relTol_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='relTol_RGC')
if (absMax_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='absMax_RGC')
if (relMax_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='relMax_RGC')
if (pPert_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='pPert_RGC')
if (xSmoo_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='xSmoo_RGC')
if (viscPower_RGC < 0.0_pReal) call IO_error(301_pInt,ext_msg='viscPower_RGC')
if (viscModus_RGC < 0.0_pReal) call IO_error(301_pInt,ext_msg='viscModus_RGC')
if (refRelaxRate_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='refRelaxRate_RGC')
if (maxdRelax_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='maxdRelax_RGC')
if (maxVolDiscr_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='maxVolDiscr_RGC')
if (volDiscrMod_RGC < 0.0_pReal) call IO_error(301_pInt,ext_msg='volDiscrMod_RGC')
if (volDiscrPow_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='volDiscrPw_RGC')
!* spectral parameters
if (err_div_tol <= 0.0_pReal) call IO_error(49_pInt)
if (err_stress_tolrel <= 0.0_pReal) call IO_error(49_pInt)
if (itmax <= 1.0_pInt) call IO_error(49_pInt)
!* spectral parameters
if (err_div_tol <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_div_tol')
if (err_stress_tolrel <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_stress_tolrel')
if (itmax <= 1.0_pInt) call IO_error(301_pInt,ext_msg='itmax')
if (fixedSeed <= 0_pInt) then
!$OMP CRITICAL (write2out)