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

made DAMASK to work with gfortran: 

-removed to long lines
-restructured f2py modules and merged make_DAMASK2Python into setup processing
-setup_code.py now sets library path in makefile and asks for compile switches for spectral code
-substituted \ in format strings with $

restructured DAMASK_spectral:
-more logical output and structure of code
-better input for spectral debug parameters
This commit is contained in:
Martin Diehl 2011-12-06 16:58:17 +00:00
parent b752becbea
commit 5ebeb96e85
15 changed files with 458 additions and 368 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

0
processing/post/DAMASK.pyf → code/DAMASK.pyf
View File

470
code/DAMASK_spectral.f90
View File

@ -52,7 +52,7 @@ program DAMASK_spectral
use math use math
use mesh, only: mesh_ipCenterOfGravity use mesh, only: mesh_ipCenterOfGravity
use CPFEM, only: CPFEM_general, CPFEM_initAll use CPFEM, only: CPFEM_general, CPFEM_initAll
use FEsolving, only: restartWrite, restartReadSpectral, restartReadStep use FEsolving, only: restartWrite, restartReadStep
use numerics, only: err_div_tol, err_stress_tolrel , rotation_tol,& use numerics, only: err_div_tol, err_stress_tolrel , rotation_tol,&
itmax, memory_efficient, DAMASK_NumThreadsInt, divergence_correction, & itmax, memory_efficient, DAMASK_NumThreadsInt, divergence_correction, &
fftw_planner_flag, fftw_timelimit fftw_planner_flag, fftw_timelimit
@ -103,10 +103,10 @@ program DAMASK_spectral
logical :: spectralPictureMode = .false. ! indicating 1 to 1 mapping of FP to microstructure logical :: spectralPictureMode = .false. ! indicating 1 to 1 mapping of FP to microstructure
! stress, stiffness and compliance average etc. ! stress, stiffness and compliance average etc.
real(pReal), dimension(3,3) :: pstress, pstress_av, defgrad_av, & real(pReal), dimension(3,3) :: pstress, pstress_av, defgrad_av_lab, &
defgradAim = math_I3, defgradAimOld= math_I3, defgradAimCorr= math_I3,& defgradAim = math_I3, defgradAimOld= math_I3, defgradAimCorr= math_I3,&
mask_stress, mask_defgrad, fDot, & mask_stress, mask_defgrad, fDot, &
pstress_av_load, defgradAim_lab ! quantities rotated to other coordinate system pstress_av_lab, defgradAim_lab ! quantities rotated to other coordinate system
real(pReal), dimension(3,3,3,3) :: dPdF, c0_reference, c_current = 0.0_pReal, s_prev, c_prev ! stiffness and compliance real(pReal), dimension(3,3,3,3) :: dPdF, c0_reference, c_current = 0.0_pReal, s_prev, c_prev ! stiffness and compliance
real(pReal), dimension(6) :: cstress ! cauchy stress real(pReal), dimension(6) :: cstress ! cauchy stress
real(pReal), dimension(6,6) :: dsde ! small strain stiffness real(pReal), dimension(6,6) :: dsde ! small strain stiffness
@ -149,7 +149,7 @@ program DAMASK_spectral
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
bc_init%outputfrequency = 1_pInt; bc_init%restartfrequency = 1_pInt bc_init%outputfrequency = 1_pInt; bc_init%restartfrequency = 0_pInt
bc_init%maskDeformation = .false.; bc_init%maskStress = .false. bc_init%maskDeformation = .false.; bc_init%maskStress = .false.
bc_init%maskStressVector = .false.; bc_init%velGradApplied = .false. bc_init%maskStressVector = .false.; bc_init%velGradApplied = .false.
bc_init%followFormerTrajectory = .true. bc_init%followFormerTrajectory = .true.
@ -242,7 +242,7 @@ program DAMASK_spectral
case('f','freq','frequency','outputfreq') ! frequency of result writings case('f','freq','frequency','outputfreq') ! frequency of result writings
bc(loadcase)%outputfrequency = IO_intValue(line,posLoadcase,j+1_pInt) bc(loadcase)%outputfrequency = IO_intValue(line,posLoadcase,j+1_pInt)
case('r','restart','restartwrite') ! frequency of writing restart information case('r','restart','restartwrite') ! frequency of writing restart information
bc(loadcase)%restartfrequency = IO_intValue(line,posLoadcase,j+1_pInt) bc(loadcase)%restartfrequency = max(0_pInt,IO_intValue(line,posLoadcase,j+1_pInt))
case('guessreset','dropguessing') case('guessreset','dropguessing')
bc(loadcase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory bc(loadcase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory
case('euler') ! rotation of loadcase given in euler angles case('euler') ! rotation of loadcase given in euler angles
@ -367,27 +367,28 @@ program DAMASK_spectral
if (.not. bc(loadcase)%followFormerTrajectory) print '(a)', 'drop guessing along trajectory' if (.not. bc(loadcase)%followFormerTrajectory) print '(a)', 'drop guessing along trajectory'
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) == .true.) .and. & if (any(bc(loadcase)%maskDeformation(j,1:3) .eqv. .true.) .and. &
any(bc(loadcase)%maskDeformation(j,1:3) == .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 = 32_pInt ! each row should be either fully or not at all defined
enddo enddo
print '(a)','velocity gradient:' print '(a)','velocity gradient:'
else else
print '(a)','deformation gradient rate:' print '(a)','deformation gradient rate:'
endif endif
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 / GPa:',1e-9*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)) transpose(bc(loadcase)%maskStress))
if (any(bc(loadcase)%rotation /= math_I3)) & if (any(bc(loadcase)%rotation /= math_I3)) &
print '(a,3(3(f12.6,x)/)\)','rotation of loadframe:',math_transpose3x3(bc(loadcase)%rotation) print '(a,3(3(f12.6,x)/)$)','rotation of loadframe:',math_transpose3x3(bc(loadcase)%rotation)
print '(a,f12.6)','temperature:',bc(loadcase)%temperature print '(a,f12.6)','temperature:',bc(loadcase)%temperature
print '(a,f12.6)','time: ',bc(loadcase)%timeIncrement print '(a,f12.6)','time: ',bc(loadcase)%timeIncrement
print '(a,i5)','steps: ',bc(loadcase)%steps print '(a,i5)','steps: ',bc(loadcase)%steps
print '(a,i5)','output frequency: ',bc(loadcase)%outputfrequency print '(a,i5)','output frequency: ',bc(loadcase)%outputfrequency
print '(a,i5)','restart frequency: ',bc(loadcase)%restartfrequency print '(a,i5)','restart frequency: ',bc(loadcase)%restartfrequency
if (any(bc(loadcase)%maskStress == bc(loadcase)%maskDeformation)) errorID = 31 ! exclusive or masking only
if (any(bc(loadcase)%maskStress .eqv. bc(loadcase)%maskDeformation)) errorID = 31 ! exclusive or masking only
if (any(bc(loadcase)%maskStress .and. transpose(bc(loadcase)%maskStress) .and. & if (any(bc(loadcase)%maskStress .and. transpose(bc(loadcase)%maskStress) .and. &
reshape((/.false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false./),(/3,3/)))) & reshape((/.false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false./),(/3,3/)))) &
errorID = 38_pInt ! no rotation is allowed by stress BC errorID = 38_pInt ! no rotation is allowed by stress BC
@ -398,8 +399,6 @@ program DAMASK_spectral
if (bc(loadcase)%timeIncrement < 0.0_pReal) errorID = 34_pInt ! negative time increment if (bc(loadcase)%timeIncrement < 0.0_pReal) errorID = 34_pInt ! negative time increment
if (bc(loadcase)%steps < 1_pInt) errorID = 35_pInt ! non-positive increment count if (bc(loadcase)%steps < 1_pInt) errorID = 35_pInt ! non-positive increment count
if (bc(loadcase)%outputfrequency < 1_pInt) errorID = 36_pInt ! non-positive result frequency if (bc(loadcase)%outputfrequency < 1_pInt) errorID = 36_pInt ! non-positive result frequency
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) if (errorID > 0_pInt) call IO_error(error_ID = errorID, ext_msg = loadcase_string)
enddo enddo
@ -455,9 +454,9 @@ program DAMASK_spectral
!************************************************************* !*************************************************************
! Initialization Start ! Initialization Start
!************************************************************* !*************************************************************
if(totalStepsCounter >= restartReadStep) then ! Do calculations (otherwise just forwarding)
if(totalStepsCounter == restartReadStep) then ! Initialize values
if (regrid==.true. ) then ! 'DeInitialize' the values changing in case of regridding if (regrid .eqv. .true. ) then ! 'DeInitialize' the values changing in case of regridding
regrid = .false. regrid = .false.
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))
if(debugDivergence) call dfftw_destroy_plan(fftw_plan(3)) if(debugDivergence) call dfftw_destroy_plan(fftw_plan(3))
@ -467,147 +466,153 @@ program DAMASK_spectral
deallocate (temperature) deallocate (temperature)
deallocate (xi) deallocate (xi)
deallocate (workfft) deallocate (workfft)
!ToDo: here we have to create the new geometry and assign the values from the previous step !ToDo: here we have to create the new geometry and assign the values from the previous step
endif endif
if(totalStepsCounter == restartReadStep) then ! Initialize values 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 allocate (defgrad ( res(1),res(2),res(3),3,3)); defgrad = 0.0_pReal
allocate (defgrad ( res(1),res(2),res(3),3,3)); defgrad = 0.0_pReal allocate (defgradold ( res(1),res(2),res(3),3,3)); defgradold = 0.0_pReal
allocate (defgradold ( res(1),res(2),res(3),3,3)); defgradold = 0.0_pReal allocate (coordinates(3,res(1),res(2),res(3))); coordinates = 0.0_pReal
allocate (coordinates(3,res(1),res(2),res(3))); coordinates = 0.0_pReal allocate (temperature( res(1),res(2),res(3))); temperature = bc(1)%temperature ! start out isothermally
allocate (temperature( res(1),res(2),res(3))); temperature = bc(1)%temperature ! start out isothermally allocate (xi (3,res(1)/2+1,res(2),res(3))); xi =0.0_pReal
allocate (xi (3,res(1)/2+1,res(2),res(3))); xi =0.0_pReal 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(Npoints, 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)
call dfftw_plan_many_dft_c2r(fftw_plan(2),3,(/res(1),res(2),res(3)/),9,& call dfftw_plan_many_dft_c2r(fftw_plan(2),3,(/res(1),res(2),res(3)/),9,&
workfft,(/res(1)/2_pInt+1_pInt,res(2),res(3)/),1,(res(1)/2_pInt+1_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),&
workfft,(/res(1) +2_pInt,res(2),res(3)/),1,(res(1) +2_pInt)*res(2)*res(3),fftw_planner_flag) workfft,(/res(1) +2_pInt,res(2),res(3)/),1,(res(1) +2_pInt)*res(2)*res(3),fftw_planner_flag)
if (debugDivergence) & if (debugDivergence) &
call dfftw_plan_many_dft_c2r(fftw_plan(3),3,(/res(1),res(2),res(3)/),3,& call dfftw_plan_many_dft_c2r(fftw_plan(3),3,(/res(1),res(2),res(3)/),3,&
divergence,(/res(1)/2_pInt+1_pInt,res(2),res(3)/),1,(res(1)/2_pInt+1_pInt)*res(2)*res(3),& divergence,(/res(1)/2_pInt+1_pInt,res(2),res(3)/),1,(res(1)/2_pInt+1_pInt)*res(2)*res(3),&
divergence,(/res(1) +2_pInt,res(2),res(3)/),1,(res(1) +2_pInt)*res(2)*res(3),fftw_planner_flag) divergence,(/res(1) +2_pInt,res(2),res(3)/),1,(res(1) +2_pInt)*res(2)*res(3),fftw_planner_flag)
if (debugGeneral) then if (debugGeneral) then
!$OMP CRITICAL (write2out)
write (6,*) 'FFTW initialized'
!$OMP END CRITICAL (write2out)
endif
if (restartReadStep==1_pInt) then ! not restarting, no deformation at the beginning
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
defgrad(i,j,k,1:3,1:3) = math_I3
defgradold(i,j,k,1:3,1:3) = math_I3
enddo; enddo; enddo
else ! using old values
if (IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',trim(getSolverJobName()),size(defgrad))) then
read (777,rec=1) defgrad
close (777)
endif
defgradold = defgrad
defgradAim = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
defgradAim = defgradAim + defgrad(i,j,k,1:3,1:3) ! calculating old average deformation
enddo; enddo; enddo
defgradAim = defgradAim * wgt
defgradAimOld = defgradAim
guessmode=0.0_pInt
endif
ielem = 0_pInt
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
ielem = ielem + 1_pInt
coordinates(1:3,i,j,k) = mesh_ipCenterOfGravity(1:3,1,ielem) ! set to initial coordinates ToDo: SHOULD BE UPDATED TO CURRENT POSITION IN FUTURE REVISIONS!!! But do we know them? I don't think so. Otherwise we don't need geometry reconstruction
call CPFEM_general(2_pInt,coordinates(1:3,i,j,k),math_I3,math_I3,temperature(i,j,k),0.0_pReal,ielem,1_pInt,cstress,dsde,pstress,dPdF)
c_current = c_current + dPdF
enddo; enddo; enddo
c0_reference = c_current * wgt ! linear reference material stiffness
c_prev = math_rotate_forward3x3x3x3(c0_reference,bc(loadcase)%rotation) ! rotate_forward: lab -> load system
if (debugGeneral) then
!$OMP CRITICAL (write2out)
write (6,*) 'first call to CPFEM_general finished'
!$OMP END CRITICAL (write2out)
endif
do k = 1_pInt, res(3) ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around)
k_s(3) = k - 1_pInt
if(k > res(3)/2_pInt + 1_pInt) k_s(3) = k_s(3) - res(3)
do j = 1_pInt, res(2)
k_s(2) = j - 1_pInt
if(j > res(2)/2_pInt + 1_pInt) k_s(2) = k_s(2) - res(2)
do i = 1, res(1)/2_pInt + 1_pInt
k_s(1) = i - 1_pInt
xi(3,i,j,k) = 0.0_pReal ! 2D case
if(res(3) > 1_pInt) xi(3,i,j,k) = real(k_s(3), pReal)/geomdimension(3) ! 3D case
xi(2,i,j,k) = real(k_s(2), pReal)/geomdimension(2)
xi(1,i,j,k) = real(k_s(1), pReal)/geomdimension(1)
enddo; enddo; enddo
! 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
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(i==res(1)/2_pInt+1_pInt) xi(1,i,j,k)= 0.0_pReal
enddo; enddo; enddo
if(memory_efficient) then ! allocate just single fourth order tensor
allocate (gamma_hat(1,1,1,3,3,3,3)); gamma_hat = 0.0_pReal
else ! precalculation of gamma_hat field
allocate (gamma_hat(res(1)/2_pInt + 1_pInt ,res(2),res(3),3,3,3,3)); gamma_hat = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)/2_pInt + 1_pInt
if (any(xi(:,i,j,k) /= 0.0_pReal)) then
do l = 1_pInt ,3_pInt; do m = 1_pInt,3_pInt
xiDyad(l,m) = xi(l,i,j,k)*xi(m,i,j,k)
enddo; enddo
temp33_Real = math_inv3x3(math_mul3333xx33(c0_reference, xiDyad))
else
xiDyad = 0.0_pReal
temp33_Real = 0.0_pReal
endif
do l=1_pInt,3_pInt; do m=1_pInt,3_pInt; do n=1_pInt,3_pInt; do p=1_pInt,3_pInt
gamma_hat(i,j,k, l,m,n,p) = - 0.25*(temp33_Real(l,n)+temp33_Real(n,l)) *&
(xiDyad(m,p)+xiDyad(p,m))
enddo; enddo; enddo; enddo
enddo; enddo; enddo
endif
! write header of output file
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
open(538,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())& write (6,*) 'FFTW initialized'
//'.spectralOut',form='UNFORMATTED',status='REPLACE') !$OMP END CRITICAL (write2out)
write(538), 'load', trim(getLoadcaseName())
write(538), 'workingdir', trim(getSolverWorkingDirectoryName())
write(538), 'geometry', trim(getSolverJobName())//InputFileExtension
write(538), 'resolution', res
write(538), 'dimension', geomdimension
write(538), 'materialpoint_sizeResults', materialpoint_sizeResults
write(538), 'loadcases', N_Loadcases
write(538), 'logscale', bc(1:N_Loadcases)%logscale ! one entry per loadcase (0: linear, 1: log)
write(538), 'frequencies', bc(1:N_Loadcases)%outputfrequency ! one entry per loadcase
write(538), 'times', bc(1:N_Loadcases)%timeIncrement ! one entry per loadcase
bc(1)%steps= bc(1)%steps + 1_pInt
write(538), 'increments', bc(1:N_Loadcases)%steps ! one entry per loadcase ToDo: rename keyword to steps
bc(1)%steps= bc(1)%steps - 1_pInt
write(538), 'startingIncrement', restartReadStep -1_pInt ! start with writing out the previous step
write(538), 'eoh' ! end of header
write(538), materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! initial (non-deformed) results
!$OMP END CRITICAL (write2out)
endif endif
if (restartReadStep==1_pInt) then ! not restarting, no deformation at the beginning
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
defgrad(i,j,k,1:3,1:3) = math_I3
defgradold(i,j,k,1:3,1:3) = math_I3
enddo; enddo; enddo
else ! using old values
if (IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',trim(getSolverJobName()),size(defgrad))) then
read (777,rec=1) defgrad
close (777)
endif
defgradold = defgrad
defgradAim = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
defgradAim = defgradAim + defgrad(i,j,k,1:3,1:3) ! calculating old average deformation
enddo; enddo; enddo
defgradAim = defgradAim * wgt
defgradAimOld = defgradAim
guessmode=0.0_pInt
endif
ielem = 0_pInt
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
ielem = ielem + 1_pInt
coordinates(1:3,i,j,k) = mesh_ipCenterOfGravity(1:3,1,ielem) ! set to initial coordinates ToDo: SHOULD BE UPDATED TO CURRENT POSITION IN FUTURE REVISIONS!!! But do we know them? I don't think so. Otherwise we don't need geometry reconstruction
call CPFEM_general(2_pInt,coordinates(1:3,i,j,k),math_I3,math_I3,temperature(i,j,k),&
0.0_pReal,ielem,1_pInt,cstress,dsde,pstress,dPdF)
c_current = c_current + dPdF
enddo; enddo; enddo
c0_reference = c_current * wgt ! linear reference material stiffness
if (debugGeneral) then
!$OMP CRITICAL (write2out)
write (6,*) 'first call to CPFEM_general finished'
!$OMP END CRITICAL (write2out)
endif
do k = 1_pInt, res(3) ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around)
k_s(3) = k - 1_pInt
if(k > res(3)/2_pInt + 1_pInt) k_s(3) = k_s(3) - res(3)
do j = 1_pInt, res(2)
k_s(2) = j - 1_pInt
if(j > res(2)/2_pInt + 1_pInt) k_s(2) = k_s(2) - res(2)
do i = 1, res(1)/2_pInt + 1_pInt
k_s(1) = i - 1_pInt
xi(3,i,j,k) = 0.0_pReal ! 2D case
if(res(3) > 1_pInt) xi(3,i,j,k) = real(k_s(3), pReal)/geomdimension(3) ! 3D case
xi(2,i,j,k) = real(k_s(2), pReal)/geomdimension(2)
xi(1,i,j,k) = real(k_s(1), pReal)/geomdimension(1)
enddo; enddo; enddo
!remove highest frequencies for calculation of divergence (CAREFULL, 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
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(i==res(1)/2_pInt+1_pInt) xi(1,i,j,k)= 0.0_pReal
enddo; enddo; enddo
if(memory_efficient) then ! allocate just single fourth order tensor
allocate (gamma_hat(1,1,1,3,3,3,3)); gamma_hat = 0.0_pReal
else ! precalculation of gamma_hat field
allocate (gamma_hat(res(1)/2_pInt + 1_pInt ,res(2),res(3),3,3,3,3)); gamma_hat = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)/2_pInt + 1_pInt
if (any(xi(:,i,j,k) /= 0.0_pReal)) then
do l = 1_pInt ,3_pInt; do m = 1_pInt,3_pInt
xiDyad(l,m) = xi(l,i,j,k)*xi(m,i,j,k)
enddo; enddo
temp33_Real = math_inv3x3(math_mul3333xx33(c0_reference, xiDyad))
else
xiDyad = 0.0_pReal
temp33_Real = 0.0_pReal
endif
do l=1_pInt,3_pInt; do m=1_pInt,3_pInt; do n=1_pInt,3_pInt; do p=1_pInt,3_pInt
gamma_hat(i,j,k, l,m,n,p) = - 0.25*(temp33_Real(l,n)+temp33_Real(n,l)) *&
(xiDyad(m,p)+xiDyad(p,m))
enddo; enddo; enddo; enddo
enddo; enddo; enddo
endif
if (divergence_correction) then
if (res(3) == 1_pInt) then
correctionFactor = minval(geomdimension(1:2))*wgt**(-1.0_pReal/4.0_pReal) ! 2D case, ToDo: correct?
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
! write header of output file
!$OMP CRITICAL (write2out)
open(538,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())&
//'.spectralOut',form='UNFORMATTED',status='REPLACE')
write(538), 'load', trim(getLoadcaseName())
write(538), 'workingdir', trim(getSolverWorkingDirectoryName())
write(538), 'geometry', trim(getSolverJobName())//InputFileExtension
write(538), 'resolution', res
write(538), 'dimension', geomdimension
write(538), 'materialpoint_sizeResults', materialpoint_sizeResults
write(538), 'loadcases', N_Loadcases
write(538), 'logscale', bc(1:N_Loadcases)%logscale ! one entry per loadcase (0: linear, 1: log)
write(538), 'frequencies', bc(1:N_Loadcases)%outputfrequency ! one entry per loadcase
write(538), 'times', bc(1:N_Loadcases)%timeIncrement ! one entry per loadcase
bc(1)%steps= bc(1)%steps + 1_pInt
write(538), 'increments', bc(1:N_Loadcases)%steps ! one entry per loadcase ToDo: rename keyword to steps
bc(1)%steps= bc(1)%steps - 1_pInt
write(538), 'startingIncrement', restartReadStep -1_pInt ! start with writing out the previous step
write(538), 'eoh' ! end of header
write(538), materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! initial (non-deformed) results
!$OMP END CRITICAL (write2out)
endif
!************************************************************* !*************************************************************
! Initialization End ! Initialization End
!************************************************************* !*************************************************************
if (mod(step - 1_pInt,bc(loadcase)%restartFrequency)==0_pInt) then ! at frequency of writing restart information if(totalStepsCounter >= restartReadStep) then ! Do calculations (otherwise just forwarding)
restartWrite = .true. ! setting restart parameter for FEsolving (first call to CPFEM_general will write ToDo: true?) if(bc(loadcase)%restartFrequency>0_pInt) &
else restartWrite = ( mod(step - 1_pInt,bc(loadcase)%restartFrequency)==0_pInt) ! at frequency of writing restart information
restartWrite = .false. ! setting restart parameter for FEsolving (first call to CPFEM_general will write ToDo: true?)
endif if (bc(loadcase)%velGradApplied) & ! calculate fDot from given L and current F
fDot = math_mul33x33(bc(loadcase)%deformation, defgradAim)
if (bc(loadcase)%velGradApplied) & ! calculate fDot from given L and current F
fDot = math_mul33x33(bc(loadcase)%deformation, defgradAim)
!winding forward of deformation aim in loadcase system !winding forward of deformation aim in loadcase system
temp33_Real = defgradAim temp33_Real = defgradAim
@ -645,7 +650,8 @@ program DAMASK_spectral
CPFEM_mode = 1_pInt ! winding forward CPFEM_mode = 1_pInt ! winding forward
iter = 0_pInt iter = 0_pInt
err_div = 2.0_pReal * err_div_tol ! go into loop err_div = 2.0_pReal * err_div_tol ! go into loop
c_prev = math_rotate_forward3x3x3x3(c_current*wgt,bc(loadcase)%rotation) ! calculate stiffness from former step
if(size_reduced > 0_pInt) then ! calculate compliance in case stress BC is applied if(size_reduced > 0_pInt) then ! calculate compliance in case stress BC is applied
c_prev99 = math_Plain3333to99(c_prev) c_prev99 = math_Plain3333to99(c_prev)
k = 0_pInt ! build reduced stiffness k = 0_pInt ! build reduced stiffness
@ -674,17 +680,17 @@ program DAMASK_spectral
s_prev = (math_Plain99to3333(s_prev99)) s_prev = (math_Plain99to3333(s_prev99))
endif endif
!$OMP CRITICAL (write2out)
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)
endif endif
endif endif
!$OMP END CRITICAL (write2out)
!************************************************************* !*************************************************************
! convergence loop ! convergence loop
do while(iter < itmax .and. & do while(iter < itmax .and. &
@ -692,34 +698,27 @@ 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 '(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_lab(m,n) = sum(defgrad(1:res(1),1:res(2),1:res(3),m,n)) * wgt
enddo; enddo enddo; enddo
!$OMP CRITICAL (write2out) print '(a,/,3(3(f12.7,x)/)$)', 'deformation gradient:',&
print '(a,/,3(3(f12.7,x)/)\)', 'deformation gradient:',math_transpose3x3(defgrad_av) math_transpose3x3(math_rotate_forward3x3(defgrad_av_lab,bc(loadcase)%rotation))
print '(l)', restartWrite print '(a)', ''
print '(a)', '... update stress field P(F) ................................' print '(a)', '... update stress field P(F) ................................'
!$OMP END CRITICAL (write2out)
defgradDetMax = -huge(1.0_pReal)
defgradDetMin = +huge(1.0_pReal)
ielem = 0_pInt ielem = 0_pInt
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
defgradDet = math_det3x3(defgrad(i,j,k,1:3,1:3))
defgradDetMax = max(defgradDetMax,defgradDet)
defgradDetMin = min(defgradDetMin,defgradDet)
ielem = ielem + 1_pInt ielem = ielem + 1_pInt
call CPFEM_general(3_pInt,& ! collect cycle call CPFEM_general(3_pInt,& ! collect cycle
coordinates(1:3,i,j,k), defgradold(i,j,k,1:3,1:3), defgrad(i,j,k,1:3,1:3),& coordinates(1:3,i,j,k), defgradold(i,j,k,1:3,1:3), defgrad(i,j,k,1:3,1:3),&
temperature(i,j,k),timeinc,ielem,1_pInt,& temperature(i,j,k),timeinc,ielem,1_pInt,&
cstress,dsde, pstress, dPdF) cstress,dsde, pstress, dPdF)
enddo; enddo; enddo enddo; enddo; enddo
print '(a,x,es10.4)' , 'max determinant of deformation:', defgradDetMax
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
ielem = 0_pInt ielem = 0_pInt
@ -731,43 +730,24 @@ program DAMASK_spectral
temperature(i,j,k),timeinc,ielem,1_pInt,& temperature(i,j,k),timeinc,ielem,1_pInt,&
cstress,dsde, pstress,dPdF) cstress,dsde, pstress,dPdF)
CPFEM_mode = 2_pInt CPFEM_mode = 2_pInt
workfft(i,j,k,1:3,1:3) = pstress ! build up average P-K stress
workfft(i,j,k,1:3,1:3) = pstress ! build up average P-K stress
c_current = c_current + dPdF c_current = c_current + dPdF
enddo; enddo; enddo enddo; enddo; enddo
restartWrite = .false. ! ToDo: don't know if we need it. Depends on how CPFEM_general is writing results
do n = 1_pInt,3_pInt; do m = 1_pInt,3_pInt do n = 1_pInt,3_pInt; do m = 1_pInt,3_pInt
pstress_av(m,n) = sum(workfft(1:res(1),1:res(2),1:res(3),m,n)) * wgt pstress_av_lab(m,n) = sum(workfft(1:res(1),1:res(2),1:res(3),m,n)) * wgt
enddo; enddo enddo; enddo
!$OMP CRITICAL (write2out)
print '(a,/,3(3(f12.7,x)/)\)', 'Piola-Kirchhoff stress / MPa: ',math_transpose3x3(pstress_av)/1.e6
err_stress_tol = 0.0_pReal
pstress_av_load = math_rotate_forward3x3(pstress_av,bc(loadcase)%rotation)
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_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)', ''
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
defgradAim = defgradAim + defgradAimCorr
print '(a,/,3(3(f12.7,x)/)\)', 'new deformation aim: ',&
math_transpose3x3(math_rotate_backward3x3(defgradAim,bc(loadcase)%rotation))
print '(a,x,es10.4)' , 'with determinant: ', math_det3x3(defgradAim)
endif
print '(a)', '' print '(a)', ''
print '(a)', '... calculating equilibrium with spectral method ............' print '(a)', '... calculating equilibrium with spectral method ............'
!$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
p_hat_avg = sqrt(maxval (math_eigenvalues3x3(math_mul33x33(workfft(1,1,1,1:3,1:3),& ! L_2 norm of average stress in fourier space, p_hat_avg = sqrt(maxval (math_eigenvalues3x3(math_mul33x33(workfft(1,1,1,1:3,1:3),& ! L_2 norm of average stress in fourier space,
math_transpose3x3(workfft(1,1,1,1:3,1:3)))))) ! ignore imaginary part as it is always zero for real only input)) math_transpose3x3(workfft(1,1,1,1:3,1:3)))))) ! ignore imaginary part as it is always zero for real only input))
err_div = 0.0_pReal err_div = 0.0_pReal
err_div_max = 0.0_pReal err_div_max = 0.0_pReal ! only if if(debugDivergence) == .true. of importace
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
err_div = err_div + sqrt(sum((& ! avg of L_2 norm of div(stress) in fourier space (Suquet small strain) err_div = err_div + sqrt(sum((& ! avg of L_2 norm of div(stress) in fourier space (Suquet small strain)
math_mul33x3_complex(workfft(i*2_pInt-1_pInt,j,k,1:3,1:3) + & math_mul33x3_complex(workfft(i*2_pInt-1_pInt,j,k,1:3,1:3) + &
@ -781,17 +761,9 @@ 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
if (divergence_correction) then
if (res(3) == 1_pInt) then
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, only if if(debugDivergence) == .true. of importace
if(memory_efficient) then ! memory saving version, on-the-fly calculation of gamma_hat if(memory_efficient) then ! memory saving version, on-the-fly calculation of gamma_hat
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
@ -825,6 +797,8 @@ program DAMASK_spectral
workfft(i*2_pInt ,j,k,1:3,1:3) = aimag(temp33_Complex) workfft(i*2_pInt ,j,k,1:3,1:3) = aimag(temp33_Complex)
enddo; enddo; enddo enddo; enddo; enddo
endif endif
! calculate additional divergence criteria and report -------------
if(debugDivergence) then if(debugDivergence) then
divergence = 0.0_pReal divergence = 0.0_pReal
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
@ -845,27 +819,11 @@ program DAMASK_spectral
err_real_div_avg = err_real_div_avg + sqrt(sum((divergence(i,j,k,1:3))**2.0_pReal)) ! avg of L_2 norm of div(stress) in fourier space (Suquet small strain) err_real_div_avg = err_real_div_avg + sqrt(sum((divergence(i,j,k,1:3))**2.0_pReal)) ! avg of L_2 norm of div(stress) in fourier space (Suquet small strain)
err_real_div_max = max(err_real_div_max,abs(sqrt(sum((divergence(i,j,k,1:3))**2.0_pReal)))) ! maximum of L two norm of div(stress) in fourier space (Suquet large strain) err_real_div_max = max(err_real_div_max,abs(sqrt(sum((divergence(i,j,k,1:3))**2.0_pReal)))) ! maximum of L two norm of div(stress) in fourier space (Suquet large strain)
enddo; enddo; enddo enddo; enddo; enddo
p_real_avg = sqrt(maxval (math_eigenvalues3x3(math_mul33x33(pstress_av,& ! L_2 norm of average stress in fourier space, p_real_avg = sqrt(maxval (math_eigenvalues3x3(math_mul33x33(pstress_av_lab,& ! L_2 norm of average stress in real space,
math_transpose3x3(pstress_av))))) ! ignore imaginary part as it is always zero for real only input)) math_transpose3x3(pstress_av_lab))))) ! ignore imaginary part as it is always zero for real only input))
err_real_div_avg = err_real_div_avg*wgt/p_real_avg err_real_div_avg = err_real_div_avg*wgt/p_real_avg
err_real_div_max = err_real_div_max/p_real_avg err_real_div_max = err_real_div_max/p_real_avg
endif
! average strain
workfft(1,1,1,1:3,1:3) = defgrad_av - math_I3 ! zero frequency (real part)
workfft(2,1,1,1:3,1:3) = 0.0_pReal ! zero frequency (imaginary part)
call dfftw_execute_dft_c2r(fftw_plan(2),workfft,workfft)
defgrad = defgrad + workfft(1:res(1),:,:,:,:)*wgt
do m = 1_pInt,3_pInt; do n = 1_pInt,3_pInt
defgrad_av(m,n) = sum(defgrad(:,:,:,m,n))*wgt
enddo; enddo
defgradAim_lab = math_rotate_backward3x3(defgradAim,bc(loadcase)%rotation)
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
enddo; enddo
!$OMP CRITICAL (write2out)
if(debugDivergence) then
print '(a,es10.4,a,f6.2)', 'error divergence FT avg = ',err_div, ', ', err_div/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 FT max = ',err_div_max
print '(a,es10.4)', 'error divergence Real avg = ',err_real_div_avg print '(a,es10.4)', 'error divergence Real avg = ',err_real_div_avg
@ -873,12 +831,63 @@ program DAMASK_spectral
else else
print '(a,es10.4,a,f6.2)', 'error divergence = ',err_div, ', ', err_div/err_div_tol print '(a,es10.4,a,f6.2)', 'error divergence = ',err_div, ', ', err_div/err_div_tol
endif endif
!$OMP END CRITICAL (write2out) ! --------------------------
workfft(1,1,1,1:3,1:3) = defgrad_av_lab - math_I3 ! assign zero frequency (real part) with average displacement gradient
workfft(2,1,1,1:3,1:3) = 0.0_pReal ! zero frequency (imaginary part)
call dfftw_execute_dft_c2r(fftw_plan(2),workfft,workfft) ! backtransformation of fluct deformation gradient
defgrad = defgrad + workfft(1:res(1),1:res(2),1:res(3),1:3,1:3)*wgt ! F(x)^(n+1) = F(x)^(n) + correction; *wgt: correcting for missing normalization
do m = 1_pInt,3_pInt; do n = 1_pInt,3_pInt
defgrad_av_lab(m,n) = sum(defgrad(1:res(1),1:res(2),1:res(3),m,n))*wgt ! ToDo: check whether this needs recalculation or is equivalent to former defgrad_av
enddo; enddo
! stress boundary condition check -------------
pstress_av = math_rotate_forward3x3(pstress_av_lab,bc(loadcase)%rotation)
print '(a,/,3(3(f12.7,x)/)$)', 'Piola-Kirchhoff stress / MPa: ',math_transpose3x3(pstress_av)/1.e6
if(size_reduced > 0_pInt) then ! calculate stress BC if applied
err_stress = maxval(abs(mask_stress * (pstress_av - bc(loadcase)%stress))) ! maximum deviaton (tensor norm not applicable)
err_stress_tol = maxval(abs(mask_defgrad * pstress_av)) * err_stress_tolrel ! don't use any tensor norm because the comparison should be coherent
print '(a)', ''
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 - bc(loadcase)%stress))) ! residual on given stress components
defgradAim = defgradAim + defgradAimCorr
print '(a,/,3(3(f12.7,x)/)$)', 'new deformation aim: ', math_transpose3x3(defgradAim)
print '(a,x,es10.4)' , 'with determinant: ', math_det3x3(defgradAim)
else
err_stress_tol = 0.0_pReal
endif
! ------------------------------
! homogeneous correction towards avg deformation gradient -------------
defgradAim_lab = math_rotate_backward3x3(defgradAim,bc(loadcase)%rotation) ! boundary conditions from load frame into lab (Fourier) frame
do m = 1_pInt,3_pInt; do n = 1_pInt,3_pInt
defgrad(1:res(1),1:res(2),1:res(3),m,n) = &
defgrad(1:res(1),1:res(2),1:res(3),m,n) + (defgradAim_lab(m,n) - defgrad_av_lab(m,n)) ! anticipated target minus current state
enddo; enddo
! ------------------------------
! bounds of det(F) -------------
defgradDetMax = -huge(1.0_pReal)
defgradDetMin = +huge(1.0_pReal)
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
defgradDet = math_det3x3(defgrad(i,j,k,1:3,1:3))
defgradDetMax = max(defgradDetMax,defgradDet)
defgradDetMin = min(defgradDetMin,defgradDet)
enddo; enddo; enddo
print '(a,x,es10.4)' , 'max determinant of deformation:', defgradDetMax
print '(a,x,es10.4)' , 'min determinant of deformation:', defgradDetMin
! ------------------------------
enddo ! end looping when convergency is achieved enddo ! end looping when convergency is achieved
c_prev = math_rotate_forward3x3x3x3(c_current*wgt,bc(loadcase)%rotation) ! calculate stiffness for next step
!ToDo: Incfluence for next loadcase
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
print '(a)', '' print '(a)', ''
print '(a)', '=============================================================' print '(a)', '============================================================='
@ -894,7 +903,7 @@ program DAMASK_spectral
write(538), materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! write result 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 ! end calculation/forwarding
enddo ! end looping over steps in current loadcase enddo ! end looping over steps in current loadcase
deallocate(c_reduced) deallocate(c_reduced)
deallocate(s_reduced) deallocate(s_reduced)
@ -902,7 +911,8 @@ program DAMASK_spectral
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
print '(a)', '' 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))

32
code/FEsolving.f90
View File

@ -62,28 +62,22 @@
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
call get_command(commandLine) ! may contain uppercase call get_command(commandLine) ! may contain uppercase
do i=1,len(commandLine) do i=1,len(commandLine)
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) ! make lowercase 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 if (index(commandLine,'-r ',.true.)>0) & ! look for -r
if (index(commandLine,'--restart ',.true.)>0) then ! if '--restart' is found, use that (contains '-r') start = index(commandLine,'-r ',.true.) + 3_pInt ! set to position after trailing space
start = index(commandLine,'--restart ',.true.) + 10_pInt if (index(commandLine,'--restart ',.true.)>0) & ! look for --restart
endif start = index(commandLine,'--restart ',.true.) + 10_pInt ! set to position after trailing space
length = index(commandLine(start:len(commandLine)),' ',.false.)
if(start /= 3_pInt) then ! found at least -r if(start /= 0_pInt) then ! found something
read(commandLine(start:start+length),'(I12)') restartReadStep length = verify(commandLine(start:len(commandLine)),'0123456789',.false.) ! where is first non number after argument?
read(commandLine(start:start+length),'(I12)') restartReadStep ! read argument
restartRead = max(0_pInt,restartReadStep) > 0_pInt
if(restartReadStep < 0_pInt) call IO_warning(warning_ID=34_pInt)
endif endif
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
@ -142,7 +136,7 @@
endif endif
enddo enddo
elseif (FEsolver == 'Spectral') then elseif (FEsolver == 'Spectral') then
restartReadSpectral = .true. !do nothing
else else
call IO_error(106) ! cannot open file for old job info call IO_error(106) ! cannot open file for old job info
endif endif

4
code/IO.f90
View File

@ -1166,8 +1166,6 @@ endfunction
msg = 'incomplete information in spectral mesh header' msg = 'incomplete information in spectral mesh header'
case (46) case (46)
msg = 'not a rotation defined for loadcase rotation' msg = 'not a rotation defined for loadcase rotation'
case (47)
msg = 'invalid restart increment given'
case (50) case (50)
msg = 'writing constitutive output description' msg = 'writing constitutive output description'
case (100) case (100)
@ -1416,6 +1414,8 @@ endfunction
select case (warning_ID) select case (warning_ID)
case (33_pInt) case (33_pInt)
msg = 'cannot guess along trajectory for first step of first loadcase' msg = 'cannot guess along trajectory for first step of first loadcase'
case (34)
msg = 'invalid restart increment given'
case (47_pInt) case (47_pInt)
msg = 'No valid parameter for FFTW given, using FFTW_PATIENT' msg = 'No valid parameter for FFTW given, using FFTW_PATIENT'
case (101_pInt) case (101_pInt)

6
code/constitutive_dislotwin.f90
View File

@ -1057,7 +1057,8 @@ if(constitutive_dislotwin_sbVelocity(myInstance) /= 0.0_pReal) then
!* Stress ratios !* Stress ratios
StressRatio_p = (abs(tau_sb(j))/constitutive_dislotwin_sbResistance(myInstance))**constitutive_dislotwin_p(myInstance) StressRatio_p = (abs(tau_sb(j))/constitutive_dislotwin_sbResistance(myInstance))**constitutive_dislotwin_p(myInstance)
StressRatio_pminus1 = (abs(tau_sb(j))/constitutive_dislotwin_sbResistance(myInstance))**(constitutive_dislotwin_p(myInstance)-1.0_pReal) StressRatio_pminus1 = (abs(tau_sb(j))/constitutive_dislotwin_sbResistance(myInstance))&
**(constitutive_dislotwin_p(myInstance)-1.0_pReal)
!* Boltzmann ratio !* Boltzmann ratio
BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(f,myInstance)/(kB*Temperature) BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(f,myInstance)/(kB*Temperature)
!* Initial shear rates !* Initial shear rates
@ -1463,7 +1464,8 @@ do o = 1,phase_Noutput(material_phase(g,ip,el))
tau = dot_product(Tstar_v,constitutive_dislotwin_sbSv(1:6,j,g,ip,el)) tau = dot_product(Tstar_v,constitutive_dislotwin_sbSv(1:6,j,g,ip,el))
!* Stress ratios !* Stress ratios
StressRatio_p = (abs(tau)/constitutive_dislotwin_sbResistance(myInstance))**constitutive_dislotwin_p(myInstance) StressRatio_p = (abs(tau)/constitutive_dislotwin_sbResistance(myInstance))**constitutive_dislotwin_p(myInstance)
StressRatio_pminus1 = (abs(tau)/constitutive_dislotwin_sbResistance(myInstance))**(constitutive_dislotwin_p(myInstance)-1.0_pReal) StressRatio_pminus1 = (abs(tau)/constitutive_dislotwin_sbResistance(myInstance))&
**(constitutive_dislotwin_p(myInstance)-1.0_pReal)
!* Boltzmann ratio !* Boltzmann ratio
BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(f,myInstance)/(kB*Temperature) BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(f,myInstance)/(kB*Temperature)
!* Initial shear rates !* Initial shear rates

3
code/crystallite.f90
View File

@ -3065,7 +3065,8 @@ logical error
! --- UPDATE SOME ADDITIONAL VARIABLES THAT ARE NEEDED FOR NONLOCAL MATERIAL --- ! --- UPDATE SOME ADDITIONAL VARIABLES THAT ARE NEEDED FOR NONLOCAL MATERIAL ---
! --- we use crystallite_orientation from above, so need a seperate loop ! --- we use crystallite_orientation from above, so need a seperate loop
!$OMP PARALLEL DO PRIVATE(myPhase,myInstance,myStructure,neighboring_e,neighboring_i,neighboringPhase,neighboringInstance,neighboringStructure) !$OMP PARALLEL DO PRIVATE(myPhase,myInstance,myStructure,neighboring_e,neighboring_i,neighboringPhase,&
!$OMP neighboringInstance,neighboringStructure)
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
myPhase = material_phase(1,i,e) ! get my phase myPhase = material_phase(1,i,e) ! get my phase

17
code/debug.f90
View File

@ -115,7 +115,7 @@ subroutine debug_init()
select case(tag) select case(tag)
case ('element','e','el') case ('element','e','el')
debug_e = IO_intValue(line,positions,2) debug_e = IO_intValue(line,positions,2)
case ('itegrationpoint','i','ip') case ('integrationpoint','i','ip')
debug_i = IO_intValue(line,positions,2) debug_i = IO_intValue(line,positions,2)
case ('grain','g','gr') case ('grain','g','gr')
debug_g = IO_intValue(line,positions,2) debug_g = IO_intValue(line,positions,2)
@ -123,12 +123,15 @@ subroutine debug_init()
debug_selectiveDebugger = IO_intValue(line,positions,2) > 0_pInt debug_selectiveDebugger = IO_intValue(line,positions,2) > 0_pInt
case ('verbosity') case ('verbosity')
debug_verbosity = IO_intValue(line,positions,2) debug_verbosity = IO_intValue(line,positions,2)
case ('(generaldebugspectral)') ! use bitwise logical and, continue with +8_pInt case ('(spectral)')
if(IO_intValue(line,positions,2)) spectral_debug_verbosity = spectral_debug_verbosity + 1_pInt select case(IO_lc(IO_stringValue(line,positions,2)))
case ('(divergencedebugspectral)') case('general')
if(IO_intValue(line,positions,2)) spectral_debug_verbosity = spectral_debug_verbosity + 2_pInt spectral_debug_verbosity = ior(spectral_debug_verbosity, 1)
case ('(restartdebugspectral)') case('divergence')
if(IO_intValue(line,positions,2)) spectral_debug_verbosity = spectral_debug_verbosity + 4_pInt spectral_debug_verbosity = ior(spectral_debug_verbosity, 2)
case('restart')
spectral_debug_verbosity = ior(spectral_debug_verbosity, 4)
endselect
endselect endselect
enddo enddo
100 close(fileunit) 100 close(fileunit)

56
code/makefile
View File

@ -25,10 +25,13 @@
# PORTABLE = TRUE (FALSE): decision, if executable is optimized for the machine on which it was built. # PORTABLE = TRUE (FALSE): decision, if executable is optimized for the machine on which it was built.
# OPTIMIZATION = DEFENSIVE (OFF,AGGRESSIVE,ULTRA): Optimization mode: O2, O0, O3 + further options for most files, 03 + further options for all files # OPTIMIZATION = DEFENSIVE (OFF,AGGRESSIVE,ULTRA): Optimization mode: O2, O0, O3 + further options for most files, 03 + further options for all files
# OPENMP = TRUE (FALSE): OpenMP multiprocessor support # OPENMP = TRUE (FALSE): OpenMP multiprocessor support
# ACMLPATH =/opt/acml4.4.0/ifort64/lib (...) Path to ACML Library, choose according to your system # FFTWPATH =TAKE_FFTW_PATH, will be adjusted by setup_code.py
# ACMLPATH =/opt/acml4.4.0/ifort64_mp/lib (...) Path to ACML Library with multicore support, choose according to your system # ACMLROOT =TAKE_ACLM_ROOT, will be adjusted by setup_code.py
# FFTWOPTIONS =./../lib/libfftw3.a include/libfftw3_threads.a -lpthread (...) Path to FFTW library files with Linux threads (multicore) support # ACMLPATH =$(ACMLROOT)/"compilerdir"/lib (...) Path to ACML Library, choose according to your system
# FFTWOPTIONS =./../lib/libfftw3.a (...) Path to FFTW library files without Linux threads (multicore) support # ACMLPATH =$(ACMLROOT/"compilerdir"_mp/lib (...) Path to ACML Library with multicore support, choose according to your system
# "compilerdir" is "intel64" for ifort
# FFTWOPTIONS =$(FFTWPATH)/libfftw3.a $(FFTWPATH)/libfftw3_threads.a -lpthread (...) Path to FFTW library files with Linux threads (multicore) support
# FFTWOPTIONS =$(FFTWPATH)/libfftw3.a (...) Path to FFTW library files without Linux threads (multicore) support
# FFTWOPTIONS is different for single and double precision. Choose the options to use OpenMP instead of pthreads support or change the directory # FFTWOPTIONS is different for single and double precision. Choose the options to use OpenMP instead of pthreads support or change the directory
# PREFIX: specify an arbitrary prefix # PREFIX: specify an arbitrary prefix
# SUFFIX: specify an arbitrary suffix # SUFFIX: specify an arbitrary suffix
@ -38,17 +41,29 @@
# information on http://software.intel.com/en-us/articles/determining-root-cause-of-sigsegv-or-sigbus-errors/ # information on http://software.intel.com/en-us/articles/determining-root-cause-of-sigsegv-or-sigbus-errors/
# check if an array index is too small (<1) or too large! # check if an array index is too small (<1) or too large!
DEBUG1 =-check bounds -g DEBUG1 =-check bounds -g
#will cause a lot of warnings because we create a bunch of temporary arrays #will cause a lot of warnings because we create a bunch of temporary arrays
DEBUG2 =-check arg_temp_created DEBUG2 =-check arg_temp_created
#check from time to time #check from time to time
DEBUG3 =-fp-stack-check -g -traceback -gen-interfaces -warn interfaces DEBUG3 =-fp-stack-check -g -traceback -gen-interfaces -warn interfaces
#should not be done for OpenMP, but set "ulimit -s unlimited" on shell. Problably it helps also to unlimit other limits #should not be done for OpenMP, but set "ulimit -s unlimited" on shell. Problably it helps also to unlimit other limits
DEBUG4 =-heap-arrays DEBUG4 =-heap-arrays
#checks for standard #checks for standard
DEBUG5 =stand std03/std95 DEBUG5 =stand std03/std95
#SUFFIX =$(DEBUG1) $(DEBUG2) $(DEBUG3)
#SUFFIX =$(DEBUG1) $(DEBUG2) $(DEBUG3)
######################################################################################## ########################################################################################
#default values below will be set by setup_code.py
#FFTWPATH =../lib
#ACMLROOT =/opt/acml4.4.0
FFTWPATH =/nethome/m.diehl/DAMASK/lib
ACMLROOT =/opt/acml4.4.0
ifndef F90 ifndef F90
F90 =ifort F90 =ifort
endif endif
@ -75,34 +90,35 @@ ifndef OPENMP
OPENMP=ON OPENMP=ON
endif endif
# setting defaults in case of multicore support
ifeq ($(OPENMP),ON) ifeq ($(OPENMP),ON)
OPENMP_FLAG_ifort =-openmp -openmp-report0 -parallel OPENMP_FLAG_ifort =-openmp -openmp-report0 -parallel
OPENMP_FLAG_gfortran =-fopenmp OPENMP_FLAG_gfortran =-fopenmp
ifndef ACMLPATH ifndef ACMLPATH
ACMLPATH =/opt/acml4.4.0/ifort64_mp/lib ACMLPATH =$(ACMLROOT)/$(F90)64_mp/lib
endif endif
ifndef FFTWOPTIONS ifndef FFTWOPTIONS
ifeq ($(PRECISION),single) ifeq ($(PRECISION),single)
FFTWOPTIONS =./../lib/libfftw3f_threads.a ./../lib/libfftw3f.a -lpthread FFTWOPTIONS =$(FFTWPATH)/libfftw3f_threads.a $(FFTWPATH)/libfftw3f.a -lpthread
else else
FFTWOPTIONS =./../lib/libfftw3_threads.a ./../lib/libfftw3.a -lpthread FFTWOPTIONS =$(FFTWPATH)/libfftw3_threads.a $(FFTWPATH)/libfftw3.a -lpthread
endif endif
endif endif
BLAS_ifort =-L $(ACMLPATH) -lacml_mp BLAS=-L $(ACMLPATH) -lacml_mp
BLAS_gfortran =
#setting defaults in case of single core compilation
else else
ifndef ACMLPATH ifndef ACMLPATH
ACMLPATH =/opt/acml4.4.0/ifort64/lib ACMLPATH=$(ACMLROOT)/$(F90)64/lib
endif endif
ifndef FFTWOPTIONS ifndef FFTWOPTIONS
ifeq ($(PRECISION),single) ifeq ($(PRECISION),single)
FFTWOPTIONS =./../lib/libfftw3f.a FFTWOPTIONS =$(FFTWPATH)/libfftw3f.a
else else
FFTWOPTIONS =./../lib/libfftw3.a FFTWOPTIONS =$(FFTWPATH)/libfftw3.a
endif endif
endif endif
BLAS_ifort =-L $(ACMLPATH) -lacml BLAS=-L $(ACMLPATH) -lacml
BLAS_gfortran =
endif endif
@ -111,10 +127,10 @@ OPTIMIZATION_OFF_gfortran =-O0
OPTIMIZATION_DEFENSIVE_ifort =-O2 OPTIMIZATION_DEFENSIVE_ifort =-O2
OPTIMIZATION_DEFENSIVE_gfortran =-O2 OPTIMIZATION_DEFENSIVE_gfortran =-O2
OPTIMIZATION_AGGRESSIVE_ifort =-O3 $(PORTABLE_SWITCH) -ip -static -fp-model fast=2 -no-prec-div OPTIMIZATION_AGGRESSIVE_ifort =-O3 $(PORTABLE_SWITCH) -ip -static -fp-model fast=2 -no-prec-div
OPTIMIZATION_AGGRESSIVE_gfortran =-O3 $(PORTABLE_SWITCH) -march=opteron -ffast-math -funroll-loops -ftree-vectorize -ftree-loop-linear OPTIMIZATION_AGGRESSIVE_gfortran =-O3 $(PORTABLE_SWITCH) -march=opteron -ffast-math -funroll-loops -ftree-vectorize -ftree-loop-linear
COMPILE_OPTIONS_ifort =-fpp -diag-disable 8291,8290 COMPILE_OPTIONS_ifort =-fpp -diag-disable 8291,8290
COMPILE_OPTIONS_gfortran =-xf95-cpp-input -ffree-line-length-none -fno-range-check COMPILE_OPTIONS_gfortran =-xf95-cpp-input -fno-range-check
COMPILE =${OPENMP_FLAG_${F90}} ${COMPILE_OPTIONS_${F90}} ${OPTIMIZATION_${OPTIMIZATION}_${F90}} -c COMPILE =${OPENMP_FLAG_${F90}} ${COMPILE_OPTIONS_${F90}} ${OPTIMIZATION_${OPTIMIZATION}_${F90}} -c
COMPILE_MAXOPTI =${OPENMP_FLAG_${F90}} ${COMPILE_OPTIONS_${F90}} ${OPTIMIZATION_${MAXOPTI}_${F90}} -c COMPILE_MAXOPTI =${OPENMP_FLAG_${F90}} ${COMPILE_OPTIONS_${F90}} ${OPTIMIZATION_${MAXOPTI}_${F90}} -c
@ -125,14 +141,14 @@ endif
ifeq ($(PRECISION),single) ifeq ($(PRECISION),single)
DAMASK_spectral_single.exe: DAMASK_spectral_single.o CPFEM.a DAMASK_spectral_single.exe: DAMASK_spectral_single.o CPFEM.a
$(PREFIX) $(COMPILERNAME) ${OPENMP_FLAG_${F90}} -o DAMASK_spectral_single.exe DAMASK_spectral_single.o CPFEM.a $(FFTWOPTIONS)\ $(PREFIX) $(COMPILERNAME) $(OPENMP_FLAG_$(F90)) -o DAMASK_spectral_single.exe DAMASK_spectral_single.o CPFEM.a $(FFTWOPTIONS)\
constitutive.a advanced.a basics.a ${BLAS_${F90}} constitutive.a advanced.a basics.a $(BLAS)
DAMASK_spectral_single.o: DAMASK_spectral_single.f90 CPFEM.o DAMASK_spectral_single.o: DAMASK_spectral_single.f90 CPFEM.o
$(PREFIX) $(COMPILERNAME) $(COMPILE_MAXOPTI) DAMASK_spectral_single.f90 $(SUFFIX) $(PREFIX) $(COMPILERNAME) $(COMPILE_MAXOPTI) DAMASK_spectral_single.f90 $(SUFFIX)
else else
DAMASK_spectral.exe: DAMASK_spectral.o CPFEM.a DAMASK_spectral.exe: DAMASK_spectral.o CPFEM.a
$(PREFIX) $(COMPILERNAME) ${OPENMP_FLAG_${F90}} -o DAMASK_spectral.exe DAMASK_spectral.o CPFEM.a $(FFTWOPTIONS)\ $(PREFIX) $(COMPILERNAME) ${OPENMP_FLAG_${F90}} -o DAMASK_spectral.exe DAMASK_spectral.o CPFEM.a $(FFTWOPTIONS)\
constitutive.a advanced.a basics.a ${BLAS_${F90}} constitutive.a advanced.a basics.a $(BLAS)
DAMASK_spectral.o: DAMASK_spectral.f90 CPFEM.o DAMASK_spectral.o: DAMASK_spectral.f90 CPFEM.o
$(PREFIX) $(COMPILERNAME) $(COMPILE_MAXOPTI) DAMASK_spectral.f90 $(SUFFIX) $(PREFIX) $(COMPILERNAME) $(COMPILE_MAXOPTI) DAMASK_spectral.f90 $(SUFFIX)
endif endif

31
code/math.f90
View File

@ -3048,7 +3048,8 @@ subroutine mesh_regular_grid(res,geomdim,defgrad_av,centroids,nodes)
me = (/i,j,k/) ! me on skin me = (/i,j,k/) ! me on skin
shift = sign(abs(res+diag-2_pInt*me)/(res+diag),res+diag-2_pInt*me) shift = sign(abs(res+diag-2_pInt*me)/(res+diag),res+diag-2_pInt*me)
lookup = me-diag+shift*res lookup = me-diag+shift*res
wrappedCentroids(i+1_pInt,j+1_pInt,k+1_pInt,1:3) = centroids(lookup(1)+1_pInt,lookup(2)+1_pInt,lookup(3)+1_pInt,1:3) - & wrappedCentroids(i+1_pInt,j+1_pInt,k+1_pInt,1:3) = &
centroids(lookup(1)+1_pInt,lookup(2)+1_pInt,lookup(3)+1_pInt,1:3) - &
matmul(defgrad_av, shift*geomdim) matmul(defgrad_av, shift*geomdim)
endif endif
enddo; enddo; enddo enddo; enddo; enddo
@ -3056,7 +3057,8 @@ subroutine mesh_regular_grid(res,geomdim,defgrad_av,centroids,nodes)
do j = 0_pInt,res(2) do j = 0_pInt,res(2)
do i = 0_pInt,res(1) do i = 0_pInt,res(1)
do n = 1_pInt,8_pInt do n = 1_pInt,8_pInt
nodes(i+1_pInt,j+1_pInt,k+1_pInt,1:3) = nodes(i+1_pInt,j+1_pInt,k+1_pInt,3) + wrappedCentroids(i+1_pInt+neighbor(1_pInt,n), & nodes(i+1_pInt,j+1_pInt,k+1_pInt,1:3) = &
nodes(i+1_pInt,j+1_pInt,k+1_pInt,3) + wrappedCentroids(i+1_pInt+neighbor(1_pInt,n), &
j+1_pInt+neighbor(2,n), & j+1_pInt+neighbor(2,n), &
k+1_pInt+neighbor(3,n),1:3) k+1_pInt+neighbor(3,n),1:3)
enddo; enddo; enddo; enddo enddo; enddo; enddo; enddo
@ -3197,6 +3199,7 @@ subroutine deformed_fft(res,geomdim,defgrad_av,scaling,defgrad,coords)
! other variables ! other variables
integer(pInt) :: i, j, k integer(pInt) :: i, j, k
integer(pInt), dimension(3) :: k_s integer(pInt), dimension(3) :: k_s
complex(pReal), parameter :: integration_factor = cmplx(0.0_pReal,pi*2.0_pReal)
real(pReal), dimension(3) :: step, offset_coords real(pReal), dimension(3) :: step, offset_coords
integer*8, dimension(2) :: plan_fft integer*8, dimension(2) :: plan_fft
@ -3228,11 +3231,11 @@ subroutine deformed_fft(res,geomdim,defgrad_av,scaling,defgrad,coords)
do i = 1_pInt, res(1)/2_pInt+1_pInt do i = 1_pInt, res(1)/2_pInt+1_pInt
k_s(1) = i-1_pInt k_s(1) = i-1_pInt
if(i/=1_pInt) coords_fft(i,j,k,1:3) = coords_fft(i,j,k,1:3)& if(i/=1_pInt) coords_fft(i,j,k,1:3) = coords_fft(i,j,k,1:3)&
+ defgrad_fft(i,j,k,1:3,1)*geomdim(1)/(real(k_s(1),pReal)*cmplx(0.0_pReal,1.0_pReal)*pi*2.0_pReal) + defgrad_fft(i,j,k,1:3,1)*geomdim(1)/(real(k_s(1),pReal)*integration_factor)
if(j/=1_pInt) coords_fft(i,j,k,1:3) = coords_fft(i,j,k,1:3)& if(j/=1_pInt) coords_fft(i,j,k,1:3) = coords_fft(i,j,k,1:3)&
+ defgrad_fft(i,j,k,1:3,2)*geomdim(2)/(real(k_s(2),pReal)*cmplx(0.0_pReal,1.0_pReal)*pi*2.0_pReal) + defgrad_fft(i,j,k,1:3,2)*geomdim(2)/(real(k_s(2),pReal)*integration_factor)
if(k/=1_pInt) coords_fft(i,j,k,1:3) = coords_fft(i,j,k,1:3)& if(k/=1_pInt) coords_fft(i,j,k,1:3) = coords_fft(i,j,k,1:3)&
+ defgrad_fft(i,j,k,1:3,3)*geomdim(3)/(real(k_s(3),pReal)*cmplx(0.0_pReal,1.0_pReal)*pi*2.0_pReal) + defgrad_fft(i,j,k,1:3,3)*geomdim(3)/(real(k_s(3),pReal)*integration_factor)
enddo; enddo; enddo enddo; enddo; enddo
call dfftw_execute_dft_c2r(plan_fft(2), coords_fft, coords) call dfftw_execute_dft_c2r(plan_fft(2), coords_fft, coords)
@ -3432,12 +3435,18 @@ subroutine divergence_fft(res,geomdim,vec_tens,field,divergence_field)
order = order + 1_pInt order = order + 1_pInt
do k = 0_pInt, res(3)-1_pInt; do j = 0_pInt, res(2)-1_pInt; do i = 0_pInt, res(1)-1_pInt do k = 0_pInt, res(3)-1_pInt; do j = 0_pInt, res(2)-1_pInt; do i = 0_pInt, res(1)-1_pInt
do m = 1_pInt, order do m = 1_pInt, order
coordinates(1,1:3) = mesh_location(mesh_index((/i+m,j,k/),(/res(1),res(2),res(3)/)),(/res(1),res(2),res(3)/)) + (/1_pInt,1_pInt,1_pInt/) coordinates(1,1:3) = mesh_location(mesh_index((/i+m,j,k/),(/res(1),res(2),res(3)/)),(/res(1),res(2),res(3)/))&
coordinates(2,1:3) = mesh_location(mesh_index((/i-m,j,k/),(/res(1),res(2),res(3)/)),(/res(1),res(2),res(3)/)) + (/1_pInt,1_pInt,1_pInt/) + (/1_pInt,1_pInt,1_pInt/)
coordinates(3,1:3) = mesh_location(mesh_index((/i,j+m,k/),(/res(1),res(2),res(3)/)),(/res(1),res(2),res(3)/)) + (/1_pInt,1_pInt,1_pInt/) coordinates(2,1:3) = mesh_location(mesh_index((/i-m,j,k/),(/res(1),res(2),res(3)/)),(/res(1),res(2),res(3)/))&
coordinates(4,1:3) = mesh_location(mesh_index((/i,j-m,k/),(/res(1),res(2),res(3)/)),(/res(1),res(2),res(3)/)) + (/1_pInt,1_pInt,1_pInt/) + (/1_pInt,1_pInt,1_pInt/)
coordinates(5,1:3) = mesh_location(mesh_index((/i,j,k+m/),(/res(1),res(2),res(3)/)),(/res(1),res(2),res(3)/)) + (/1_pInt,1_pInt,1_pInt/) coordinates(3,1:3) = mesh_location(mesh_index((/i,j+m,k/),(/res(1),res(2),res(3)/)),(/res(1),res(2),res(3)/))&
coordinates(6,1:3) = mesh_location(mesh_index((/i,j,k-m/),(/res(1),res(2),res(3)/)),(/res(1),res(2),res(3)/)) + (/1_pInt,1_pInt,1_pInt/) + (/1_pInt,1_pInt,1_pInt/)
coordinates(4,1:3) = mesh_location(mesh_index((/i,j-m,k/),(/res(1),res(2),res(3)/)),(/res(1),res(2),res(3)/))&
+ (/1_pInt,1_pInt,1_pInt/)
coordinates(5,1:3) = mesh_location(mesh_index((/i,j,k+m/),(/res(1),res(2),res(3)/)),(/res(1),res(2),res(3)/))&
+ (/1_pInt,1_pInt,1_pInt/)
coordinates(6,1:3) = mesh_location(mesh_index((/i,j,k-m/),(/res(1),res(2),res(3)/)),(/res(1),res(2),res(3)/))&
+ (/1_pInt,1_pInt,1_pInt/)
do l = 1_pInt, vec_tens do l = 1_pInt, vec_tens
divergence_field(i+1_pInt,j+1_pInt,k+1_pInt,l) = divergence_field(i+1_pInt,j+1_pInt,k+1_pInt,l) + FDcoefficient(m,order) * & divergence_field(i+1_pInt,j+1_pInt,k+1_pInt,l) = divergence_field(i+1_pInt,j+1_pInt,k+1_pInt,l) + FDcoefficient(m,order) * &
((field(coordinates(1,1),coordinates(1,2),coordinates(1,3),l,1)- & ((field(coordinates(1,1),coordinates(1,2),coordinates(1,3),l,1)- &

26
code/setup/setup_code.py
View File

@ -1,7 +1,7 @@
#!/usr/bin/env python #!/usr/bin/env python
# $Id$ # $Id$
# Writes version specific files for different MARC releases # Writes version specific files for different MARC releases
import os,sys import os,sys,string,damask_tools
architectures = { architectures = {
'marc': { 'marc': {
@ -29,8 +29,28 @@ for arch in architectures:
childFile.write(line.replace(me['versions'][0],version)) childFile.write(line.replace(me['versions'][0],version))
childFile.close() childFile.close()
if raw_input("Do you want to compile the spectral code now? (y/n)") is 'y': # changing dirs in make file
os.system('make --directory %s'%wd) damask_variables = damask_tools.DAMASK_TOOLS()
makefile = open(os.path.join(damask_variables.rootDir(),'code/makefile'))
content = makefile.readlines()
makefile.close()
makefile = open(os.path.join(damask_variables.rootDir(),'code/makefile'),'w')
for line in content:
if line.startswith('FFTWPATH'):
line='FFTWPATH =%s\n'%(damask_variables.pathInfo['fftw'])
print line
if line.startswith('ACMLROOT'):
line='ACMLROOT =%s\n'%(damask_variables.pathInfo['acml'])
print line
makefile.writelines(line)
makefile.close()
# compiling spectral code
if raw_input("Do you want to compile the spectral code now? (y/n) ") is 'y' or 'Y':
compiler_switches = raw_input("Please give compiling switches (Enter to use default) ")
os.system('make --directory %s clean'%(wd))
print compiler_switches
os.system('make --directory %s %s'%(wd,compiler_switches))
if '--clean' in [s.lower() for s in sys.argv]: if '--clean' in [s.lower() for s in sys.argv]:
os.system('make --directory %s clean'%wd) os.system('make --directory %s clean'%wd)

4
processing/post/3Dvisualize.py
View File

@ -310,7 +310,7 @@ parser.add_option('--box', dest='output_box', action='store_true', \
help='produce VTK box file') help='produce VTK box file')
parser.add_option('--nobox', dest='output_box', action='store_false', \ parser.add_option('--nobox', dest='output_box', action='store_false', \
help='omit VTK box file') help='omit VTK box file')
parser.add_option('--scaling', dest='scaling', type='int', \ parser.add_option('--scaling', dest='scaling', type='float', \
help='scaling of fluctuation [%default]') help='scaling of fluctuation [%default]')
parser.add_option('-u', '--unitlength', dest='unitlength', type='float', \ parser.add_option('-u', '--unitlength', dest='unitlength', type='float', \
help='set unit length for 2D model [%default]') help='set unit length for 2D model [%default]')
@ -321,7 +321,7 @@ parser.set_defaults(tensor = [])
parser.set_defaults(output_mesh = True) parser.set_defaults(output_mesh = True)
parser.set_defaults(output_points = False) parser.set_defaults(output_points = False)
parser.set_defaults(output_box = False) parser.set_defaults(output_box = False)
parser.set_defaults(scaling = 1) parser.set_defaults(scaling = 1.0)
parser.set_defaults(undeformed = False) parser.set_defaults(undeformed = False)
parser.set_defaults(unitlength = 0.0) parser.set_defaults(unitlength = 0.0)
parser.set_defaults(cell = True) parser.set_defaults(cell = True)

26
processing/post/make_DAMASK2Python
View File

@ -1,26 +0,0 @@
#!/bin/bash
# $Id: make_postprocessingMath 1106 2011-11-17 21:36:56Z MPIE\m.diehl $
# This script is used to compile math.f90 and make the functions defined in DAMASK_math.pyf
# avialable for python in the module DAMASK_math.so
# It uses the fortran wrapper f2py that is included in the numpy package to construct the
# module postprocessingMath.so out of the fortran code postprocessingMath.f90
# for the generation of the pyf file:
#f2py -m DAMASK -h DAMASK.pyf --overwrite-signature ../../code/math.f90 \
if [[ $# -eq 0 ]]; then
wd='.'
else
wd=$1
fi
# use env. variables here!
cd $wd
rm ../../lib/DAMASK.so
f2py -c \
DAMASK.pyf \
../../code/DAMASK2Python_helper.f90 \
../../code/math.f90 \
../../lib/libfftw3.a \
/opt/acml4.4.0/ifort64/lib/libacml.a
mv DAMASK.so ../../lib/.

4
processing/pre/voronoi_randomSeeding.f90
View File

@ -81,8 +81,8 @@ program voronoi
open(21, file = trim(filename)//('.seeds')) open(21, file = trim(filename)//('.seeds'))
write(21, '(i1,a1,a6)') 4,achar(9),'header' write(21, '(i1,a1,a6)') 4,achar(9),'header'
write(21, '(A, I8, A, I8, A, I8)') 'resolution a ', a, ' b ', b, ' c ', c write(21, '(A, I8, A, I8, A, I8)') 'resolution a ', a, ' b ', b, ' c ', c
write(21,*) 'grains', N_Seeds write(21, '(A, I8)') 'grains', N_Seeds
write(21,*) 'random seed ',rndInit(1) write(21, '(A, I8)') 'random seed ',rndInit(1)
write(21,'(6(a,a1))') 'x',achar(9),'y',achar(9),'z',achar(9),'phi1',achar(9),'Phi',achar(9),'phi2',achar(9) write(21,'(6(a,a1))') 'x',achar(9),'y',achar(9),'z',achar(9),'phi1',achar(9),'Phi',achar(9),'phi2',achar(9)
do i = 1, n_Seeds do i = 1, n_Seeds
write(21, '(6(F10.6,a1))'),seeds(i,1), achar(9), seeds(i,2), achar(9), seeds(i,3), achar(9), & write(21, '(6(F10.6,a1))'),seeds(i,1), achar(9), seeds(i,2), achar(9), seeds(i,3), achar(9), &

10
processing/pre/voronoi_tessellation.f90
View File

@ -169,7 +169,7 @@ program voronoi
use IO use IO
implicit none implicit none
logical gotN_Seeds, gotResolution logical :: gotN_Seeds=.false., gotResolution = .false.
logical, dimension(:), allocatable :: grainCheck logical, dimension(:), allocatable :: grainCheck
character(len=1024) input_name, output_name, format1, format2, N_Digits, line, key character(len=1024) input_name, output_name, format1, format2, N_Digits, line, key
integer(pInt) N_Seeds, theGrain, i, j, k, l, m integer(pInt) N_Seeds, theGrain, i, j, k, l, m
@ -273,14 +273,14 @@ program voronoi
if(geomdim(i) .gt. 0.0) validDim(i) = .true. if(geomdim(i) .gt. 0.0) validDim(i) = .true.
enddo enddo
if(all(validDim == .false.)) then if(all(validDim .eqv. .false.)) then
geomdim(maxval(maxloc(resolution))) = 1.0 geomdim(maxval(maxloc(resolution))) = 1.0
validDim(maxval(maxloc(resolution))) = .true. validDim(maxval(maxloc(resolution))) = .true.
print*, 'no valid dimension specified, using automated setting' print*, 'no valid dimension specified, using automated setting'
endif endif
do i=1,3 do i=1,3
if (validDim(i) == .false.) then if (validDim(i) .eqv. .false.) then
print*, 'rescaling ivalid dimension' , i print*, 'rescaling ivalid dimension' , i
geomdim(i) = maxval(geomdim/real(resolution),validDim)*real(resolution(i)) geomdim(i) = maxval(geomdim/real(resolution),validDim)*real(resolution(i))
endif endif
@ -317,11 +317,11 @@ program voronoi
open(21, file = ((trim(output_name))//'.spectral')) open(21, file = ((trim(output_name))//'.spectral'))
write(20, '(A)'), '3 header' write(20, '(A)'), '3 header'
write(20, '(A, I8, A, I8, A, I8)'), 'resolution a ', resolution(1), ' b ', resolution(2), ' c ', resolution(3) write(20, '(A, I8, A, I8, A, I8)'), 'resolution a ', resolution(1), ' b ', resolution(2), ' c ', resolution(3)
write(20, '(A, g15.10, A, g15.10, A, g15.10)'), 'dimension x ', geomdim(1), ' y ', geomdim(2), ' z ', geomdim(3) write(20, '(A, g17.10, A, g17.10, A, g17.10)'), 'dimension x ', geomdim(1), ' y ', geomdim(2), ' z ', geomdim(3)
write(20, '(A)'), 'homogenization 1' write(20, '(A)'), 'homogenization 1'
format1 = '(I'//trim(N_Digits)//'.'//trim(N_Digits)//')' ! geom format format1 = '(I'//trim(N_Digits)//'.'//trim(N_Digits)//')' ! geom format
format2 = '(3(tr2, f6.2), 3(tr2,g10.5), I10, a)' ! spectral (Lebensohn) format format2 = '(3(tr2, f6.2), 3(tr2,g10.5), I10, a)' ! spectral (Lebensohn) format
! perform voronoi tessellation and write result to files ! perform voronoi tessellation and write result to files

137
processing/setup/setup_processing.py
View File

@ -2,11 +2,62 @@
# Makes postprocessing routines acessible from everywhere. # Makes postprocessing routines acessible from everywhere.
import os,sys,glob import os,sys,glob,string
from optparse import OptionParser, Option
import damask_tools; reload(damask_tools) import damask_tools
damask_tools.DAMASK_TOOLS().check_env()
# -----------------------------
class extendableOption(Option):
# -----------------------------
# used for definition of new option parser action 'extend', which enables to take multiple option arguments
# taken from online tutorial http://docs.python.org/library/optparse.html
ACTIONS = Option.ACTIONS + ("extend",)
STORE_ACTIONS = Option.STORE_ACTIONS + ("extend",)
TYPED_ACTIONS = Option.TYPED_ACTIONS + ("extend",)
ALWAYS_TYPED_ACTIONS = Option.ALWAYS_TYPED_ACTIONS + ("extend",)
def take_action(self, action, dest, opt, value, values, parser):
if action == "extend":
lvalue = value.split(",")
values.ensure_value(dest, []).extend(lvalue)
else:
Option.take_action(self, action, dest, opt, value, values, parser)
########################################################
# MAIN
########################################################
parser = OptionParser(option_class=extendableOption, usage='%prog options', description = """
Sets up the pre and post processing tools of DAMASK
""" + string.replace('$Id: addDivergence.py 1129 2011-12-01 12:01:13Z MPIE\m.diehl $','\n','\\n')
)
parser.add_option('--F90', dest='compiler', type='string', \
help='name of F90 compiler [%default]')
parser.set_defaults(compiler = 'ifort')
(options,filenames) = parser.parse_args()
#translating name of compiler for use with f2py and setting subdirname of acml
if options.compiler == 'gfortran':
f2py_compiler='gnu95 --f90flags="-fno-range-check"'
acml_subdir='ifort64/lib'
else:
f2py_compiler='intelem'
acml_subdir='ifort64/lib'
#getting pathinfo
damask_variables = damask_tools.DAMASK_TOOLS()
baseDir = os.path.join(damask_variables.rootDir(),'processing/')
codeDir = os.path.join(damask_variables.rootDir(),'code/')
#define ToDo list
bin_link = { \ bin_link = { \
'pre' : [ 'pre' : [
'marc_addUserOutput.py', 'marc_addUserOutput.py',
@ -16,8 +67,8 @@ bin_link = { \
'spectral_geomCheck.py', 'spectral_geomCheck.py',
'spectral_minimalSurface.py', 'spectral_minimalSurface.py',
'spectral_vicinityOffset.py', 'spectral_vicinityOffset.py',
'voronoi_randomSeeding.py', 'voronoi_randomSeeding.exe',
'voronoi_tessellation.py', 'voronoi_tessellation.exe',
], ],
'post' : [ 'post' : [
'3Dvisualize.py', '3Dvisualize.py',
@ -46,28 +97,58 @@ compile = { \
'post' : [ 'post' : [
], ],
} }
execute = { \ execute = { \
'pre' : [ 'postMath' : [
], 'rm %slib/DAMASK.so' %(os.path.join(damask_variables.rootDir(),'lib/')),
'post' : [ # The following command is used to compile math.f90 and make the functions defined in DAMASK_math.pyf
'make_DAMASK2Python', # avialable for python in the module DAMASK_math.so
], # It uses the fortran wrapper f2py that is included in the numpy package to construct the
# module postprocessingMath.so out of the fortran code postprocessingMath.f90
# for the generation of the pyf file:
#f2py -m DAMASK -h DAMASK.pyf --overwrite-signature ../../code/math.f90 \
'f2py %sDAMASK.pyf '%(codeDir) +\
'-c --fcompiler=%s '%(f2py_compiler) +\
'%sDAMASK2Python_helper.f90 ' %(codeDir)+\
'%smath.f90 ' %(codeDir)+\
'%s/libfftw3.a ' %(damask_variables.pathInfo['fftw'])+\
'%s/%s/libacml.a' %(damask_variables.pathInfo['acml'],acml_subdir),
'mv %sDAMASK.so %s/.' %(codeDir,os.path.join(damask_variables.rootDir(),'lib/')),
]
} }
damask_root = os.getenv('DAMASK_ROOT')
if damask_root == '' or damask_root == None: damask_root = os.path.join(os.path.dirname(sys.argv[0]),'../../')
damask_bin = os.getenv('DAMASK_BIN')
if damask_bin == '' or damask_bin == None: damask_bin = os.path.join(damask_root,'bin/')
basedir = os.path.join(damask_root,'processing/')
for dir in compile:
for file in compile[dir]:
src = os.path.abspath(os.path.join(baseDir,dir,file))
if os.path.isfile(src):
try:
os.system('rm %s.exe'%(os.path.splitext(src)[0]))
print 'removing %s.exe '%(os.path.splitext(src)[0])
except:
pass
print 'compiling ',src,'using',options.compiler
os.system('%s -O2 -o%s.exe %s'%(options.compiler,os.path.splitext(src)[0],src))
os.chdir(codeDir) # needed for compilation with gfortran and f2py
for tasks in execute:
for cmd in execute[tasks]:
os.system(cmd)
os.chdir(os.path.join(damask_variables.rootDir(),'processing/setup/'))
modules = glob.glob('*.mod')
for module in modules:
print 'removing', module
os.remove(module)
for dir in bin_link: for dir in bin_link:
for file in bin_link[dir]: for file in bin_link[dir]:
src = os.path.abspath(os.path.join(basedir,dir,file)) src = os.path.abspath(os.path.join(baseDir,dir,file))
if (file == ''): if (file == ''):
sym_link = os.path.abspath(os.path.join(damask_bin,dir)) sym_link = os.path.abspath(os.path.join(damask_variables.binDir(),dir))
else: else:
sym_link = os.path.abspath(os.path.join(damask_bin,os.path.splitext(file)[0])) sym_link = os.path.abspath(os.path.join(damask_variables.binDir(),os.path.splitext(file)[0]))
print sym_link,'-->',src print sym_link,'-->',src
if os.path.lexists(sym_link): if os.path.lexists(sym_link):
os.remove(sym_link) os.remove(sym_link)
@ -78,23 +159,3 @@ for dir in bin_link:
#if os.path.lexists(old_link): #if os.path.lexists(old_link):
# os.remove(old_link) # os.remove(old_link)
for dir in compile:
for file in compile[dir]:
src = os.path.abspath(os.path.join(basedir,dir,file))
if os.path.isfile(src):
print file
os.system('rm %s.exe'%(os.path.splitext(src)[0]))
os.system('ifort -O3 -parallel -o%s %s'%(os.path.splitext(src)[0],src))
modules = glob.glob('*.mod')
for module in modules:
print module
os.remove(module)
for dir in execute:
for file in execute[dir]:
cmd = os.path.abspath(os.path.join(basedir,dir,file))
if os.access(cmd,os.X_OK):
print cmd
os.system('%s %s'%(cmd,os.path.join(basedir,dir)))