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

mainly fixed error in output of spectral results (1:N,…) instead of (N,…) 

rearranged some logic here and there.

(hopefully) improved readability of debug/standard output.

restarting logic would need some discussion with Martin/Krishna still…
This commit is contained in:
Philip Eisenlohr 2011-12-04 10:01:32 +00:00
parent 9d3f7b8d3d
commit efadf9f728
2 changed files with 165 additions and 155 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

203
code/DAMASK_spectral.f90
View File

@ -96,8 +96,9 @@ program DAMASK_spectral
! variables storing information from geom file ! variables storing information from geom file
real(pReal) :: wgt real(pReal) :: wgt
real(pReal), dimension(3) :: geomdimension = 0.0_pReal ! physical dimension of volume element in each direction real(pReal), dimension(3) :: geomdimension = 0.0_pReal ! physical dimension of volume element per direction
integer(pInt) :: homog ! homogenization scheme used integer(pInt) :: Npoints,& ! number of Fourier points
homog ! homogenization scheme used
integer(pInt), dimension(3) :: res = 1_pInt ! resolution (number of Fourier points) in each direction integer(pInt), dimension(3) :: res = 1_pInt ! resolution (number of Fourier points) in each direction
logical :: spectralPictureMode = .false. ! indicating 1 to 1 mapping of FP to microstructure logical :: spectralPictureMode = .false. ! indicating 1 to 1 mapping of FP to microstructure
@ -132,19 +133,19 @@ program DAMASK_spectral
real(pReal), dimension(3,3), parameter :: ones = 1.0_pReal, zeroes = 0.0_pReal real(pReal), dimension(3,3), parameter :: ones = 1.0_pReal, zeroes = 0.0_pReal
complex(pReal), dimension(3,3) :: temp33_Complex complex(pReal), dimension(3,3) :: temp33_Complex
real(pReal), dimension(3,3) :: temp33_Real real(pReal), dimension(3,3) :: temp33_Real
integer(pInt) :: i, j, k, l, m, n, p integer(pInt) :: i, j, k, l, m, n, p, errorID
integer(pInt) :: N_Loadcases, loadcase, step, iter, ielem, CPFEM_mode, & integer(pInt) :: N_Loadcases, loadcase, step, iter, ielem, CPFEM_mode, &
ierr, notConvergedCounter = 0_pInt, totalStepsCounter = 0_pInt ierr, notConvergedCounter = 0_pInt, totalStepsCounter = 0_pInt
logical :: errmatinv, regrid = .false. logical :: errmatinv, regrid = .false.
real(pReal) :: defgradDet, defgradDetMax, defgradDetMin real(pReal) :: defgradDet, defgradDetMax, defgradDetMin
real(pReal) :: correctionFactor real(pReal) :: correctionFactor
! debuging variables ! --- debugging variables
real(pReal), dimension(:,:,:,:), allocatable :: divergence real(pReal), dimension(:,:,:,:), allocatable :: divergence
real(pReal) :: p_real_avg, err_div_max, err_real_div_avg, err_real_div_max real(pReal) :: p_real_avg, err_div_max, err_real_div_avg, err_real_div_max
logical :: debugGeneral = .false., debugDivergence = .false., debugRestart = .false. logical :: debugGeneral = .false., debugDivergence = .false., debugRestart = .false.
! initialize default value for loadcase ! --- initialize default value for loadcase
bc_init%deformation = zeroes; bc_init%stress = zeroes; bc_init%rotation = zeroes bc_init%deformation = zeroes; bc_init%stress = zeroes; bc_init%rotation = zeroes
bc_init%timeIncrement = 0.0_pReal; bc_init%temperature = 300.0_pReal bc_init%timeIncrement = 0.0_pReal; bc_init%temperature = 300.0_pReal
bc_init%steps = 0_pInt; bc_init%logscale = 0_pInt bc_init%steps = 0_pInt; bc_init%logscale = 0_pInt
@ -154,14 +155,13 @@ program DAMASK_spectral
bc_init%followFormerTrajectory = .true. bc_init%followFormerTrajectory = .true.
bc_init%rotation = math_I3 ! assume no rotation bc_init%rotation = math_I3 ! assume no rotation
! Initializing model size independed parameters ! --- initializing model size independed parameters
!$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS !$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS
if (.not.(command_argument_count()==4 .or. command_argument_count()==6)) &! check for correct number of given arguments if (.not.(command_argument_count()==4 .or. command_argument_count()==6)) &! check for correct number of given arguments
call IO_error(error_ID=102_pInt) call IO_error(error_ID=102_pInt)
call DAMASK_interface_init() call DAMASK_interface_init()
!$OMP CRITICAL (write2out)
print '(a)', '' print '(a)', ''
print '(a,a)', ' <<<+- DAMASK_spectral init -+>>>' print '(a,a)', ' <<<+- DAMASK_spectral init -+>>>'
print '(a,a)', ' $Id$' print '(a,a)', ' $Id$'
@ -169,7 +169,6 @@ program DAMASK_spectral
print '(a,a)', ' Working Directory: ',trim(getSolverWorkingDirectoryName()) print '(a,a)', ' Working Directory: ',trim(getSolverWorkingDirectoryName())
print '(a,a)', ' Solver Job Name: ',trim(getSolverJobName()) print '(a,a)', ' Solver Job Name: ',trim(getSolverJobName())
print '(a)', '' print '(a)', ''
!$OMP END CRITICAL (write2out)
! Reading the loadcase file and allocate variables for loadcases ! Reading the loadcase file and allocate variables for loadcases
path = getLoadcaseName() path = getLoadcaseName()
@ -199,7 +198,7 @@ program DAMASK_spectral
allocate (bc(N_Loadcases)) allocate (bc(N_Loadcases))
! Reading the loadcase and assign values to the allocated data structure ! --- reading the loadcase and assign values to the allocated data structure
rewind(myUnit) rewind(myUnit)
loadcase = 0_pInt loadcase = 0_pInt
do do
@ -267,7 +266,7 @@ program DAMASK_spectral
101 close(myUnit) 101 close(myUnit)
!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()
if (.not. IO_open_file(myUnit,trim(path)//InputFileExtension))& if (.not. IO_open_file(myUnit,trim(path)//InputFileExtension))&
@ -325,94 +324,83 @@ program DAMASK_spectral
mod(res(2),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) (mod(res(3),2_pInt)/=0_pInt .and. res(3)/= 1_pInt)) call IO_error(error_ID = 103_pInt)
! Initialization of CPFEM_general (= constitutive law) Npoints = res(1)*res(2)*res(3)
! --- initialization of CPFEM_general (= constitutive law)
call CPFEM_initAll(bc(1)%temperature,1_pInt,1_pInt) call CPFEM_initAll(bc(1)%temperature,1_pInt,1_pInt)
! Get debugging parameters ! --- debugging parameters
if (iand(spectral_debug_verbosity,1_pInt)==1_pInt) debugGeneral = .true. if (iand(spectral_debug_verbosity,1_pInt)==1_pInt) debugGeneral = .true.
if (iand(spectral_debug_verbosity,2_pInt)==2_pInt) debugDivergence = .true. if (iand(spectral_debug_verbosity,2_pInt)==2_pInt) debugDivergence = .true.
if (iand(spectral_debug_verbosity,4_pInt)==4_pInt) debugRestart = .true. if (iand(spectral_debug_verbosity,4_pInt)==4_pInt) debugRestart = .true.
!Output of geometry ! --- output of geometry
!$OMP CRITICAL (write2out)
print '(a)', '' print '(a)', ''
print '(a)', '#############################################################' print '(a)', '#############################################################'
print '(a)', 'DAMASK spectral:' print '(a)', 'DAMASK spectral:'
print '(a)', 'The spectral method boundary value problem solver for' print '(a)', 'The spectral method boundary value problem solver for'
print '(a)', 'the Duesseldorf Advanced Material Simulation Kit' print '(a)', 'the Duesseldorf Advanced Material Simulation Kit'
print '(a)', '#############################################################' print '(a)', '#############################################################'
print '(a,a)', 'Geom File Name: ',trim(path)//'.geom' print '(a,a)', 'geometry file: ',trim(path)//'.geom'
print '(a)', '=============================================================' print '(a)', '============================================================='
print '(a,i12,i12,i12)','resolution a b c:', res print '(a,i12,i12,i12)','resolution a b c:', res
print '(a,f12.5,f12.5,f12.5)','dimension x y z:', geomdimension print '(a,f12.5,f12.5,f12.5)','dimension x y z:', geomdimension
print '(a,i5)','homogenization: ',homog print '(a,i5)','homogenization: ',homog
print '(a,L)', 'spectralPictureMode: ',spectralPictureMode print '(a,L)', 'spectralPictureMode: ',spectralPictureMode
print '(a)', '#############################################################' print '(a)', '#############################################################'
print '(a,a)','Loadcase File Name: ',trim(getLoadcaseName()) print '(a,a)', 'loadcase file: ',trim(getLoadcaseName())
!$OMP END CRITICAL (write2out)
if (bc(1)%followFormerTrajectory) then if (bc(1)%followFormerTrajectory) then
call IO_warning(warning_ID = 33_pInt) ! cannot guess along trajectory for first step of first loadcase call IO_warning(warning_ID = 33_pInt) ! cannot guess along trajectory for first step of first loadcase
bc(1)%followFormerTrajectory = .false. bc(1)%followFormerTrajectory = .false.
endif endif
! consistency checks and output of loadcase ! --- consistency checks and output of loadcase
errorID = 0_pInt
do loadcase = 1_pInt, N_Loadcases do loadcase = 1_pInt, N_Loadcases
!$OMP CRITICAL (write2out)
print '(a)', '============================================================='
print '(a,i5)', 'Loadcase: ', loadcase
write (loadcase_string, '(i3)' ) loadcase write (loadcase_string, '(i3)' ) loadcase
if (.not. bc(loadcase)%followFormerTrajectory) &
print '(a)', 'Drop Guessing Along Trajectory' print '(a)', '============================================================='
!$OMP END CRITICAL (write2out) print '(a,i5)', 'loadcase: ', loadcase
if (any(bc(loadcase)%maskStress .eqv. bc(loadcase)%maskDeformation))& ! exclusive or masking only
call IO_error(error_ID=31_pInt,ext_msg=loadcase_string) if (.not. bc(loadcase)%followFormerTrajectory) print '(a)', 'drop guessing along trajectory'
if (any(bc(loadcase)%maskStress.and.transpose(bc(loadcase)%maskStress).and.& !checking if no rotation is allowed by stress BC
reshape((/.false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false./),(/3,3/))))&
call IO_error(error_ID=38_pInt,ext_msg=loadcase_string)
if (bc(loadcase)%velGradApplied) then if (bc(loadcase)%velGradApplied) then
do j = 1_pInt, 3_pInt do j = 1_pInt, 3_pInt
if (any(bc(loadcase)%maskDeformation(j,1:3) .eqv. .true.) .and.& if (any(bc(loadcase)%maskDeformation(j,1:3) == .true.) .and. &
any(bc(loadcase)%maskDeformation(j,1:3) .eqv. .false.)) call IO_error(error_ID=32_pInt,ext_msg=loadcase_string) ! each line should be either fully or not at all defined any(bc(loadcase)%maskDeformation(j,1:3) == .false.)) errorID = 32_pInt ! each row should be either fully or not at all defined
enddo enddo
!$OMP CRITICAL (write2out) print '(a)','velocity gradient:'
print '(a)','Velocity Gradient:'
!$OMP END CRITICAL (write2out)
else else
!$OMP CRITICAL (write2out) print '(a)','deformation gradient rate:'
print '(a)','Change of Deformation Gradient:'
!$OMP END CRITICAL (write2out)
endif endif
!$OMP CRITICAL (write2out)
print '(3(3(f12.6,x)/)\)', merge(math_transpose3x3(bc(loadcase)%deformation),& print '(3(3(f12.6,x)/)\)', merge(math_transpose3x3(bc(loadcase)%deformation),&
reshape(spread(DAMASK_NaN,1,9),(/3,3/)),transpose(bc(loadcase)%maskDeformation)) reshape(spread(DAMASK_NaN,1,9),(/3,3/)),transpose(bc(loadcase)%maskDeformation))
print '(a,/,3(3(f12.6,x)/)\)','Stress Boundary Condition/MPa:',merge(math_transpose3x3(bc(loadcase)%stress),& print '(a,/,3(3(f12.6,x)/)\)','stress / GPa:',1e-9*merge(math_transpose3x3(bc(loadcase)%stress),&
reshape(spread(DAMASK_NaN,1,9),(/3,3/)),& reshape(spread(DAMASK_NaN,1,9),(/3,3/)),&
transpose(bc(loadcase)%maskStress))*1e-6 transpose(bc(loadcase)%maskStress))
!$OMP END CRITICAL (write2out)
if (any(abs(math_mul33x33(bc(loadcase)%rotation,math_transpose3x3(bc(loadcase)%rotation))-math_I3)& ! given rotation matrix contains strain
>reshape(spread(rotation_tol,1,9),(/3,3/)))&
.or. abs(math_det3x3(bc(loadcase)%rotation))>1.0_pReal + rotation_tol) call IO_error(error_ID=46_pInt,ext_msg=loadcase_string)
!$OMP CRITICAL (write2out)
if (any(bc(loadcase)%rotation /= math_I3)) & if (any(bc(loadcase)%rotation /= math_I3)) &
print '(a,3(3(f12.6,x)/)\)','Rotation of BCs:',math_transpose3x3(bc(loadcase)%rotation) print '(a,3(3(f12.6,x)/)\)','rotation of loadframe:',math_transpose3x3(bc(loadcase)%rotation)
!$OMP END CRITICAL (write2out) print '(a,f12.6)','temperature:',bc(loadcase)%temperature
if (bc(loadcase)%timeIncrement < 0.0_pReal) call IO_error(error_ID=34_pInt,ext_msg=loadcase_string) ! negative time increment print '(a,f12.6)','time: ',bc(loadcase)%timeIncrement
!$OMP CRITICAL (write2out) print '(a,i5)','steps: ',bc(loadcase)%steps
print '(a,f12.6)','Temperature:',bc(loadcase)%temperature print '(a,i5)','output frequency: ',bc(loadcase)%outputfrequency
print '(a,f12.6)','Time: ',bc(loadcase)%timeIncrement print '(a,i5)','restart frequency: ',bc(loadcase)%restartfrequency
!$OMP END CRITICAL (write2out)
if (bc(loadcase)%steps < 1_pInt) call IO_error(error_ID=35_pInt,ext_msg=loadcase_string) ! non-positive increment count if (any(bc(loadcase)%maskStress == bc(loadcase)%maskDeformation)) errorID = 31 ! exclusive or masking only
!$OMP CRITICAL (write2out) if (any(bc(loadcase)%maskStress .and. transpose(bc(loadcase)%maskStress) .and. &
print '(a,i5)','Steps: ',bc(loadcase)%steps reshape((/.false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false./),(/3,3/)))) &
!$OMP END CRITICAL (write2out) errorID = 38_pInt ! no rotation is allowed by stress BC
if (bc(loadcase)%outputfrequency < 1_pInt) call IO_error(error_ID=36_pInt,ext_msg=loadcase_string) ! non-positive result frequency if (any(abs(math_mul33x33(bc(loadcase)%rotation,math_transpose3x3(bc(loadcase)%rotation))-math_I3)&
!$OMP CRITICAL (write2out) > reshape(spread(rotation_tol,1,9),(/3,3/)))&
print '(a,i5)','Freq. of Results Output: ',bc(loadcase)%outputfrequency .or. abs(math_det3x3(bc(loadcase)%rotation)) > 1.0_pReal + rotation_tol) &
!$OMP END CRITICAL (write2out) errorID = 46_pInt ! given rotation matrix contains strain
if (bc(loadcase)%restartfrequency < 1_pInt) call IO_error(error_ID=39_pInt,ext_msg=loadcase_string) ! non-positive restart frequency if (bc(loadcase)%timeIncrement < 0.0_pReal) errorID = 34_pInt ! negative time increment
!$OMP CRITICAL (write2out) if (bc(loadcase)%steps < 1_pInt) errorID = 35_pInt ! non-positive increment count
print '(a,i5)','Freq. of Restart Information Output: ',bc(loadcase)%restartfrequency if (bc(loadcase)%outputfrequency < 1_pInt) errorID = 36_pInt ! non-positive result frequency
!$OMP END CRITICAL (write2out) if (bc(loadcase)%restartfrequency < 1_pInt) errorID = 39_pInt ! non-positive restart frequency
if (errorID > 0_pInt) call IO_error(error_ID = errorID, ext_msg = loadcase_string)
enddo enddo
! Initialization of fftw (see manual on fftw.org for more details) ! Initialization of fftw (see manual on fftw.org for more details)
@ -492,7 +480,7 @@ program DAMASK_spectral
allocate (workfft(res(1)+2,res(2),res(3),3,3)); workfft = 0.0_pReal allocate (workfft(res(1)+2,res(2),res(3),3,3)); workfft = 0.0_pReal
if (debugDivergence) allocate (divergence(res(1)+2,res(2),res(3),3)); divergence = 0.0_pReal if (debugDivergence) allocate (divergence(res(1)+2,res(2),res(3),3)); divergence = 0.0_pReal
wgt = 1.0_pReal/real(res(1)*res(2)*res(3), pReal) wgt = 1.0_pReal/real(Npoints, pReal)
call dfftw_plan_many_dft_r2c(fftw_plan(1),3,(/res(1),res(2),res(3)/),9,& call dfftw_plan_many_dft_r2c(fftw_plan(1),3,(/res(1),res(2),res(3)/),9,&
workfft,(/res(1) +2_pInt,res(2),res(3)/),1,(res(1) +2_pInt)*res(2)*res(3),& workfft,(/res(1) +2_pInt,res(2),res(3)/),1,(res(1) +2_pInt)*res(2)*res(3),&
workfft,(/res(1)/2_pInt+1_pInt,res(2),res(3)/),1,(res(1)/2_pInt+1_pInt)*res(2)*res(3),fftw_planner_flag) workfft,(/res(1)/2_pInt+1_pInt,res(2),res(3)/),1,(res(1)/2_pInt+1_pInt)*res(2)*res(3),fftw_planner_flag)
@ -541,7 +529,7 @@ program DAMASK_spectral
if (debugGeneral) then if (debugGeneral) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write (6,*) 'First Call to CPFEM_general finished' write (6,*) 'first call to CPFEM_general finished'
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
@ -558,7 +546,7 @@ program DAMASK_spectral
xi(2,i,j,k) = real(k_s(2), pReal)/geomdimension(2) xi(2,i,j,k) = real(k_s(2), pReal)/geomdimension(2)
xi(1,i,j,k) = real(k_s(1), pReal)/geomdimension(1) xi(1,i,j,k) = real(k_s(1), pReal)/geomdimension(1)
enddo; enddo; enddo enddo; enddo; enddo
! remove highest frequencies for calculation of divergence (CAREFULL, they will be used for pre calculatet gamma operator!) ! remove highest frequencies for calculation of divergence (CAREFUL, they will be used for pre calculatet gamma operator!)
do k = 1_pInt ,res(3); do j = 1_pInt ,res(2); do i = 1_pInt,res(1)/2_pInt + 1_pInt do k = 1_pInt ,res(3); do j = 1_pInt ,res(2); do i = 1_pInt,res(1)/2_pInt + 1_pInt
if(k==res(3)/2_pInt+1_pInt) xi(3,i,j,k)= 0.0_pReal if(k==res(3)/2_pInt+1_pInt) xi(3,i,j,k)= 0.0_pReal
if(j==res(2)/2_pInt+1_pInt) xi(2,i,j,k)= 0.0_pReal if(j==res(2)/2_pInt+1_pInt) xi(2,i,j,k)= 0.0_pReal
@ -605,7 +593,7 @@ program DAMASK_spectral
bc(1)%steps= bc(1)%steps - 1_pInt bc(1)%steps= bc(1)%steps - 1_pInt
write(538), 'startingIncrement', restartReadStep -1_pInt ! start with writing out the previous step write(538), 'startingIncrement', restartReadStep -1_pInt ! start with writing out the previous step
write(538), 'eoh' ! end of header write(538), 'eoh' ! end of header
write(538), materialpoint_results(materialpoint_sizeResults,1,res(1)*res(2)*res(3)) ! initial (non-deformed) results write(538), materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! initial (non-deformed) results
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
!************************************************************* !*************************************************************
@ -690,7 +678,7 @@ program DAMASK_spectral
print '(a)', '#############################################################' print '(a)', '#############################################################'
print '(A,I5.5,A,es12.6)', 'Increment ', totalStepsCounter, ' Time ',time print '(A,I5.5,A,es12.6)', 'Increment ', totalStepsCounter, ' Time ',time
if (restartWrite ) then if (restartWrite ) then
print '(A)', 'Writing converged Results of previous Step for Restart' print '(A)', 'writing converged results of previous step for restart'
if(IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(defgrad))) then ! and writing deformation gradient field to file if(IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(defgrad))) then ! and writing deformation gradient field to file
write (777,rec=1) defgrad write (777,rec=1) defgrad
close (777) close (777)
@ -704,14 +692,16 @@ program DAMASK_spectral
err_stress > err_stress_tol)) err_stress > err_stress_tol))
iter = iter + 1_pInt iter = iter + 1_pInt
!************************************************************* !*************************************************************
print '(a)', ''
print '(a)', '=============================================================' print '(a)', '============================================================='
print '(5(A,I5.5))', 'Loadcase ',loadcase,' Step ',step,'/',bc(loadcase)%steps,'@Iteration ',iter,'/',itmax print '(5(a,i5.5))', 'Loadcase ',loadcase,' Step ',step,'/',bc(loadcase)%steps,'@Iteration ',iter,'/',itmax
do n = 1_pInt,3_pInt; do m = 1_pInt,3_pInt do n = 1_pInt,3_pInt; do m = 1_pInt,3_pInt
defgrad_av(m,n) = sum(defgrad(1:res(1),1:res(2),1:res(3),m,n)) * wgt defgrad_av(m,n) = sum(defgrad(1:res(1),1:res(2),1:res(3),m,n)) * wgt
enddo; enddo enddo; enddo
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
print '(a,/,3(3(f12.7,x)/)\)', 'Deformation Gradient:',math_transpose3x3(defgrad_av) print '(a,/,3(3(f12.7,x)/)\)', 'deformation gradient:',math_transpose3x3(defgrad_av)
print '(A)', '... Update Stress Field (Constitutive Evaluation P(F)) ......' print '(l)', restartWrite
print '(a)', '... update stress field P(F) ................................'
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
defgradDetMax = -huge(1.0_pReal) defgradDetMax = -huge(1.0_pReal)
defgradDetMin = +huge(1.0_pReal) defgradDetMin = +huge(1.0_pReal)
@ -727,8 +717,8 @@ program DAMASK_spectral
cstress,dsde, pstress, dPdF) cstress,dsde, pstress, dPdF)
enddo; enddo; enddo enddo; enddo; enddo
print '(a,x,es10.4)' , 'Maximum Determinant of Deformation:', defgradDetMax print '(a,x,es10.4)' , 'max determinant of deformation:', defgradDetMax
print '(a,x,es10.4)' , 'Minimum Determinant of Deformation:', defgradDetMin print '(a,x,es10.4)' , 'min determinant of deformation:', defgradDetMin
workfft = 0.0_pReal ! needed because of the padding for FFTW workfft = 0.0_pReal ! needed because of the padding for FFTW
c_current = 0.0_pReal c_current = 0.0_pReal
@ -750,22 +740,27 @@ program DAMASK_spectral
enddo; enddo enddo; enddo
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
print '(a,/,3(3(f12.7,x)/)\)', 'Piola-Kirchhoff Stress / MPa: ',math_transpose3x3(pstress_av)/1.e6 print '(a,/,3(3(f12.7,x)/)\)', 'Piola-Kirchhoff stress / MPa: ',math_transpose3x3(pstress_av)/1.e6
err_stress_tol = 0.0_pReal err_stress_tol = 0.0_pReal
pstress_av_load = math_rotate_forward3x3(pstress_av,bc(loadcase)%rotation) pstress_av_load = math_rotate_forward3x3(pstress_av,bc(loadcase)%rotation)
if(size_reduced > 0_pInt) then ! calculate stress BC if applied if(size_reduced > 0_pInt) then ! calculate stress BC if applied
err_stress = maxval(abs(mask_stress * (pstress_av_load - bc(loadcase)%stress))) ! maximum deviaton (tensor norm not applicable) err_stress = maxval(abs(mask_stress * (pstress_av_load - bc(loadcase)%stress))) ! maximum deviaton (tensor norm not applicable)
err_stress_tol = maxval(abs(mask_defgrad * pstress_av_load)) * err_stress_tolrel ! don't use any tensor norm because the comparison should be coherent err_stress_tol = maxval(abs(mask_defgrad * pstress_av_load)) * err_stress_tolrel ! don't use any tensor norm because the comparison should be coherent
print '(A)', '... Correcting Deformation Gradient to Fullfill BCs .........' print '(a)', ''
print '(2(a,es10.4))', 'Error Stress = ',err_stress, ', Tol. = ', err_stress_tol print '(a,es10.4,a,f6.2)', 'error stress = ',err_stress, ', ', err_stress/err_stress_tol
print '(a)', '... correcting deformation gradient to fulfill BCs ..........'
defgradAimCorr = - math_mul3333xx33(s_prev, ((pstress_av_load - bc(loadcase)%stress))) ! residual on given stress components defgradAimCorr = - math_mul3333xx33(s_prev, ((pstress_av_load - bc(loadcase)%stress))) ! residual on given stress components
defgradAim = defgradAim + defgradAimCorr defgradAim = defgradAim + defgradAimCorr
print '(a,/,3(3(f12.7,x)/)\)', 'New Deformation Aim: ',math_transpose3x3(math_rotate_backward3x3(& print '(a,/,3(3(f12.7,x)/)\)', 'new deformation aim: ',&
defgradAim,bc(loadcase)%rotation)) math_transpose3x3(math_rotate_backward3x3(defgradAim,bc(loadcase)%rotation))
print '(a,x,es10.4)' , 'Determinant of New Deformation Aim:', math_det3x3(defgradAim) print '(a,x,es10.4)' , 'with determinant: ', math_det3x3(defgradAim)
endif endif
print '(A)', '... Calculating Equilibrium Using Spectral Method ...........'
print '(a)', ''
print '(a)', '... calculating equilibrium with spectral method ............'
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
call dfftw_execute_dft_r2c(fftw_plan(1),workfft,workfft) ! FFT of pstress call dfftw_execute_dft_r2c(fftw_plan(1),workfft,workfft) ! FFT of pstress
@ -786,10 +781,15 @@ program DAMASK_spectral
xi(1:3,i,j,k))& xi(1:3,i,j,k))&
)**2.0_pReal)))) )**2.0_pReal))))
enddo; enddo; enddo enddo; enddo; enddo
correctionFactor = minval(geomdimension)*wgt**(-1.0_pReal/4.0_pReal) ! multiplying by minimum dimension to get rid of dimension dependency and phenomenologigal factor wgt**(-1/4) to get rid of resolution dependency if (divergence_correction) then
if (res(3)==1_pInt) correctionFactor = minval(geomdimension(1:2))*wgt**(-1.0_pReal/4.0_pReal) ! 2D case, ToDo: correct? if (res(3) == 1_pInt) then
if (.not. divergence_correction) correctionFactor = 1.0_pReal correctionFactor = minval(geomdimension(1:2))*wgt**(-1.0_pReal/4.0_pReal) ! 2D case, ToDo: correct? PE: Do we need this in the loop or can be pre-calculated?
else
correctionFactor = minval(geomdimension(1:3))*wgt**(-1.0_pReal/4.0_pReal) ! multiplying by minimum dimension to get rid of dimension dependency and phenomenologigal factor wgt**(-1/4) to get rid of resolution dependency
endif
else
correctionFactor = 1.0_pReal
endif
err_div = err_div*wgt/p_hat_avg*correctionFactor ! weighting by points and average stress and multiplying with correction factor err_div = err_div*wgt/p_hat_avg*correctionFactor ! weighting by points and average stress and multiplying with correction factor
err_div_max = err_div_max/p_hat_avg*correctionFactor ! weighting by average stress and multiplying with correction factor err_div_max = err_div_max/p_hat_avg*correctionFactor ! weighting by average stress and multiplying with correction factor
@ -857,21 +857,21 @@ program DAMASK_spectral
call dfftw_execute_dft_c2r(fftw_plan(2),workfft,workfft) call dfftw_execute_dft_c2r(fftw_plan(2),workfft,workfft)
defgrad = defgrad + workfft(1:res(1),:,:,:,:)*wgt defgrad = defgrad + workfft(1:res(1),:,:,:,:)*wgt
do m = 1,3; do n = 1,3 do m = 1_pInt,3_pInt; do n = 1_pInt,3_pInt
defgrad_av(m,n) = sum(defgrad(:,:,:,m,n))*wgt defgrad_av(m,n) = sum(defgrad(:,:,:,m,n))*wgt
enddo; enddo enddo; enddo
defgradAim_lab = math_rotate_backward3x3(defgradAim,bc(loadcase)%rotation) defgradAim_lab = math_rotate_backward3x3(defgradAim,bc(loadcase)%rotation)
do m = 1,3; do n = 1,3 do m = 1_pInt,3_pInt; do n = 1_pInt,3_pInt
defgrad(:,:,:,m,n) = defgrad(:,:,:,m,n) + (defgradAim_lab(m,n) - defgrad_av(m,n)) ! anticipated target minus current state defgrad(:,:,:,m,n) = defgrad(:,:,:,m,n) + (defgradAim_lab(m,n) - defgrad_av(m,n)) ! anticipated target minus current state
enddo; enddo enddo; enddo
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
if(.not. debugDivergence) then if(debugDivergence) then
print '(2(a,es10.4))', 'Error Divergence = ',err_div, ', Tol. = ', err_div_tol print '(a,es10.4,a,f6.2)', 'error divergence FT avg = ',err_div, ', ', err_div/err_div_tol
print '(a,es10.4)', 'error divergence FT max = ',err_div_max
print '(a,es10.4)', 'error divergence Real avg = ',err_real_div_avg
print '(a,es10.4)', 'error divergence Real max = ',err_real_div_max
else else
print '(2(a,es10.4))', 'Error Divergence FT avg= ',err_div, ', Tol. = ', err_div_tol print '(a,es10.4,a,f6.2)', 'error divergence = ',err_div, ', ', err_div/err_div_tol
print '(a,es10.4)', 'Error Divergence FT max= ',err_div_max
print '(a,es10.4)', 'Error Divergence Real avg= ',err_real_div_avg
print '(a,es10.4)', 'Error Divergence Real max= ',err_real_div_max
endif endif
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
@ -880,16 +880,18 @@ program DAMASK_spectral
c_prev = math_rotate_forward3x3x3x3(c_current*wgt,bc(loadcase)%rotation) ! calculate stiffness for next step c_prev = math_rotate_forward3x3x3x3(c_current*wgt,bc(loadcase)%rotation) ! calculate stiffness for next step
!ToDo: Incfluence for next loadcase !ToDo: Incfluence for next loadcase
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
print '(a)', ''
print '(a)', '=============================================================' print '(a)', '============================================================='
if(err_div<=err_div_tol .and. err_stress<=err_stress_tol) then if(err_div > err_div_tol .or. err_stress > err_stress_tol) then
print '(A,I5.5,A)', 'Increment ', totalStepsCounter, ' Converged' print '(A,I5.5,A)', 'increment ', totalStepsCounter, ' NOT converged'
notConvergedCounter = notConvergedCounter + 1_pInt
else else
print '(A,I5.5,A)', 'Increment ', totalStepsCounter, ' NOT Converged' print '(A,I5.5,A)', 'increment ', totalStepsCounter, ' converged'
notConvergedCounter = notConvergedCounter + 1
endif endif
if (mod(totalStepsCounter -1_pInt,bc(loadcase)%outputfrequency) == 0_pInt) then ! at output frequency if (mod(totalStepsCounter -1_pInt,bc(loadcase)%outputfrequency) == 0_pInt) then ! at output frequency
print '(A)', '... Writing Results to File .................................' print '(a)', ''
write(538), materialpoint_results(materialpoint_sizeResults,1,res(1)*res(2)*res(3)) ! write result to file print '(a)', '... writing results to file .................................'
write(538), materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! write result to file
endif endif
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
@ -898,8 +900,9 @@ program DAMASK_spectral
deallocate(s_reduced) deallocate(s_reduced)
enddo ! end looping over loadcases enddo ! end looping over loadcases
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
print '(a)', ''
print '(a)', '#############################################################' print '(a)', '#############################################################'
print '(a,i5.5,a,i5.5,a)', 'Of ', totalStepsCounter - restartReadStep + 1_pInt, ' Calculated Steps, ', notConvergedCounter, ' Steps did not Converge!' print '(a,i5.5,a,i5.5,a)', 'of ', totalStepsCounter - restartReadStep + 1_pInt, ' calculated steps, ', notConvergedCounter, ' steps did not converge!'
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
close(538) close(538)
call dfftw_destroy_plan(fftw_plan(1)); call dfftw_destroy_plan(fftw_plan(2)) call dfftw_destroy_plan(fftw_plan(1)); call dfftw_destroy_plan(fftw_plan(2))

21
code/FEsolving.f90
View File

@ -62,10 +62,10 @@
FEmodelGeometry = getModelName() FEmodelGeometry = getModelName()
if (IO_open_inputFile(fileunit,FEmodelGeometry)) then if (IO_open_inputFile(fileunit,FEmodelGeometry)) then
if (trim(FEsolver) == 'Spectral') then if (trim(FEsolver) == 'Spectral') then
restartWrite = .true. call get_command(commandLine) ! may contain uppercase
call get_command(commandLine) ! may contain capitals do i=1,len(commandLine)
do i=1,len(commandLine) ! remove capitals if(64 < iachar(commandLine(i:i)) .and. iachar(commandLine(i:i)) < 91) &
if(64<iachar(commandLine(i:i)) .and. iachar(commandLine(i:i))<91) commandLine(i:i) =achar(iachar(commandLine(i:i))+32) commandLine(i:i) = achar(iachar(commandLine(i:i))+32) ! make lowercase
enddo enddo
start = index(commandLine,'-r ',.true.) + 3_pInt ! search for '-r' and jump forward to given name start = index(commandLine,'-r ',.true.) + 3_pInt ! search for '-r' and jump forward to given name
if (index(commandLine,'--restart ',.true.)>0) then ! if '--restart' is found, use that (contains '-r') if (index(commandLine,'--restart ',.true.)>0) then ! if '--restart' is found, use that (contains '-r')
@ -74,9 +74,16 @@
length = index(commandLine(start:len(commandLine)),' ',.false.) length = index(commandLine(start:len(commandLine)),' ',.false.)
if(start /= 3_pInt) then ! found at least -r if(start /= 3_pInt) then ! found at least -r
read(commandLine(start:start+length),'(I12)') restartReadStep read(commandLine(start:start+length),'(I12)') restartReadStep
restartRead = .true.
endif endif
if(restartReadStep<0_pInt .and. restartRead .eq. .true.) call IO_error(error_ID=47) if (restartReadStep > 0_pInt) then
restartRead = .true.
restartWrite = .true.
endif
if (restartReadStep == 0_pInt) then
restartRead = .false.
restartWrite = .false.
endif
if(restartReadStep < 0_pInt) call IO_error(error_ID=47)
else else
rewind(fileunit) rewind(fileunit)
do do
@ -107,7 +114,7 @@
enddo enddo
endif endif
else else
call IO_error(101) ! cannot open input file call IO_error(101, ext_msg=FEmodelGeometry) ! cannot open input file
endif endif
100 close(fileunit) 100 close(fileunit)