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

added and rearranged error messages, polished output and simplified calculation of f depending on fdot or L 

guessing along former trajectory is now on per default, 'guessreset' and 'dropguessing' switch it off.
This commit is contained in:
Martin Diehl 2011-07-14 09:37:31 +00:00
parent 56340fd487
commit 09ba92c26e
2 changed files with 85 additions and 82 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

129
code/DAMASK_spectral.f90
View File

@ -91,9 +91,9 @@ program DAMASK_spectral
real(pReal), dimension(3,3) :: ones, zeroes, temp33_Real, damper,& real(pReal), dimension(3,3) :: ones, zeroes, temp33_Real, damper,&
pstress, pstress_av, cstress_av, defgrad_av,& pstress, pstress_av, cstress_av, defgrad_av,&
defgradAim, defgradAimOld, defgradAimCorr, defgradAimCorrPrev,& defgradAim, defgradAimOld, defgradAimCorr, defgradAimCorrPrev,&
mask_stress, mask_defgrad mask_stress, mask_defgrad, deltaF
real(pReal), dimension(3,3,3,3) :: dPdF, c0, s0, c0_temp real(pReal), dimension(3,3,3,3) :: dPdF, c0, s0, c0_temp
real(pReal), dimension(9,9) :: s099 ! compliance in matrix notation ! real(pReal), dimension(9,9) :: s099 ! compliance in matrix notation
real(pReal), dimension(6) :: cstress ! cauchy stress in Mandel notation real(pReal), dimension(6) :: cstress ! cauchy stress in Mandel notation
real(pReal), dimension(6,6) :: dsde, c066, s066 ! Mandel notation of 4th order tensors real(pReal), dimension(6,6) :: dsde, c066, s066 ! Mandel notation of 4th order tensors
real(pReal), dimension(:,:,:,:,:), allocatable :: workfft, defgrad, defgradold real(pReal), dimension(:,:,:,:,:), allocatable :: workfft, defgrad, defgradold
@ -145,7 +145,7 @@ program DAMASK_spectral
print '(a,/,a)', 'Workingdir: ',trim(getSolverWorkingDirectoryName()) print '(a,/,a)', 'Workingdir: ',trim(getSolverWorkingDirectoryName())
print '(a,/,a)', 'SolverJobName: ',trim(getSolverJobName()) print '(a,/,a)', 'SolverJobName: ',trim(getSolverJobName())
if (.not. IO_open_file(unit,path)) call IO_error(45,ext_msg = path) if (.not. IO_open_file(unit,path)) call IO_error(30,ext_msg = path)
rewind(unit) rewind(unit)
do do
@ -183,89 +183,94 @@ if ((N_l + N_Fdot /= N_n).or.(N_n /= N_t)) & ! sanity check
allocate (bc_steps(N_Loadcases)); bc_steps = 0_pInt allocate (bc_steps(N_Loadcases)); bc_steps = 0_pInt
allocate (bc_logscale(N_Loadcases)); bc_logscale = .false. allocate (bc_logscale(N_Loadcases)); bc_logscale = .false.
allocate (bc_frequency(N_Loadcases)); bc_frequency = 1_pInt allocate (bc_frequency(N_Loadcases)); bc_frequency = 1_pInt
allocate (followFormeTrajectory(N_Loadcases)); followFormeTrajectory = .false. allocate (followFormeTrajectory(N_Loadcases)); followFormeTrajectory = .true.
rewind(unit) rewind(unit)
i = 0_pInt loadcase = 0_pInt
do do
read(unit,'(a1024)',END = 200) line read(unit,'(a1024)',END = 200) line
if (IO_isBlank(line)) cycle ! skip empty lines if (IO_isBlank(line)) cycle ! skip empty lines
i = i + 1 loadcase = loadcase + 1
posInput = IO_stringPos(line,maxNchunksInput) posInput = IO_stringPos(line,maxNchunksInput)
do j = 1,maxNchunksInput,2 do j = 1,maxNchunksInput,2
select case (IO_lc(IO_stringValue(line,posInput,j))) select case (IO_lc(IO_stringValue(line,posInput,j)))
case('l','velocitygrad') case('fdot') ! assign values for the deformation BC matrix (in case of given fdot)
velGradApplied(i) = .true.
valuevector = 0.0_pReal valuevector = 0.0_pReal
forall (k = 1:9) bc_maskvector(k) = IO_stringValue(line,posInput,j+k) /= '*' forall (k = 1:9) bc_maskvector(k) = IO_stringValue(line,posInput,j+k) /= '*'
do k = 1,9 do k = 1,9
if (bc_maskvector(k)) valuevector(k) = IO_floatValue(line,posInput,j+k) ! assign values for the velocity gradient matrix if (bc_maskvector(k)) valuevector(k) = IO_floatValue(line,posInput,j+k)
enddo enddo
bc_mask(:,:,1,i) = transpose(reshape(bc_maskvector,(/3,3/))) bc_mask(:,:,1,loadcase) = transpose(reshape(bc_maskvector,(/3,3/)))
bc_deformation(:,:,i) = math_transpose3x3(reshape(valuevector,(/3,3/))) bc_deformation(:,:,loadcase) = math_transpose3x3(reshape(valuevector,(/3,3/)))
case('fdot') case('l','velocitygrad') ! assign values for the deformation BC matrix (in case of given L)
velGradApplied(i) = .false. velGradApplied(loadcase) = .true. ! in case of given L, set flag to true
valuevector = 0.0_pReal valuevector = 0.0_pReal
forall (k = 1:9) bc_maskvector(k) = IO_stringValue(line,posInput,j+k) /= '*' forall (k = 1:9) bc_maskvector(k) = IO_stringValue(line,posInput,j+k) /= '*'
do k = 1,9 do k = 1,9
if (bc_maskvector(k)) valuevector(k) = IO_floatValue(line,posInput,j+k) ! assign values for the matrix of time derivative of defgrad if (bc_maskvector(k)) valuevector(k) = IO_floatValue(line,posInput,j+k)
enddo enddo
bc_mask(:,:,1,i) = transpose(reshape(bc_maskvector,(/3,3/))) bc_mask(:,:,1,loadcase) = transpose(reshape(bc_maskvector,(/3,3/)))
bc_deformation(:,:,i) = math_transpose3x3(reshape(valuevector,(/3,3/))) bc_deformation(:,:,loadcase) = math_transpose3x3(reshape(valuevector,(/3,3/)))
case('s', 'stress', 'pk1', 'piolakirchhoff') case('s', 'stress', 'pk1', 'piolakirchhoff')
valuevector = 0.0_pReal valuevector = 0.0_pReal
forall (k = 1:9) bc_maskvector(k) = IO_stringValue(line,posInput,j+k) /= '*' forall (k = 1:9) bc_maskvector(k) = IO_stringValue(line,posInput,j+k) /= '*'
do k = 1,9 do k = 1,9
if (bc_maskvector(k)) valuevector(k) = IO_floatValue(line,posInput,j+k) ! assign values for the bc_stress matrix if (bc_maskvector(k)) valuevector(k) = IO_floatValue(line,posInput,j+k) ! assign values for the bc_stress matrix
enddo enddo
bc_mask(:,:,2,i) = transpose(reshape(bc_maskvector,(/3,3/))) bc_mask(:,:,2,loadcase) = transpose(reshape(bc_maskvector,(/3,3/)))
bc_stress(:,:,i) = math_transpose3x3(reshape(valuevector,(/3,3/))) bc_stress(:,:,loadcase) = math_transpose3x3(reshape(valuevector,(/3,3/)))
case('t','time','delta') ! increment time case('t','time','delta') ! increment time
bc_timeIncrement(i) = IO_floatValue(line,posInput,j+1) bc_timeIncrement(loadcase) = IO_floatValue(line,posInput,j+1)
case('n','incs','increments','steps') ! bc_steps case('n','incs','increments','steps') ! bc_steps
bc_steps(i) = IO_intValue(line,posInput,j+1) bc_steps(loadcase) = IO_intValue(line,posInput,j+1)
case('logincs','logsteps') ! true, if log scale case('logincs','logsteps') ! true, if log scale
bc_steps(i) = IO_intValue(line,posInput,j+1) bc_steps(loadcase) = IO_intValue(line,posInput,j+1)
bc_logscale(i) = .true. bc_logscale(loadcase) = .true.
case('f','freq','frequency') ! frequency of result writings case('f','freq','frequency') ! frequency of result writings
bc_frequency(i) = IO_intValue(line,posInput,j+1) bc_frequency(loadcase) = IO_intValue(line,posInput,j+1)
case('continue','keepguessing','followtrajectory','followformertrajectory') case('guessreset','dropguessing')
followFormeTrajectory(i) = .true. ! continue to predict deformation along former trajectory followFormeTrajectory(loadcase) = .false. ! do not continue to predict deformation along former trajectory
end select end select
enddo; enddo enddo; enddo
200 close(unit) 200 close(unit)
do i = 1, N_Loadcases ! consistency checks if (followFormeTrajectory(1)) then
if (any(bc_mask(:,:,1,i) == bc_mask(:,:,2,i))) call IO_error(46,i) ! exclusive or masking only call IO_warning(33,loadcase) ! cannot guess along trajectory for first step of first loadcase
if (velGradApplied(i) == .true.) then followFormeTrajectory(1) = .false.
endif
do loadcase = 1, N_Loadcases ! consistency checks and output
print *, '------------------------------------------------------'
print '(a,i5)', 'Loadcase:', loadcase
if (.not. followFormeTrajectory(loadcase)) &
print '(a)', 'drop guessing along trajectory'
if (any(bc_mask(:,:,1,loadcase) == bc_mask(:,:,2,loadcase)))&
call IO_error(31,loadcase) ! exclusive or masking only
if (velGradApplied(loadcase)) then
do j = 1, 3 do j = 1, 3
if (any(bc_mask(j,:,1,i) == .true.) .and. any(bc_mask(j,:,1,i) == .false.)) & if (any(bc_mask(j,:,1,loadcase) == .true.) .and.&
call IO_error(46,i) ! ToDo: change message any(bc_mask(j,:,1,loadcase) == .false.)) call IO_error(32,loadcase) ! each line should be either fully or not at all defined
enddo enddo
endif print '(a,/,3(3(f12.6,x)/))','L:' ,math_transpose3x3(bc_deformation(:,:,loadcase))
if (followFormeTrajectory(1) == .true.) call IO_error(47,i) ! ToDo: Change message print '(a,/,3(3(l,x)/))', 'bc_mask for L:',transpose(bc_mask(:,:,1,loadcase))
if (bc_timeIncrement(i) < 0.0_pReal) call IO_error(47,i) ! negative time increment
if (bc_steps(i) < 1_pInt) call IO_error(48,i) ! non-positive increment count
if (bc_frequency(i) < 1_pInt) call IO_error(49,i) ! non-positive result frequency
print '(a,i5)', 'Loadcase:', i
if (velGradApplied(i) == .true.) then
print '(a,/,3(3(f12.6,x)/))','L:' ,math_transpose3x3(bc_deformation(:,:,i))
print '(a,/,3(3(l,x)/))', 'bc_mask for L:',transpose(bc_mask(:,:,1,i))
else else
print '(a,/,3(3(f12.6,x)/))','Fdot:' ,math_transpose3x3(bc_deformation(:,:,i)) print '(a,/,3(3(f12.6,x)/))','Fdot:' ,math_transpose3x3(bc_deformation(:,:,loadcase))
print '(a,/,3(3(l,x)/))', 'bc_mask for Fdot:',transpose(bc_mask(:,:,1,i)) print '(a,/,3(3(l,x)/))', 'bc_mask for Fdot:',transpose(bc_mask(:,:,1,loadcase))
endif endif
print '(a,/,3(3(f12.6,x)/))','bc_stress:',math_transpose3x3(bc_stress(:,:,i)) print '(a,/,3(3(f12.6,x)/))','bc_stress:',math_transpose3x3(bc_stress(:,:,loadcase))
print '(a,/,3(3(l,x)/))', 'bc_mask for stress:' ,transpose(bc_mask(:,:,2,i)) print '(a,/,3(3(l,x)/))', 'bc_mask for stress:' ,transpose(bc_mask(:,:,2,loadcase))
print '(a,f12.6)','time: ',bc_timeIncrement(i) if (bc_timeIncrement(loadcase) < 0.0_pReal) call IO_error(34,loadcase) ! negative time increment
print '(a,i6)','incs: ',bc_steps(i) print '(a,f12.6)','time: ',bc_timeIncrement(loadcase)
print '(a,i6)','freq: ',bc_frequency(i) if (bc_steps(loadcase) < 1_pInt) call IO_error(35,loadcase) ! non-positive increment count
print *, '' print '(a,i6)','incs: ',bc_steps(loadcase)
if (bc_frequency(loadcase) < 1_pInt) call IO_error(36,loadcase) ! non-positive result frequency
print '(a,i6)','freq: ',bc_frequency(loadcase)
enddo enddo
!read header of geom file to get the information needed before the complete geom file is intepretated by mesh.f90 !read header of geom file to get the information needed before the complete geom file is intepretated by mesh.f90
path = getModelName() path = getModelName()
print *, '------------------------------------------------------'
print '(a,a)', 'GeomName: ',trim(path) print '(a,a)', 'GeomName: ',trim(path)
if (.not. IO_open_file(unit,trim(path)//InputFileExtension)) call IO_error(101,ext_msg = trim(path)//InputFileExtension) if (.not. IO_open_file(unit,trim(path)//InputFileExtension)) call IO_error(101,ext_msg = trim(path)//InputFileExtension)
@ -417,7 +422,7 @@ if ((N_l + N_Fdot /= N_n).or.(N_n /= N_t)) & ! sanity check
!************************************************************* !*************************************************************
time0 = time ! loadcase start time time0 = time ! loadcase start time
if (followFormeTrajectory(loadcase) == .true.) then if (followFormeTrajectory(loadcase)) then
guessmode = 1.0_pReal guessmode = 1.0_pReal
else else
guessmode = 0.0_pReal ! change of load case, homogeneous guess for the first step guessmode = 0.0_pReal ! change of load case, homogeneous guess for the first step
@ -426,6 +431,7 @@ if ((N_l + N_Fdot /= N_n).or.(N_n /= N_t)) & ! sanity check
mask_defgrad = merge(ones,zeroes,bc_mask(:,:,1,loadcase)) mask_defgrad = merge(ones,zeroes,bc_mask(:,:,1,loadcase))
mask_stress = merge(ones,zeroes,bc_mask(:,:,2,loadcase)) mask_stress = merge(ones,zeroes,bc_mask(:,:,2,loadcase))
damper = ones/10 damper = ones/10
deltaF = bc_deformation(:,:,loadcase) ! only valid for given fDot. will be overwritten later in case L is given
!************************************************************* !*************************************************************
! loop oper steps defined in input file for current loadcase ! loop oper steps defined in input file for current loadcase
do step = 1, bc_steps(loadcase) do step = 1, bc_steps(loadcase)
@ -448,32 +454,25 @@ if ((N_l + N_Fdot /= N_n).or.(N_n /= N_t)) & ! sanity check
time = time + timeinc time = time + timeinc
! update macroscopic deformation gradient (defgrad BC) ! update macroscopic deformation gradient (defgrad BC)
if (velGradApplied(loadcase)) & ! calculate deltaF from given L and current F
deltaF = math_mul33x33(bc_deformation(:,:,loadcase), defgradAim)
temp33_Real = defgradAim temp33_Real = defgradAim
if (velGradApplied(loadcase) == .true.) then ! calculate from given L and current F
defgradAim = defgradAim & defgradAim = defgradAim &
+ guessmode * mask_stress * (defgradAim - defgradAimOld) & + guessmode * mask_stress * (defgradAim - defgradAimOld) &
+ math_mul33x33(bc_deformation(:,:,loadcase), defgradAim)*timeinc + mask_defgrad * deltaF * timeinc
else ! calculate from given Fdot and current F
defgradAim = defgradAim &
+ guessmode * mask_stress * (defgradAim - defgradAimOld) &
+ mask_defgrad * bc_deformation(:,:,loadcase) * timeinc
endif
defgradAimOld = temp33_Real defgradAimOld = temp33_Real
! update local deformation gradient ! update local deformation gradient
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)
if (velGradApplied(loadcase) == .true.) then ! using velocity gradient to calculate new deformation gradient (if not guessing)
temp33_Real = defgrad(i,j,k,:,:) temp33_Real = defgrad(i,j,k,:,:)
defgrad(i,j,k,:,:) = defgrad(i,j,k,:,:) & ! decide if guessing along former trajectory or apply homogeneous addon if (velGradApplied(loadcase)) & ! using velocity gradient to calculate new deformation gradient (if not guessing)
deltaF = math_mul33x33(bc_deformation(:,:,loadcase),defgradold(i,j,k,:,:))
defgrad(i,j,k,:,:) = defgrad(i,j,k,:,:) & ! decide if guessing along former trajectory or apply homogeneous addon (addon only for applied deformation)
+ guessmode * (defgrad(i,j,k,:,:) - defgradold(i,j,k,:,:))& + guessmode * (defgrad(i,j,k,:,:) - defgradold(i,j,k,:,:))&
+ (1.0_pReal-guessmode) * math_mul33x33(bc_deformation(:,:,loadcase),defgradold(i,j,k,:,:))*timeinc + (1.0_pReal-guessmode) * mask_defgrad * deltaF *timeinc
defgradold(i,j,k,:,:) = temp33_Real defgradold(i,j,k,:,:) = temp33_Real
else ! using time derivative of deformation gradient to calculate new deformation gradient (if not guessing)
temp33_Real = defgrad(i,j,k,:,:)
defgrad(i,j,k,:,:) = defgrad(i,j,k,:,:) & ! decide if guessing along former trajectory or apply homogeneous addon
+ guessmode * (defgrad(i,j,k,:,:) - defgradold(i,j,k,:,:))&
+ (1.0_pReal-guessmode) * bc_deformation(:,:,loadcase) * timeinc
defgradold(i,j,k,:,:) = temp33_Real
endif
enddo; enddo; enddo enddo; enddo; enddo
guessmode = 1.0_pReal ! keep guessing along former trajectory during same loadcase guessmode = 1.0_pReal ! keep guessing along former trajectory during same loadcase

24
code/IO.f90
View File

@ -1154,6 +1154,18 @@ endfunction
character(len=120) msg character(len=120) msg
select case (ID) select case (ID)
case (30)
msg = 'error opening spectral loadcase'
case (31)
msg = 'mask consistency violated in spectral loadcase'
case (32)
msg = 'ill defined L (each line should be either fully or not at all defined) in spectral loadcase'
case (34)
msg = 'negative time increment in spectral loadcase'
case (35)
msg = 'Non-positive increments in spectral loadcase'
case (36)
msg = 'Non-positive result frequency in spectral loadcase'
case (40) case (40)
msg = 'path rectification error' msg = 'path rectification error'
case (41) case (41)
@ -1164,16 +1176,6 @@ endfunction
msg = 'resolution error in spectral mesh' msg = 'resolution error in spectral mesh'
case (44) case (44)
msg = 'dimension error in spectral mesh' msg = 'dimension error in spectral mesh'
case (45)
msg = 'error opening spectral loadcase'
case (46)
msg = 'mask consistency violated in spectral loadcase'
case (47)
msg = 'negative time increment in spectral loadcase'
case (48)
msg = 'Non-positive increments in spectral loadcase'
case (49)
msg = 'Non-positive result frequency in spectral loadcase'
case (50) case (50)
msg = 'Error writing constitutive output description' msg = 'Error writing constitutive output description'
case (100) case (100)
@ -1417,6 +1419,8 @@ endfunction
character(len=80) msg character(len=80) msg
select case (ID) select case (ID)
case (33)
msg = 'cannot guess along trajectory for first step of first loadcase'
case (101) case (101)
msg = '+ crystallite debugging off... +' msg = '+ crystallite debugging off... +'
case (600) case (600)