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

polishing: removed variable names like 'unit' and 'data' that are keywords of fortran and ensured that integer and real precision matches independent of machine standard. 

removed cut_off parameter for damask_spectral
removed outpot of derived divergence measures and added RMS output in brackets
added comments and options to the makefile
This commit is contained in:
Martin Diehl 2012-02-15 18:58:38 +00:00
parent d642730776
commit 6c0f9d163b
16 changed files with 1120 additions and 1129 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

69
code/DAMASK_spectral.f90
View File

@ -53,7 +53,6 @@ program DAMASK_spectral
use numerics, only: err_div_tol, err_stress_tolrel, rotation_tol, itmax, & use numerics, only: err_div_tol, err_stress_tolrel, rotation_tol, itmax, &
memory_efficient, update_gamma, & memory_efficient, update_gamma, &
simplified_algorithm, divergence_correction, & simplified_algorithm, divergence_correction, &
cut_off_value, &
DAMASK_NumThreadsInt, & DAMASK_NumThreadsInt, &
fftw_planner_flag, fftw_timelimit fftw_planner_flag, fftw_timelimit
use homogenization, only: materialpoint_sizeResults, materialpoint_results use homogenization, only: materialpoint_sizeResults, materialpoint_results
@ -121,17 +120,8 @@ program DAMASK_spectral
s0_reference s0_reference
real(pReal), dimension(6) :: cstress ! cauchy stress real(pReal), dimension(6) :: cstress ! cauchy stress
real(pReal), dimension(6,6) :: dsde, c0_66, s0_66 ! small strain stiffness real(pReal), dimension(6,6) :: dsde, c0_66, s0_66 ! small strain stiffness
real(pReal), dimension(9,9) :: s_prev99, c_prev99, c0_99, s0_99 ! compliance and stiffness in matrix notation real(pReal), dimension(9,9) :: s_prev99, c_prev99 ! compliance and stiffness in matrix notation
real(pReal), dimension(:,:), allocatable :: s_reduced, c_reduced ! reduced compliance and stiffness (only for stress BC) real(pReal), dimension(:,:), allocatable :: s_reduced, c_reduced ! reduced compliance and stiffness (only for stress BC)
real(pReal), dimension(6,6) :: mask_inversion = reshape([&
1.0_pReal, 1.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal,&
1.0_pReal, 1.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal,&
1.0_pReal, 1.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal,&
0.0_pReal, 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal,&
0.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal,&
0.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal, 1.0_pReal],&
[ 6_pInt, 6_pInt])
real(pReal), dimension(3,3,3,3) :: temp_3333 = 0.0_pReal
integer(pInt) :: size_reduced = 0_pInt ! number of stress BCs integer(pInt) :: size_reduced = 0_pInt ! number of stress BCs
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -151,7 +141,7 @@ program DAMASK_spectral
real(pReal), dimension(3,3) :: xiDyad ! product of wave vectors real(pReal), dimension(3,3) :: xiDyad ! product of wave vectors
real(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat ! gamma operator (field) for spectral method real(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat ! gamma operator (field) for spectral method
real(pReal), dimension(:,:,:,:), allocatable :: xi ! wave vector field for divergence and for gamma operator real(pReal), dimension(:,:,:,:), allocatable :: xi ! wave vector field for divergence and for gamma operator
integer(pInt), dimension(3) :: k_s, cutting_freq integer(pInt), dimension(3) :: k_s
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! loop variables, convergence etc. ! loop variables, convergence etc.
@ -174,7 +164,7 @@ program DAMASK_spectral
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!variables for additional output due to general debugging !variables for additional output due to general debugging
real(pReal) :: defgradDetMax, defgradDetMin, maxCorrectionSym, maxCorrectionSkew, max_diag, max_offdiag real(pReal) :: defgradDetMax, defgradDetMin, maxCorrectionSym, maxCorrectionSkew
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! variables for additional output of divergence calculations ! variables for additional output of divergence calculations
@ -378,8 +368,6 @@ program DAMASK_spectral
res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c) res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c)
Npoints = res(1)*res(2)*res(3) Npoints = res(1)*res(2)*res(3)
wgt = 1.0_pReal/real(Npoints, pReal) wgt = 1.0_pReal/real(Npoints, pReal)
if (cut_off_value <0.0_pReal .or. cut_off_value >0.9_pReal) stop
cutting_freq = nint(real(res,pReal)*cut_off_value,pInt) ! for cut_off_value=0.0 just the highest freq. is removed
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! output of geometry ! output of geometry
@ -394,7 +382,6 @@ program DAMASK_spectral
print '(a,3(i12 ))','resolution a b c:', res print '(a,3(i12 ))','resolution a b c:', res
print '(a,3(f12.5))','dimension x y z:', geomdim print '(a,3(f12.5))','dimension x y z:', geomdim
print '(a,i5)','homogenization: ',homog print '(a,i5)','homogenization: ',homog
if(cut_off_value/=0.0_pReal) print '(a,3(i12),a)', 'cutting away ', cutting_freq, ' frequencies'
print '(a)', '#############################################################' print '(a)', '#############################################################'
print '(a,a)', 'loadcase file: ',trim(getLoadcaseName()) print '(a,a)', 'loadcase file: ',trim(getLoadcaseName())
@ -421,8 +408,8 @@ program DAMASK_spectral
write (*,'(3(3(f12.7,1x)/))',advance='no') merge(math_transpose33(bc(loadcase)%deformation),& write (*,'(3(3(f12.7,1x)/))',advance='no') merge(math_transpose33(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))
write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') ' stress / GPa:',& write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') ' stress / GPa:',&
1e-9*merge(math_transpose33(bc(loadcase)%stress),reshape(spread(DAMASK_NaN,1,9),[ 3,3])& 1e-9_pReal*merge(math_transpose33(bc(loadcase)%stress),&
,transpose(bc(loadcase)%maskStress)) reshape(spread(DAMASK_NaN,1,9),[ 3,3]),transpose(bc(loadcase)%maskStress))
if (any(bc(loadcase)%rotation /= math_I3)) & if (any(bc(loadcase)%rotation /= math_I3)) &
write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') ' rotation of loadframe:',& write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') ' rotation of loadframe:',&
math_transpose33(bc(loadcase)%rotation) math_transpose33(bc(loadcase)%rotation)
@ -473,7 +460,7 @@ program DAMASK_spectral
ielem = ielem + 1_pInt ielem = ielem + 1_pInt
defgrad(i,j,k,1:3,1:3) = math_I3 defgrad(i,j,k,1:3,1:3) = math_I3
defgradold(i,j,k,1:3,1:3) = math_I3 defgradold(i,j,k,1:3,1:3) = math_I3
coordinates(i,j,k,1:3) = geomdim/real(res, pReal)*[i,j,k] - geomdim/real(2_pInt*res,pReal) coordinates(i,j,k,1:3) = geomdim/real(res * [i,j,k], pReal) - geomdim/real(2_pInt*res,pReal)
call CPFEM_general(2_pInt,coordinates(i,j,k,1:3),math_I3,math_I3,temperature(i,j,k),& call CPFEM_general(2_pInt,coordinates(i,j,k,1:3),math_I3,math_I3,temperature(i,j,k),&
0.0_pReal,ielem,1_pInt,cstress,dsde,pstress,dPdF) 0.0_pReal,ielem,1_pInt,cstress,dsde,pstress,dPdF)
c_current = c_current + dPdF c_current = c_current + dPdF
@ -511,7 +498,7 @@ program DAMASK_spectral
do j = 1_pInt, res(2) do j = 1_pInt, res(2)
k_s(2) = j - 1_pInt k_s(2) = j - 1_pInt
if(j > res(2)/2_pInt + 1_pInt) k_s(2) = k_s(2) - res(2) if(j > res(2)/2_pInt + 1_pInt) k_s(2) = k_s(2) - res(2)
do i = 1, res1_red do i = 1_pInt, res1_red
k_s(1) = i - 1_pInt k_s(1) = i - 1_pInt
xi(1:3,i,j,k) = real(k_s, pReal)/geomdim xi(1:3,i,j,k) = real(k_s, pReal)/geomdim
enddo; enddo; enddo enddo; enddo; enddo
@ -528,6 +515,7 @@ program DAMASK_spectral
! i==res(1)/2 .or. i==res(1)/2+2) then ! i==res(1)/2 .or. i==res(1)/2+2) then
! gamma_hat(i,j,k,1:3,1:3,1:3,1:3) = s0_reference ! gamma_hat(i,j,k,1:3,1:3,1:3,1:3) = s0_reference
! else ! else
if(any([i,j,k] /= 1_pInt)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) &
xiDyad(l,m) = xi(l, i,j,k)*xi(m, i,j,k) xiDyad(l,m) = xi(l, i,j,k)*xi(m, i,j,k)
forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) &
@ -535,7 +523,7 @@ program DAMASK_spectral
temp33_Real = math_inv33(temp33_Real) temp33_Real = math_inv33(temp33_Real)
forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, p=1_pInt:3_pInt)& forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, p=1_pInt:3_pInt)&
gamma_hat(i,j,k, l,m,n,p) = temp33_Real(l,n)*xiDyad(m,p) gamma_hat(i,j,k, l,m,n,p) = temp33_Real(l,n)*xiDyad(m,p)
! endif endif
enddo; enddo; enddo enddo; enddo; enddo
gamma_hat(1,1,1, 1:3,1:3,1:3,1:3) = 0.0_pReal ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 gamma_hat(1,1,1, 1:3,1:3,1:3,1:3) = 0.0_pReal ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
endif endif
@ -739,7 +727,7 @@ program DAMASK_spectral
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report begin of new increment ! report begin of new increment
print '(a)', '##################################################################' print '(a)', '##################################################################'
print '(A,I5.5,A,es12.6)', 'Increment ', totalIncsCounter, ' Time ',time print '(A,I5.5,A,es12.5)', 'Increment ', totalIncsCounter, ' Time ',time
guessmode = 1.0_pReal ! keep guessing along former trajectory during same loadcase guessmode = 1.0_pReal ! keep guessing along former trajectory during same loadcase
CPFEM_mode = 1_pInt ! winding forward CPFEM_mode = 1_pInt ! winding forward
@ -801,7 +789,7 @@ program DAMASK_spectral
row = (mod(totalIncsCounter+iter-2_pInt,9_pInt))/3_pInt + 1_pInt ! go through the elements of the tensors, controlled by totalIncsCounter and iter, starting at 1 row = (mod(totalIncsCounter+iter-2_pInt,9_pInt))/3_pInt + 1_pInt ! go through the elements of the tensors, controlled by totalIncsCounter and iter, starting at 1
column = (mod(totalIncsCounter+iter-2_pInt,3_pInt)) + 1_pInt column = (mod(totalIncsCounter+iter-2_pInt,3_pInt)) + 1_pInt
scalarField_real(1:res(1),1:res(2),1:res(3)) =& ! store the selected component scalarField_real(1:res(1),1:res(2),1:res(3)) =& ! store the selected component
tensorField_real(1:res(1),1:res(2),1:res(3),row,column) cmplx(tensorField_real(1:res(1),1:res(2),1:res(3),row,column),0.0_pReal,pReal)
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -829,7 +817,7 @@ program DAMASK_spectral
pstress_av_lab = real(tensorField_fourier(1,1,1,1:3,1:3),pReal)*wgt pstress_av_lab = real(tensorField_fourier(1,1,1,1:3,1:3),pReal)*wgt
pstress_av = math_rotate_forward33(pstress_av_lab,bc(loadcase)%rotation) pstress_av = math_rotate_forward33(pstress_av_lab,bc(loadcase)%rotation)
write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') 'Piola-Kirchhoff stress / MPa:',& write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') 'Piola-Kirchhoff stress / MPa:',&
math_transpose33(pstress_av)/1.e6 math_transpose33(pstress_av)/1.e6_pReal
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! comparing 1 and 3x3 FT results ! comparing 1 and 3x3 FT results
@ -940,17 +928,14 @@ program DAMASK_spectral
print '(a,es11.4)', 'error divergence FT max = ',err_div_max print '(a,es11.4)', 'error divergence FT max = ',err_div_max
print '(a,es11.4)', 'error divergence Real RMS = ',err_real_div_RMS print '(a,es11.4)', 'error divergence Real RMS = ',err_real_div_RMS
print '(a,es11.4)', 'error divergence Real max = ',err_real_div_max print '(a,es11.4)', 'error divergence Real max = ',err_real_div_max
print '(a,es11.4)', 'divergence RMS FT/real = ',err_div_RMS/err_real_div_RMS
print '(a,es11.4)', 'divergence max FT/real = ',err_div_max/err_real_div_max
print '(a,es11.4)', 'max deviat. from postProc = ',max_div_error print '(a,es11.4)', 'max deviat. from postProc = ',max_div_error
endif endif
print '(a,f6.2,a,es11.4,a)', 'error divergence = ', err_div/err_div_tol, ' (',err_div,' 1/m)' print '(a,f6.2,a,es11.4,3a)','error divergence = ', err_div/err_div_tol, &
' (',err_div_RMS,' N/m',char(179),')'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! divergence is calculated from FT(stress), depending on algorithm use field for spectral method ! divergence is calculated from FT(stress), depending on algorithm use field for spectral method
if (.not. simplified_algorithm) tensorField_fourier = tau_fourier if (.not. simplified_algorithm) tensorField_fourier = tau_fourier
max_diag = tiny(1.0_pReal)
max_offdiag = tiny(1.0_pReal)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! to the actual spectral method calculation (mechanical equilibrium) ! to the actual spectral method calculation (mechanical equilibrium)
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
@ -962,6 +947,7 @@ program DAMASK_spectral
! forall( m = 1_pInt:3_pInt, n = 1_pInt:3_pInt)& ! forall( m = 1_pInt:3_pInt, n = 1_pInt:3_pInt)&
! temp33_Complex(m,n) = sum(s0_reference(m,n, 1:3,1:3)* tensorField_fourier(i,j,k,1:3,1:3)) ! temp33_Complex(m,n) = sum(s0_reference(m,n, 1:3,1:3)* tensorField_fourier(i,j,k,1:3,1:3))
! else ! else
if(any([i,j,k] /= 1_pInt)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) &
xiDyad(l,m) = xi(l, i,j,k)*xi(m, i,j,k) xiDyad(l,m) = xi(l, i,j,k)*xi(m, i,j,k)
forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) &
@ -970,26 +956,29 @@ program DAMASK_spectral
forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, p=1_pInt:3_pInt)& forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, p=1_pInt:3_pInt)&
gamma_hat(1,1,1, l,m,n,p) = temp33_Real(l,n)*xiDyad(m,p) gamma_hat(1,1,1, l,m,n,p) = temp33_Real(l,n)*xiDyad(m,p)
forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) &
temp33_Complex(l,m) = sum(gamma_hat(1,1,1, l,m, 1:3,1:3) * tensorField_fourier(i,j,k,1:3,1:3)) temp33_Complex(l,m) = sum(gamma_hat(1,1,1, l,m, 1:3,1:3) *&
tensorField_fourier(i,j,k,1:3,1:3))
tensorField_fourier(i,j,k,1:3,1:3) = temp33_Complex tensorField_fourier(i,j,k,1:3,1:3) = temp33_Complex
! endif endif
enddo; enddo; enddo enddo; enddo; enddo
else ! use precalculated gamma-operator else ! use precalculated gamma-operator
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt,res1_red do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt,res1_red
forall( m = 1_pInt:3_pInt, n = 1_pInt:3_pInt) & forall( m = 1_pInt:3_pInt, n = 1_pInt:3_pInt) &
temp33_Complex(m,n) = sum(gamma_hat(i,j,k, m,n, 1:3,1:3) * tensorField_fourier(i,j,k,1:3,1:3)) temp33_Complex(m,n) = sum(gamma_hat(i,j,k, m,n, 1:3,1:3) *&
tensorField_fourier(i,j,k,1:3,1:3))
tensorField_fourier(i,j,k, 1:3,1:3) = temp33_Complex tensorField_fourier(i,j,k, 1:3,1:3) = temp33_Complex
enddo; enddo; enddo enddo; enddo; enddo
endif endif
if (simplified_algorithm) then ! do not use the polarization field based algorithm if (simplified_algorithm) then ! do not use the polarization field based algorithm
tensorField_fourier(1,1,1,1:3,1:3) = (defgrad_av_lab - defgradAim_lab) & ! assign (negative) average deformation gradient change to zero frequency (real part) tensorField_fourier(1,1,1,1:3,1:3) = cmplx((defgrad_av_lab - defgradAim_lab) & ! assign (negative) average deformation gradient change to zero frequency (real part)
* real(Npoints,pReal) ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 * real(Npoints,pReal),0.0_pReal,pReal) ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
else else
tensorField_fourier(1,1,1,1:3,1:3) = defgradAim_lab * real(Npoints,pReal) ! assign deformation aim to zero frequency (real part) tensorField_fourier(1,1,1,1:3,1:3) = cmplx(defgradAim_lab*real(Npoints,pReal),& ! assign deformation aim to zero frequency (real part)
0.0_pReal,pReal)
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -998,7 +987,7 @@ program DAMASK_spectral
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red
scalarField_fourier(i,j,k) = tensorField_fourier(i,j,k,row,column) scalarField_fourier(i,j,k) = tensorField_fourier(i,j,k,row,column)
enddo; enddo; enddo enddo; enddo; enddo
do i = 0_pInt, res(1)/2_pInt-2_pInt !unpack fft data for conj complex symmetric part. can be directly used in calculation of cstress_field do i = 0_pInt, res(1)/2_pInt-2_pInt !unpack fft data for conj complex symmetric part
m = 1_pInt m = 1_pInt
do k = 1_pInt, res(3) do k = 1_pInt, res(3)
n = 1_pInt n = 1_pInt
@ -1136,11 +1125,13 @@ end program DAMASK_spectral
! quit subroutine to satisfy IO_error ! quit subroutine to satisfy IO_error
! !
!******************************************************************** !********************************************************************
subroutine quit(id) subroutine quit(stop_id)
use prec use prec
implicit none implicit none
integer(pInt) id integer(pInt), intent(in) :: stop_id
print*, stop_id
stop 'abnormal termination of DAMASK_spectral'
stop
end subroutine end subroutine

78
code/DAMASK_spectral_interface.f90
View File

@ -39,7 +39,7 @@ subroutine DAMASK_interface_init()
implicit none implicit none
character(len=1024) commandLine, hostName, userName character(len=1024) commandLine, hostName, userName
integer(pInt):: i, start = 0_pInt, length=0_pInt integer :: i, start = 0, length=0
integer, dimension(8) :: date_and_time_values ! type default integer integer, dimension(8) :: date_and_time_values ! type default integer
call get_command(commandLine) call get_command(commandLine)
call DATE_AND_TIME(VALUES=date_and_time_values) call DATE_AND_TIME(VALUES=date_and_time_values)
@ -47,7 +47,7 @@ subroutine DAMASK_interface_init()
if(64<iachar(commandLine(i:i)) .and. iachar(commandLine(i:i))<91) commandLine(i:i) =achar(iachar(commandLine(i:i))+32) if(64<iachar(commandLine(i:i)) .and. iachar(commandLine(i:i))<91) commandLine(i:i) =achar(iachar(commandLine(i:i))+32)
enddo enddo
if(index(commandLine,' -h ',.true.)>0_pInt .or. index(commandLine,' --help ',.true.)>0_pInt) then ! search for ' -h ' or '--help' if(index(commandLine,' -h ',.true.)>0 .or. index(commandLine,' --help ',.true.)>0) then ! search for ' -h ' or '--help'
write(6,*) '$Id$' write(6,*) '$Id$'
#include "compilation_info.f90" #include "compilation_info.f90"
print '(a)', '#############################################################' print '(a)', '#############################################################'
@ -87,12 +87,12 @@ subroutine DAMASK_interface_init()
endif endif
if (.not.(command_argument_count()==4 .or. command_argument_count()==6)) & ! check for correct number of given arguments (no --help) if (.not.(command_argument_count()==4 .or. command_argument_count()==6)) & ! check for correct number of given arguments (no --help)
stop 'Wrong Nr. of Arguments. Run DAMASK_spectral.exe --help' ! Could not find valid keyword (position 0 +3). Functions from IO.f90 are not available stop 'Wrong Nr. of Arguments. Run DAMASK_spectral.exe --help' ! Could not find valid keyword (position 0 +3). Functions from IO.f90 are not available
start = index(commandLine,'-g',.true.) + 3_pInt ! search for '-g' and jump to first char of geometry start = index(commandLine,'-g',.true.) + 3 ! search for '-g' and jump to first char of geometry
if (index(commandLine,'--geom',.true.)>0) then ! if '--geom' is found, use that (contains '-g') if (index(commandLine,'--geom',.true.)>0) then ! if '--geom' is found, use that (contains '-g')
start = index(commandLine,'--geom',.true.) + 7_pInt start = index(commandLine,'--geom',.true.) + 7
endif endif
if (index(commandLine,'--geometry',.true.)>0) then ! again, now searching for --geometry' if (index(commandLine,'--geometry',.true.)>0) then ! again, now searching for --geometry'
start = index(commandLine,'--geometry',.true.) + 11_pInt start = index(commandLine,'--geometry',.true.) + 11
endif endif
if(start==3_pInt) stop 'No Geometry specified, terminating DAMASK'! Could not find valid keyword (position 0 +3). Functions from IO.f90 are not available if(start==3_pInt) stop 'No Geometry specified, terminating DAMASK'! Could not find valid keyword (position 0 +3). Functions from IO.f90 are not available
length = index(commandLine(start:len(commandLine)),' ',.false.) length = index(commandLine(start:len(commandLine)),' ',.false.)
@ -105,12 +105,12 @@ subroutine DAMASK_interface_init()
if(64<iachar(commandLine(i:i)) .and. iachar(commandLine(i:i))<91) commandLine(i:i) =achar(iachar(commandLine(i:i))+32) if(64<iachar(commandLine(i:i)) .and. iachar(commandLine(i:i))<91) commandLine(i:i) =achar(iachar(commandLine(i:i))+32)
enddo enddo
start = index(commandLine,'-l',.true.) + 3_pInt ! search for '-l' and jump forward iby 3 to given name start = index(commandLine,'-l',.true.) + 3 ! search for '-l' and jump forward iby 3 to given name
if (index(commandLine,'--load',.true.)>0) then ! if '--load' is found, use that (contains '-l') if (index(commandLine,'--load',.true.)>0) then ! if '--load' is found, use that (contains '-l')
start = index(commandLine,'--load',.true.) + 7_pInt start = index(commandLine,'--load',.true.) + 7
endif endif
if (index(commandLine,'--loadcase',.true.)>0) then ! again, now searching for --loadcase' if (index(commandLine,'--loadcase',.true.)>0) then ! again, now searching for --loadcase'
start = index(commandLine,'--loadcase',.true.) + 11_pInt start = index(commandLine,'--loadcase',.true.) + 11
endif endif
if(start==3_pInt) stop 'No Loadcase specified, terminating DAMASK'! Could not find valid keyword (position 0 +3). Functions from IO.f90 are not available if(start==3_pInt) stop 'No Loadcase specified, terminating DAMASK'! Could not find valid keyword (position 0 +3). Functions from IO.f90 are not available
length = index(commandLine(start:len(commandLine)),' ',.false.) length = index(commandLine(start:len(commandLine)),' ',.false.)
@ -123,9 +123,9 @@ subroutine DAMASK_interface_init()
if(64<iachar(commandLine(i:i)) .and. iachar(commandLine(i:i))<91) commandLine(i:i) =achar(iachar(commandLine(i:i))+32) if(64<iachar(commandLine(i:i)) .and. iachar(commandLine(i:i))<91) commandLine(i:i) =achar(iachar(commandLine(i:i))+32)
enddo enddo
start = index(commandLine,'-r',.true.) + 3_pInt ! search for '-r' and jump forward iby 3 to given name start = index(commandLine,'-r',.true.) + 3 ! search for '-r' and jump forward iby 3 to given name
if (index(commandLine,'--restart',.true.)>0) then ! if '--restart' is found, use that (contains '-l') if (index(commandLine,'--restart',.true.)>0) then ! if '--restart' is found, use that (contains '-l')
start = index(commandLine,'--restart',.true.) + 7_pInt start = index(commandLine,'--restart',.true.) + 7
endif endif
length = index(commandLine(start:len(commandLine)),' ',.false.) length = index(commandLine(start:len(commandLine)),' ',.false.)
@ -201,12 +201,12 @@ function getModelName()
character(1024) getModelName, cwd character(1024) getModelName, cwd
character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forwardslash, backwardslash character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forwardslash, backwardslash
integer(pInt) :: posExt,posSep integer :: posExt,posSep
posExt = scan(geometryParameter,'.',back=.true.) posExt = scan(geometryParameter,'.',back=.true.)
posSep = scan(geometryParameter,pathSep,back=.true.) posSep = scan(geometryParameter,pathSep,back=.true.)
if (posExt <= posSep) posExt = len_trim(geometryParameter)+1_pInt ! no extension present if (posExt <= posSep) posExt = len_trim(geometryParameter)+1 ! no extension present
getModelName = geometryParameter(1:posExt-1_pInt) ! path to geometry file (excl. extension) getModelName = geometryParameter(1:posExt-1_pInt) ! path to geometry file (excl. extension)
if (scan(getModelName,pathSep) /= 1) then ! relative path given as command line argument if (scan(getModelName,pathSep) /= 1) then ! relative path given as command line argument
@ -227,19 +227,17 @@ endfunction getModelName
!******************************************************************** !********************************************************************
function getLoadCase() function getLoadCase()
use prec, only: pInt
implicit none implicit none
character(1024) getLoadCase character(1024) :: getLoadCase
character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forwardslash, backwardslash character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forwardslash, backwardslash
integer(pInt) posExt,posSep integer :: posExt,posSep
posExt = scan(loadcaseParameter,'.',back=.true.) posExt = scan(loadcaseParameter,'.',back=.true.)
posSep = scan(loadcaseParameter,pathSep,back=.true.) posSep = scan(loadcaseParameter,pathSep,back=.true.)
if (posExt <= posSep) posExt = len_trim(loadcaseParameter)+1_pInt ! no extension present if (posExt <= posSep) posExt = len_trim(loadcaseParameter)+1 ! no extension present
getLoadCase = loadcaseParameter(posSep+1_pInt:posExt-1_pInt) ! name of load case file exluding extension getLoadCase = loadcaseParameter(posSep+1:posExt-1) ! name of load case file exluding extension
endfunction getLoadCase endfunction getLoadCase
@ -254,11 +252,9 @@ function getLoadcaseName()
implicit none implicit none
character(len=1024) getLoadcaseName,cwd character(len=1024) :: getLoadcaseName,cwd
character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forwardslash, backwardslash character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forwardslash, backwardslash
integer(pInt) posExt,posSep integer :: posExt = 0, posSep
posExt = 0_pInt
getLoadcaseName = loadcaseParameter getLoadcaseName = loadcaseParameter
posExt = scan(getLoadcaseName,'.',back=.true.) posExt = scan(getLoadcaseName,'.',back=.true.)
posSep = scan(getLoadcaseName,pathSep,back=.true.) posSep = scan(getLoadcaseName,pathSep,back=.true.)
@ -286,31 +282,31 @@ function rectifyPath(path)
implicit none implicit none
character(len=*) path character(len=*) :: path
character(len=len_trim(path)) rectifyPath character(len=len_trim(path)) :: rectifyPath
integer(pInt) i,j,k,l integer :: i,j,k,l !no pInt
!remove ./ from path !remove ./ from path
l = len_trim(path) l = len_trim(path)
rectifyPath = path rectifyPath = path
do i = l,3_pInt,-1_pInt do i = l,3,-1
if ( rectifyPath(i-1_pInt:i) == './' .and. rectifyPath(i-2_pInt:i-2_pInt) /= '.' ) & if ( rectifyPath(i-1:i) == './' .and. rectifyPath(i-2:i-2) /= '.' ) &
rectifyPath(i-1_pInt:l) = rectifyPath(i+1_pInt:l)//' ' rectifyPath(i-1:l) = rectifyPath(i+1:l)//' '
enddo enddo
!remove ../ and corresponding directory from rectifyPath !remove ../ and corresponding directory from rectifyPath
l = len_trim(rectifyPath) l = len_trim(rectifyPath)
i = index(rectifyPath(i:l),'../') i = index(rectifyPath(i:l),'../')
j = 0_pInt j = 0
do while (i > j) do while (i > j)
j = scan(rectifyPath(1:i-2_pInt),'/',back=.true.) j = scan(rectifyPath(1:i-2),'/',back=.true.)
rectifyPath(j+1_pInt:l) = rectifyPath(i+3_pInt:l)//repeat(' ',2_pInt+i-j) rectifyPath(j+1:l) = rectifyPath(i+3:l)//repeat(' ',2+i-j)
if (rectifyPath(j+1_pInt:j+1_pInt) == '/') then !search for '//' that appear in case of XXX/../../XXX if (rectifyPath(j+1:j+1) == '/') then !search for '//' that appear in case of XXX/../../XXX
k = len_trim(rectifyPath) k = len_trim(rectifyPath)
rectifyPath(j+1_pInt:k-1_pInt) = rectifyPath(j+2_pInt:k) rectifyPath(j+1:k-1) = rectifyPath(j+2:k)
rectifyPath(k:k) = ' ' rectifyPath(k:k) = ' '
endif endif
i = j+index(rectifyPath(j+1_pInt:l),'../') i = j+index(rectifyPath(j+1:l),'../')
enddo enddo
if(len_trim(rectifyPath) == 0) rectifyPath = '/' if(len_trim(rectifyPath) == 0) rectifyPath = '/'
@ -330,18 +326,18 @@ function makeRelativePath(a,b)
character (len=*) :: a,b character (len=*) :: a,b
character (len=1024) :: makeRelativePath character (len=1024) :: makeRelativePath
integer(pInt) i,posLastCommonSlash,remainingSlashes integer :: i,posLastCommonSlash,remainingSlashes !no pInt
posLastCommonSlash = 0_pInt posLastCommonSlash = 0
remainingSlashes = 0_pInt remainingSlashes = 0
do i = 1_pInt,min(1024,len_trim(a),len_trim(b)) do i = 1, min(1024,len_trim(a),len_trim(b))
if (a(i:i) /= b(i:i)) exit if (a(i:i) /= b(i:i)) exit
if (a(i:i) == '/') posLastCommonSlash = i if (a(i:i) == '/') posLastCommonSlash = i
enddo enddo
do i = posLastCommonSlash+1_pInt,len_trim(a) do i = posLastCommonSlash+1,len_trim(a)
if (a(i:i) == '/') remainingSlashes = remainingSlashes + 1_pInt if (a(i:i) == '/') remainingSlashes = remainingSlashes + 1
enddo enddo
makeRelativePath = repeat('../',remainingSlashes)//b(posLastCommonSlash+1_pInt:len_trim(b)) makeRelativePath = repeat('../',remainingSlashes)//b(posLastCommonSlash+1:len_trim(b))
endfunction makeRelativePath endfunction makeRelativePath

12
code/FEsolving.f90
View File

@ -70,9 +70,9 @@
commandLine(i:i) = achar(iachar(commandLine(i:i))+32) ! make lowercase commandLine(i:i) = achar(iachar(commandLine(i:i))+32) ! make lowercase
enddo enddo
if (index(commandLine,'-r ',.true.)>0) & ! look for -r if (index(commandLine,'-r ',.true.)>0) & ! look for -r
start = index(commandLine,'-r ',.true.) + 3_pInt ! set to position after trailing space start = index(commandLine,'-r ',.true.) + 3 ! set to position after trailing space
if (index(commandLine,'--restart ',.true.)>0) & ! look for --restart if (index(commandLine,'--restart ',.true.)>0) & ! look for --restart
start = index(commandLine,'--restart ',.true.) + 10_pInt ! set to position after trailing space start = index(commandLine,'--restart ',.true.) + 10 ! set to position after trailing space
if(start /= 0_pInt) then ! found something if(start /= 0_pInt) then ! found something
length = verify(commandLine(start:len(commandLine)),'0123456789',.false.) ! where is first non number after argument? length = verify(commandLine(start:len(commandLine)),'0123456789',.false.) ! where is first non number after argument?
read(commandLine(start:start+length),'(I12)') restartInc ! read argument read(commandLine(start:start+length),'(I12)') restartInc ! read argument
@ -99,12 +99,12 @@
restartWrite = iand(IO_intValue(line,positions,1_pInt),1_pInt) > 0_pInt restartWrite = iand(IO_intValue(line,positions,1_pInt),1_pInt) > 0_pInt
restartRead = iand(IO_intValue(line,positions,1_pInt),2_pInt) > 0_pInt restartRead = iand(IO_intValue(line,positions,1_pInt),2_pInt) > 0_pInt
case ('*restart') case ('*restart')
do i=2,positions(1) do j=2_pInt,positions(1)
restartWrite = (IO_lc(IO_StringValue(line,positions,i)) == 'write') .or. restartWrite restartWrite = (IO_lc(IO_StringValue(line,positions,j)) == 'write') .or. restartWrite
restartRead = (IO_lc(IO_StringValue(line,positions,i)) == 'read') .or. restartRead restartRead = (IO_lc(IO_StringValue(line,positions,j)) == 'read') .or. restartRead
enddo enddo
if(restartWrite) then if(restartWrite) then
do j=2,positions(1) do j=2_pInt,positions(1)
restartWrite = (IO_lc(IO_StringValue(line,positions,j)) /= 'frequency=0') .and. restartWrite restartWrite = (IO_lc(IO_StringValue(line,positions,j)) /= 'frequency=0') .and. restartWrite
enddo enddo
endif endif

432
code/IO.f90
View File

@ -25,20 +25,20 @@
CONTAINS CONTAINS
!--------------------------- !---------------------------
! function IO_abaqus_assembleInputFile ! function IO_abaqus_assembleInputFile
! subroutine IO_open_file(unit,relPath) ! subroutine IO_open_file(myUnit,relPath)
! subroutine IO_open_inputFile(unit, model) ! subroutine IO_open_inputFile(myUnit, model)
! subroutine IO_open_logFile(unit) ! subroutine IO_open_logFile(myUnit)
! function IO_hybridIA(Nast,ODFfileName) ! function IO_hybridIA(Nast,ODFfileName)
! private function hybridIA_reps(dV_V,steps,C) ! private function hybridIA_reps(dV_V,steps,C)
! function IO_stringPos(line,maxN) ! function IO_stringPos(line,maxN)
! function IO_stringValue(line,positions,pos) ! function IO_stringValue(line,positions,myPos)
! function IO_floatValue(line,positions,pos) ! function IO_floatValue(line,positions,myPos)
! function IO_intValue(line,positions,pos) ! function IO_intValue(line,positions,myPos)
! function IO_fixedStringValue(line,ends,pos) ! function IO_fixedStringValue(line,ends,myPos)
! function IO_fixedFloatValue(line,ends,pos) ! function IO_fixedFloatValue(line,ends,myPos)
! function IO_fixedFloatNoEValue(line,ends,pos) ! function IO_fixedFloatNoEValue(line,ends,myPos)
! function IO_fixedIntValue(line,ends,pos) ! function IO_fixedIntValue(line,ends,myPos)
! function IO_continousIntValues(unit,maxN) ! function IO_continousIntValues(myUnit,maxN)
! function IO_lc(line) ! function IO_lc(line)
! subroutine IO_lcInplace(line) ! subroutine IO_lcInplace(line)
! subroutine IO_error(ID) ! subroutine IO_error(ID)
@ -76,7 +76,7 @@ recursive function IO_abaqus_assembleInputFile(unit1,unit2) result(createSuccess
character(len=300) line,fname character(len=300) line,fname
integer(pInt), intent(in) :: unit1, unit2 integer(pInt), intent(in) :: unit1, unit2
logical createSuccess,fexist logical createSuccess,fexist
integer(pInt), parameter :: maxNchunks = 6 integer(pInt), parameter :: maxNchunks = 6_pInt
integer(pInt), dimension(1+2*maxNchunks) :: positions integer(pInt), dimension(1+2*maxNchunks) :: positions
@ -88,7 +88,7 @@ recursive function IO_abaqus_assembleInputFile(unit1,unit2) result(createSuccess
! call IO_lcInPlace(line) ! call IO_lcInPlace(line)
if (IO_lc(IO_StringValue(line,positions,1_pInt))=='*include') then if (IO_lc(IO_StringValue(line,positions,1_pInt))=='*include') then
fname = trim(getSolverWorkingDirectoryName())//trim(line(9_pInt+scan(line(9_pInt:),'='):)) fname = trim(getSolverWorkingDirectoryName())//trim(line(9+scan(line(9:),'='):))
inquire(file=fname, exist=fexist) inquire(file=fname, exist=fexist)
if (.not.(fexist)) then if (.not.(fexist)) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
@ -121,25 +121,25 @@ end function
!*********************************************************** !***********************************************************
! check if the input file for Abaqus contains part info ! check if the input file for Abaqus contains part info
!*********************************************************** !***********************************************************
function IO_abaqus_hasNoPart(unit) function IO_abaqus_hasNoPart(myUnit)
use prec, only: pInt use prec, only: pInt
implicit none implicit none
integer(pInt) unit integer(pInt) myUnit
integer(pInt), parameter :: maxNchunks = 1 integer(pInt), parameter :: maxNchunks = 1_pInt
integer(pInt), dimension(1+2*maxNchunks) :: pos integer(pInt), dimension(1+2*maxNchunks) :: myPos
logical IO_abaqus_hasNoPart logical IO_abaqus_hasNoPart
character(len=300) line character(len=300) line
IO_abaqus_hasNoPart = .true. IO_abaqus_hasNoPart = .true.
610 FORMAT(A300) 610 FORMAT(A300)
rewind(unit) rewind(myUnit)
do do
read(unit,610,END=620) line read(myUnit,610,END=620) line
pos = IO_stringPos(line,maxNchunks) myPos = IO_stringPos(line,maxNchunks)
if (IO_lc(IO_stringValue(line,pos,1_pInt)) == '*part' ) then if (IO_lc(IO_stringValue(line,myPos,1_pInt)) == '*part' ) then
IO_abaqus_hasNoPart = .false. IO_abaqus_hasNoPart = .false.
exit exit
endif endif
@ -150,22 +150,22 @@ end function
!******************************************************************** !********************************************************************
! open existing file to given unit ! open existing file to given myUnit
! path to file is relative to working directory ! path to file is relative to working directory
!******************************************************************** !********************************************************************
logical function IO_open_file_stat(unit,relPath) logical function IO_open_file_stat(myUnit,relPath)
use prec, only: pInt use prec, only: pInt
use DAMASK_interface use DAMASK_interface
implicit none implicit none
integer(pInt), intent(in) :: unit integer(pInt), intent(in) :: myUnit
character(len=*), intent(in) :: relPath character(len=*), intent(in) :: relPath
character(len=1024) path character(len=1024) path
integer(pInt) stat integer(pInt) stat
path = trim(getSolverWorkingDirectoryName())//relPath path = trim(getSolverWorkingDirectoryName())//relPath
open(unit,status='old',iostat=stat,file=path) open(myUnit,status='old',iostat=stat,file=path)
IO_open_file_stat = (stat == 0_pInt) IO_open_file_stat = (stat == 0_pInt)
endfunction endfunction
@ -175,37 +175,37 @@ end function
! open existing file to given unit ! open existing file to given unit
! path to file is relative to working directory ! path to file is relative to working directory
!******************************************************************** !********************************************************************
subroutine IO_open_file(unit,relPath) subroutine IO_open_file(myUnit,relPath)
use prec, only: pInt use prec, only: pInt
use DAMASK_interface use DAMASK_interface
implicit none implicit none
integer(pInt), intent(in) :: unit integer(pInt), intent(in) ::myUnit
character(len=*), intent(in) :: relPath character(len=*), intent(in) :: relPath
character(len=1024) path character(len=1024) path
integer(pInt) stat integer(pInt) stat
path = trim(getSolverWorkingDirectoryName())//relPath path = trim(getSolverWorkingDirectoryName())//relPath
open(unit,status='old',iostat=stat,file=path) open(myUnit,status='old',iostat=stat,file=path)
if (stat /= 0_pInt) call IO_error(100_pInt,ext_msg=path) if (stat /= 0_pInt) call IO_error(100_pInt,ext_msg=path)
endsubroutine endsubroutine
!******************************************************************** !********************************************************************
! open FEM inputfile to given unit ! open FEM inputfile to given myUnit
! AP: 12.07.10 ! AP: 12.07.10
! : changed the function to open *.inp_assembly, which is basically ! : changed the function to open *.inp_assembly, which is basically
! the input file without comment lines and possibly assembled includes ! the input file without comment lines and possibly assembled includes
!******************************************************************** !********************************************************************
subroutine IO_open_inputFile(unit,model) subroutine IO_open_inputFile(myUnit,model)
use prec, only: pReal, pInt use prec, only: pReal, pInt
use DAMASK_interface use DAMASK_interface
implicit none implicit none
integer(pInt), intent(in) :: unit integer(pInt), intent(in) :: myUnit
character(len=*), intent(in) :: model character(len=*), intent(in) :: model
character(len=1024) path character(len=1024) path
integer(pInt) stat integer(pInt) stat
@ -214,19 +214,18 @@ end function
if (FEsolver == 'Abaqus') then if (FEsolver == 'Abaqus') then
path = trim(getSolverWorkingDirectoryName())//trim(model)//InputFileExtension path = trim(getSolverWorkingDirectoryName())//trim(model)//InputFileExtension
open(unit+1,status='old',iostat=stat,file=path) open(myUnit+1,status='old',iostat=stat,file=path)
if (stat /= 0_pInt) call IO_error(100_pInt,ext_msg=path) if (stat /= 0_pInt) call IO_error(100_pInt,ext_msg=path)
path = trim(getSolverWorkingDirectoryName())//trim(model)//InputFileExtension//'_assembly' path = trim(getSolverWorkingDirectoryName())//trim(model)//InputFileExtension//'_assembly'
open(unit,iostat=stat,file=path) open(myUnit,iostat=stat,file=path)
if (stat /= 0_pInt) call IO_error(100_pInt,ext_msg=path) if (stat /= 0_pInt) call IO_error(100_pInt,ext_msg=path)
if (IO_abaqus_assembleInputFile(unit,unit+1_pInt)) call IO_error(103_pInt) ! strip comments and concatenate any "include"s if (IO_abaqus_assembleInputFile(myUnit,myUnit+1_pInt)) call IO_error(103_pInt) ! strip comments and concatenate any "include"s
close(unit+1_pInt) close(myUnit+1_pInt)
else else
path = trim(getSolverWorkingDirectoryName())//trim(model)//InputFileExtension path = trim(getSolverWorkingDirectoryName())//trim(model)//InputFileExtension
open(unit,status='old',iostat=stat,file=path) open(myUnit,status='old',iostat=stat,file=path)
if (stat /= 0_pInt) call IO_error(100_pInt,ext_msg=path) if (stat /= 0_pInt) call IO_error(100_pInt,ext_msg=path)
endif endif
@ -235,20 +234,19 @@ end function
!******************************************************************** !********************************************************************
! open FEM logfile to given unit ! open FEM logfile to given myUnit
!******************************************************************** !********************************************************************
subroutine IO_open_logFile(unit) subroutine IO_open_logFile(myUnit)
use prec, only: pReal, pInt use prec, only: pReal, pInt
use DAMASK_interface use DAMASK_interface
implicit none implicit none
integer(pInt), intent(in) :: unit integer(pInt), intent(in) :: myUnit
character(len=1024) path character(len=1024) path
integer(pInt) stat integer(pInt) stat
path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//LogFileExtension path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//LogFileExtension
open(unit,status='old',iostat=stat,file=path) open(myUnit,status='old',iostat=stat,file=path)
if (stat /= 0) call IO_error(100_pInt,ext_msg=path) if (stat /= 0) call IO_error(100_pInt,ext_msg=path)
endsubroutine endsubroutine
@ -256,21 +254,21 @@ end function
!******************************************************************** !********************************************************************
! open (write) file related to current job ! open (write) file related to current job
! but with different extension to given unit ! but with different extension to given myUnit
!******************************************************************** !********************************************************************
logical function IO_open_jobFile_stat(unit,newExt) logical function IO_open_jobFile_stat(myUnit,newExt)
use prec, only: pReal, pInt use prec, only: pReal, pInt
use DAMASK_interface use DAMASK_interface
implicit none implicit none
integer(pInt), intent(in) :: unit integer(pInt), intent(in) :: myUnit
character(len=*), intent(in) :: newExt character(len=*), intent(in) :: newExt
character(len=1024) path character(len=1024) path
integer(pInt) stat integer(pInt) stat
path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//newExt path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//newExt
open(unit,status='old',iostat=stat,file=path) open(myUnit,status='old',iostat=stat,file=path)
IO_open_jobFile_stat = (stat == 0_pInt) IO_open_jobFile_stat = (stat == 0_pInt)
endfunction endfunction
@ -280,19 +278,19 @@ end function
! open (write) file related to current job ! open (write) file related to current job
! but with different extension to given unit ! but with different extension to given unit
!******************************************************************** !********************************************************************
subroutine IO_open_jobFile(unit,newExt) subroutine IO_open_jobFile(myUnit,newExt)
use prec, only: pReal, pInt use prec, only: pReal, pInt
use DAMASK_interface use DAMASK_interface
implicit none implicit none
integer(pInt), intent(in) :: unit integer(pInt), intent(in) :: myUnit
character(len=*), intent(in) :: newExt character(len=*), intent(in) :: newExt
character(len=1024) path character(len=1024) path
integer(pInt) stat integer(pInt) stat
path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//newExt path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//newExt
open(unit,status='old',iostat=stat,file=path) open(myUnit,status='old',iostat=stat,file=path)
if (stat /= 0_pInt) call IO_error(100_pInt,ext_msg=path) if (stat /= 0_pInt) call IO_error(100_pInt,ext_msg=path)
endsubroutine endsubroutine
@ -300,21 +298,21 @@ end function
!******************************************************************** !********************************************************************
! open (write) file related to current job ! open (write) file related to current job
! but with different extension to given unit ! but with different extension to given myUnit
!******************************************************************** !********************************************************************
subroutine IO_write_jobFile(unit,newExt) subroutine IO_write_jobFile(myUnit,newExt)
use prec, only: pReal, pInt use prec, only: pReal, pInt
use DAMASK_interface use DAMASK_interface
implicit none implicit none
integer(pInt), intent(in) :: unit integer(pInt), intent(in) :: myUnit
character(len=*), intent(in) :: newExt character(len=*), intent(in) :: newExt
character(len=1024) path character(len=1024) path
integer(pInt) stat integer(pInt) stat
path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//newExt path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//newExt
open(unit,status='replace',iostat=stat,file=path) open(myUnit,status='replace',iostat=stat,file=path)
if (stat /= 0_pInt) call IO_error(100_pInt,ext_msg=path) if (stat /= 0_pInt) call IO_error(100_pInt,ext_msg=path)
endsubroutine endsubroutine
@ -322,15 +320,15 @@ end function
!******************************************************************** !********************************************************************
! open (write) binary file related to current job ! open (write) binary file related to current job
! but with different extension to given unit ! but with different extension to given myUnit
!******************************************************************** !********************************************************************
subroutine IO_write_jobBinaryFile(unit,newExt,recMultiplier) subroutine IO_write_jobBinaryFile(myUnit,newExt,recMultiplier)
use prec, only: pReal, pInt use prec, only: pReal, pInt
use DAMASK_interface use DAMASK_interface
implicit none implicit none
integer(pInt), intent(in) :: unit integer(pInt), intent(in) :: myUnit
integer(pInt), intent(in), optional :: recMultiplier integer(pInt), intent(in), optional :: recMultiplier
character(len=*), intent(in) :: newExt character(len=*), intent(in) :: newExt
character(len=1024) path character(len=1024) path
@ -338,9 +336,9 @@ end function
path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//newExt path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//newExt
if (present(recMultiplier)) then if (present(recMultiplier)) then
open(unit,status='replace',form='unformatted',access='direct',recl=pReal*recMultiplier,iostat=stat,file=path) open(myUnit,status='replace',form='unformatted',access='direct',recl=pReal*recMultiplier,iostat=stat,file=path)
else else
open(unit,status='replace',form='unformatted',access='direct',recl=pReal,iostat=stat,file=path) open(myUnit,status='replace',form='unformatted',access='direct',recl=pReal,iostat=stat,file=path)
endif endif
if (stat /= 0_pInt) call IO_error(100_pInt,ext_msg=path) if (stat /= 0_pInt) call IO_error(100_pInt,ext_msg=path)
@ -349,15 +347,15 @@ end function
!******************************************************************** !********************************************************************
! open (read) binary file related to restored job ! open (read) binary file related to restored job
! and with different extension to given unit ! and with different extension to given myUnit
!******************************************************************** !********************************************************************
subroutine IO_read_jobBinaryFile(unit,newExt,jobName,recMultiplier) subroutine IO_read_jobBinaryFile(myUnit,newExt,jobName,recMultiplier)
use prec, only: pReal, pInt use prec, only: pReal, pInt
use DAMASK_interface use DAMASK_interface
implicit none implicit none
integer(pInt), intent(in) :: unit integer(pInt), intent(in) :: myUnit
integer(pInt), intent(in), optional :: recMultiplier integer(pInt), intent(in), optional :: recMultiplier
character(len=*), intent(in) :: newExt, jobName character(len=*), intent(in) :: newExt, jobName
character(len=1024) path character(len=1024) path
@ -365,9 +363,9 @@ end function
path = trim(getSolverWorkingDirectoryName())//trim(jobName)//'.'//newExt path = trim(getSolverWorkingDirectoryName())//trim(jobName)//'.'//newExt
if (present(recMultiplier)) then if (present(recMultiplier)) then
open(unit,status='old',form='unformatted',access='direct',recl=pReal*recMultiplier,iostat=stat,file=path) open(myUnit,status='old',form='unformatted',access='direct',recl=pReal*recMultiplier,iostat=stat,file=path)
else else
open(unit,status='old',form='unformatted',access='direct',recl=pReal,iostat=stat,file=path) open(myUnit,status='old',form='unformatted',access='direct',recl=pReal,iostat=stat,file=path)
endif endif
if (stat /= 0) then if (stat /= 0) then
call IO_error(100_pInt,ext_msg=path) call IO_error(100_pInt,ext_msg=path)
@ -390,9 +388,9 @@ end function
real(pReal), intent(in) :: C real(pReal), intent(in) :: C
hybridIA_reps = 0_pInt hybridIA_reps = 0_pInt
do phi1=1,steps(1) do phi1=1_pInt,steps(1)
do Phi =1,steps(2) do Phi =1_pInt,steps(2)
do phi2=1,steps(3) do phi2=1_pInt,steps(3)
hybridIA_reps = hybridIA_reps+nint(C*dV_V(phi2,Phi,phi1), pInt) hybridIA_reps = hybridIA_reps+nint(C*dV_V(phi2,Phi,phi1), pInt)
enddo enddo
enddo enddo
@ -416,7 +414,7 @@ end function
character(len=80) line character(len=80) line
character(len=*), parameter :: fileFormat = '(A80)' character(len=*), parameter :: fileFormat = '(A80)'
integer(pInt) i,j,bin,Nast,NnonZero,Nset,Nreps,reps,phi1,Phi,phi2 integer(pInt) i,j,bin,Nast,NnonZero,Nset,Nreps,reps,phi1,Phi,phi2
integer(pInt), dimension(7) :: pos integer(pInt), dimension(7) :: myPos
integer(pInt), dimension(3) :: steps integer(pInt), dimension(3) :: steps
integer(pInt), dimension(:), allocatable :: binSet integer(pInt), dimension(:), allocatable :: binSet
real(pReal) center,sum_dV_V,prob,dg_0,C,lowerC,upperC,rnd real(pReal) center,sum_dV_V,prob,dg_0,C,lowerC,upperC,rnd
@ -429,18 +427,18 @@ end function
!--- parse header of ODF file --- !--- parse header of ODF file ---
!--- limits in phi1, Phi, phi2 --- !--- limits in phi1, Phi, phi2 ---
read(999,fmt=fileFormat,end=100) line read(999,fmt=fileFormat,end=100) line
pos = IO_stringPos(line,3_pInt) myPos = IO_stringPos(line,3_pInt)
if (pos(1).ne.3) goto 100 if (myPos(1).ne.3) goto 100
do i=1,3 do i=1_pInt,3_pInt
limits(i) = IO_floatValue(line,pos,i)*inRad limits(i) = IO_floatValue(line,myPos,i)*inRad
enddo enddo
!--- deltas in phi1, Phi, phi2 --- !--- deltas in phi1, Phi, phi2 ---
read(999,fmt=fileFormat,end=100) line read(999,fmt=fileFormat,end=100) line
pos = IO_stringPos(line,3_pInt) myPos = IO_stringPos(line,3_pInt)
if (pos(1).ne.3) goto 100 if (myPos(1).ne.3) goto 100
do i=1,3 do i=1_pInt,3_pInt
deltas(i) = IO_floatValue(line,pos,i)*inRad deltas(i) = IO_floatValue(line,myPos,i)*inRad
enddo enddo
steps = nint(limits/deltas,pInt) steps = nint(limits/deltas,pInt)
allocate(dV_V(steps(3),steps(2),steps(1))) allocate(dV_V(steps(3),steps(2),steps(1)))
@ -461,12 +459,12 @@ end function
dg_0 = deltas(1)*deltas(3)*2.0_pReal*sin(deltas(2)/2.0_pReal) dg_0 = deltas(1)*deltas(3)*2.0_pReal*sin(deltas(2)/2.0_pReal)
NnonZero = 0_pInt NnonZero = 0_pInt
do phi1=1,steps(1) do phi1=1_pInt,steps(1)
do Phi=1,steps(2) do Phi=1_pInt,steps(2)
do phi2=1,steps(3) do phi2=1_pInt,steps(3)
read(999,fmt=*,end=100) prob read(999,fmt=*,end=100) prob
if (prob > 0.0_pReal) then if (prob > 0.0_pReal) then
NnonZero = NnonZero+1 NnonZero = NnonZero+1_pInt
sum_dV_V = sum_dV_V+prob sum_dV_V = sum_dV_V+prob
else else
prob = 0.0_pReal prob = 0.0_pReal
@ -506,19 +504,19 @@ end function
allocate(binSet(Nreps)) allocate(binSet(Nreps))
bin = 0_pInt ! bin counter bin = 0_pInt ! bin counter
i = 1 ! set counter i = 1_pInt ! set counter
do phi1=1,steps(1) do phi1=1_pInt,steps(1)
do Phi=1,steps(2) do Phi=1_pInt,steps(2)
do phi2=1,steps(3) do phi2=1_pInt,steps(3)
reps = nint(C*dV_V(phi2,Phi,phi1), pInt) reps = nint(C*dV_V(phi2,Phi,phi1), pInt)
binSet(i:i+reps-1) = bin binSet(i:i+reps-1) = bin
bin = bin+1 ! advance bin bin = bin+1_pInt ! advance bin
i = i+reps ! advance set i = i+reps ! advance set
enddo enddo
enddo enddo
enddo enddo
do i=1,Nast do i=1_pInt,Nast
if (i < Nast) then if (i < Nast) then
call random_number(rnd) call random_number(rnd)
j = nint(rnd*(Nreps-i)+i+0.5_pReal,pInt) j = nint(rnd*(Nreps-i)+i+0.5_pReal,pInt)
@ -552,7 +550,7 @@ end function
character(len=*), intent(in) :: line character(len=*), intent(in) :: line
character(len=*), parameter :: blank = achar(32)//achar(9)//achar(10)//achar(13) ! whitespaces character(len=*), parameter :: blank = achar(32)//achar(9)//achar(10)//achar(13) ! whitespaces
character(len=*), parameter :: comment = achar(35) ! comment id '#' character(len=*), parameter :: comment = achar(35) ! comment id '#'
integer(pInt) posNonBlank, posComment integer :: posNonBlank, posComment !no pInt
logical IO_isBlank logical IO_isBlank
posNonBlank = verify(line,blank) posNonBlank = verify(line,blank)
@ -572,7 +570,7 @@ end function
character(len=*), intent(in) :: line,openChar,closeChar character(len=*), intent(in) :: line,openChar,closeChar
character(len=*), parameter :: sep=achar(32)//achar(9)//achar(10)//achar(13) ! whitespaces character(len=*), parameter :: sep=achar(32)//achar(9)//achar(10)//achar(13) ! whitespaces
character(len=len_trim(line)) IO_getTag character(len=len_trim(line)) IO_getTag
integer(pInt) left,right integer :: left,right !no pInt
IO_getTag = '' IO_getTag = ''
left = scan(line,openChar) left = scan(line,openChar)
@ -596,7 +594,7 @@ end function
integer(pInt) IO_countSections integer(pInt) IO_countSections
character(len=1024) line character(len=1024) line
IO_countSections = 0 IO_countSections = 0_pInt
line = '' line = ''
rewind(file) rewind(file)
@ -609,7 +607,7 @@ end function
if (IO_isBlank(line)) cycle ! skip empty lines if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') & ! found [section] identifier if (IO_getTag(line,'[',']') /= '') & ! found [section] identifier
IO_countSections = IO_countSections + 1 IO_countSections = IO_countSections + 1_pInt
enddo enddo
100 endfunction 100 endfunction
@ -627,7 +625,7 @@ end function
integer(pInt), intent(in) :: file, Nsections integer(pInt), intent(in) :: file, Nsections
character(len=*), intent(in) :: part, myTag character(len=*), intent(in) :: part, myTag
integer(pInt), dimension(Nsections) :: IO_countTagInPart, counter integer(pInt), dimension(Nsections) :: IO_countTagInPart, counter
integer(pInt), parameter :: maxNchunks = 1 integer(pInt), parameter :: maxNchunks = 1_pInt
integer(pInt), dimension(1+2*maxNchunks) :: positions integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt) section integer(pInt) section
character(len=1024) line,tag character(len=1024) line,tag
@ -646,7 +644,7 @@ end function
if (IO_isBlank(line)) cycle ! skip empty lines if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') & ! found [section] identifier if (IO_getTag(line,'[',']') /= '') & ! found [section] identifier
section = section + 1 section = section + 1_pInt
if (section > 0) then if (section > 0) then
positions = IO_stringPos(line,maxNchunks) positions = IO_stringPos(line,maxNchunks)
tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key
@ -672,7 +670,7 @@ endfunction
integer(pInt), intent(in) :: file, Nsections integer(pInt), intent(in) :: file, Nsections
character(len=*), intent(in) :: part, myTag character(len=*), intent(in) :: part, myTag
logical, dimension(Nsections) :: IO_spotTagInPart logical, dimension(Nsections) :: IO_spotTagInPart
integer(pInt), parameter :: maxNchunks = 1 integer(pInt), parameter :: maxNchunks = 1_pInt
integer(pInt), dimension(1+2*maxNchunks) :: positions integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt) section integer(pInt) section
character(len=1024) line,tag character(len=1024) line,tag
@ -691,8 +689,8 @@ endfunction
if (IO_isBlank(line)) cycle ! skip empty lines if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') & ! found [section] identifier if (IO_getTag(line,'[',']') /= '') & ! found [section] identifier
section = section + 1 section = section + 1_pInt
if (section > 0) then if (section > 0_pInt) then
positions = IO_stringPos(line,maxNchunks) positions = IO_stringPos(line,maxNchunks)
tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key
if (tag == myTag) & ! match if (tag == myTag) & ! match
@ -716,11 +714,11 @@ endfunction
character(len=*), intent(in) :: line character(len=*), intent(in) :: line
character(len=*), parameter :: sep=achar(44)//achar(32)//achar(9)//achar(10)//achar(13) ! comma and whitespaces character(len=*), parameter :: sep=achar(44)//achar(32)//achar(9)//achar(10)//achar(13) ! comma and whitespaces
integer(pInt), intent(in) :: N integer(pInt), intent(in) :: N
integer(pInt) left,right integer :: left,right !no pInt (verify and scan return default integer)
integer(pInt) IO_stringPos(1+N*2) integer(pInt) :: IO_stringPos(1_pInt+N*2_pInt)
IO_stringPos = -1 IO_stringPos = -1_pInt
IO_stringPos(1) = 0 IO_stringPos(1) = 0_pInt
right = 0 right = 0
do while (verify(line(right+1:),sep)>0) do while (verify(line(right+1:),sep)>0)
@ -730,55 +728,55 @@ endfunction
exit exit
endif endif
if ( IO_stringPos(1)<N ) then if ( IO_stringPos(1)<N ) then
IO_stringPos(1+IO_stringPos(1)*2+1) = left IO_stringPos(1_pInt+IO_stringPos(1)*2_pInt+1_pInt) = int(left, pInt)
IO_stringPos(1+IO_stringPos(1)*2+2) = right IO_stringPos(1_pInt+IO_stringPos(1)*2_pInt+2_pInt) = int(right, pInt)
endif endif
IO_stringPos(1) = IO_stringPos(1)+1 IO_stringPos(1) = IO_stringPos(1)+1_pInt
enddo enddo
endfunction endfunction
!******************************************************************** !********************************************************************
! read string value at pos from line ! read string value at myPos from line
!******************************************************************** !********************************************************************
pure function IO_stringValue (line,positions,pos) pure function IO_stringValue (line,positions,myPos)
use prec, only: pReal,pInt use prec, only: pReal,pInt
implicit none implicit none
character(len=*), intent(in) :: line character(len=*), intent(in) :: line
integer(pInt), intent(in) :: positions(*),pos integer(pInt), intent(in) :: positions(*),myPos
character(len=1+positions(pos*2+1)-positions(pos*2)) IO_stringValue character(len=1+positions(myPos*2+1)-positions(myPos*2)) :: IO_stringValue
if (positions(1) < pos) then if (positions(1) < myPos) then
IO_stringValue = '' IO_stringValue = ''
else else
IO_stringValue = line(positions(pos*2):positions(pos*2+1)) IO_stringValue = line(positions(myPos*2):positions(myPos*2+1))
endif endif
endfunction endfunction
!******************************************************************** !********************************************************************
! read string value at pos from fixed format line ! read string value at myPos from fixed format line
!******************************************************************** !********************************************************************
pure function IO_fixedStringValue (line,ends,pos) pure function IO_fixedStringValue (line,ends,myPos)
use prec, only: pReal,pInt use prec, only: pReal,pInt
implicit none implicit none
character(len=*), intent(in) :: line character(len=*), intent(in) :: line
integer(pInt), intent(in) :: ends(*),pos integer(pInt), intent(in) :: ends(*),myPos
character(len=ends(pos+1)-ends(pos)) IO_fixedStringValue character(len=ends(myPos+1)-ends(myPos)) :: IO_fixedStringValue
IO_fixedStringValue = line(ends(pos)+1:ends(pos+1)) IO_fixedStringValue = line(ends(myPos)+1:ends(myPos+1))
endfunction endfunction
!******************************************************************** !********************************************************************
! read float value at pos from line ! read float value at myPos from line
!******************************************************************** !********************************************************************
pure function IO_floatValue (line,positions,myPos) pure function IO_floatValue (line,positions,myPos)
@ -787,7 +785,7 @@ endfunction
character(len=*), intent(in) :: line character(len=*), intent(in) :: line
integer(pInt), intent(in) :: positions(*),myPos integer(pInt), intent(in) :: positions(*),myPos
real(pReal) IO_floatValue real(pReal) :: IO_floatValue
if (positions(1) < myPos) then if (positions(1) < myPos) then
IO_floatValue = 0.0_pReal IO_floatValue = 0.0_pReal
@ -801,18 +799,18 @@ endfunction
!******************************************************************** !********************************************************************
! read float value at pos from fixed format line ! read float value at myPos from fixed format line
!******************************************************************** !********************************************************************
pure function IO_fixedFloatValue (line,ends,pos) pure function IO_fixedFloatValue (line,ends,myPos)
use prec, only: pReal,pInt use prec, only: pReal,pInt
implicit none implicit none
character(len=*), intent(in) :: line character(len=*), intent(in) :: line
integer(pInt), intent(in) :: ends(*),pos integer(pInt), intent(in) :: ends(*),myPos
real(pReal) IO_fixedFloatValue real(pReal) :: IO_fixedFloatValue
read(UNIT=line(ends(pos-1)+1:ends(pos)),ERR=100,FMT=*) IO_fixedFloatValue read(UNIT=line(ends(myPos-1)+1:ends(myPos)),ERR=100,FMT=*) IO_fixedFloatValue
return return
100 IO_fixedFloatValue = huge(1.0_pReal) 100 IO_fixedFloatValue = huge(1.0_pReal)
@ -820,24 +818,25 @@ endfunction
!******************************************************************** !********************************************************************
! read float x.y+z value at pos from format line line ! read float x.y+z value at myPos from format line line
!******************************************************************** !********************************************************************
pure function IO_fixedNoEFloatValue (line,ends,pos) pure function IO_fixedNoEFloatValue (line,ends,myPos)
use prec, only: pReal,pInt use prec, only: pReal,pInt
implicit none implicit none
character(len=*), intent(in) :: line character(len=*), intent(in) :: line
integer(pInt), intent(in) :: ends(*),pos integer(pInt), intent(in) :: ends(*),myPos
integer(pInt) pos_exp,expon integer(pInt) :: expon
integer :: pos_exp
real(pReal) IO_fixedNoEFloatValue,base real(pReal) IO_fixedNoEFloatValue,base
pos_exp = scan(line(ends(pos)+1:ends(pos+1)),'+-',back=.true.) pos_exp = scan(line(ends(myPos)+1:ends(myPos+1)),'+-',back=.true.)
if (pos_exp > 1) then if (pos_exp > 1) then
read(UNIT=line(ends(pos)+1:ends(pos)+pos_exp-1),ERR=100,FMT=*) base read(UNIT=line(ends(myPos)+1:ends(myPos)+pos_exp-1),ERR=100,FMT=*) base
read(UNIT=line(ends(pos)+pos_exp:ends(pos+1)),ERR=100,FMT=*) expon read(UNIT=line(ends(myPos)+pos_exp:ends(myPos+1)),ERR=100,FMT=*) expon
else else
read(UNIT=line(ends(pos)+1:ends(pos+1)),ERR=100,FMT=*) base read(UNIT=line(ends(myPos)+1:ends(myPos+1)),ERR=100,FMT=*) base
expon = 0_pInt expon = 0_pInt
endif endif
IO_fixedNoEFloatValue = base*10.0_pReal**expon IO_fixedNoEFloatValue = base*10.0_pReal**expon
@ -848,21 +847,21 @@ endfunction
!******************************************************************** !********************************************************************
! read int value at pos from line ! read int value at myPos from line
!******************************************************************** !********************************************************************
pure function IO_intValue (line,positions,pos) pure function IO_intValue (line,positions,myPos)
use prec, only: pReal,pInt use prec, only: pReal,pInt
implicit none implicit none
character(len=*), intent(in) :: line character(len=*), intent(in) :: line
integer(pInt), intent(in) :: positions(*),pos integer(pInt), intent(in) :: positions(*),myPos
integer(pInt) IO_intValue integer(pInt) :: IO_intValue
if (positions(1) < pos) then if (positions(1) < myPos) then
IO_intValue = 0_pInt IO_intValue = 0_pInt
else else
read(UNIT=line(positions(pos*2):positions(pos*2+1)),ERR=100,FMT=*) IO_intValue read(UNIT=line(positions(myPos*2):positions(myPos*2+1)),ERR=100,FMT=*) IO_intValue
endif endif
return return
100 IO_intValue = huge(1_pInt) 100 IO_intValue = huge(1_pInt)
@ -871,18 +870,18 @@ endfunction
!******************************************************************** !********************************************************************
! read int value at pos from fixed format line ! read int value at myPos from fixed format line
!******************************************************************** !********************************************************************
pure function IO_fixedIntValue (line,ends,pos) pure function IO_fixedIntValue (line,ends,myPos)
use prec, only: pReal,pInt use prec, only: pReal,pInt
implicit none implicit none
character(len=*), intent(in) :: line character(len=*), intent(in) :: line
integer(pInt), intent(in) :: ends(*),pos integer(pInt), intent(in) :: ends(*),myPos
integer(pInt) IO_fixedIntValue integer(pInt) IO_fixedIntValue
read(UNIT=line(ends(pos)+1:ends(pos+1)),ERR=100,FMT=*) IO_fixedIntValue read(UNIT=line(ends(myPos)+1:ends(myPos+1)),ERR=100,FMT=*) IO_fixedIntValue
return return
100 IO_fixedIntValue = huge(1_pInt) 100 IO_fixedIntValue = huge(1_pInt)
@ -899,7 +898,7 @@ endfunction
character (len=*), intent(in) :: line character (len=*), intent(in) :: line
character (len=len(line)) IO_lc character (len=len(line)) IO_lc
integer(pInt) i integer :: i !no pInt (len returns default integer)
IO_lc = line IO_lc = line
do i=1,len(line) do i=1,len(line)
@ -919,7 +918,7 @@ endfunction
character (len=*) line character (len=*) line
character (len=len(line)) IO_lc character (len=len(line)) IO_lc
integer(pInt) i integer :: i !no pInt (len returns default integer)
IO_lc = line IO_lc = line
do i=1,len(line) do i=1,len(line)
@ -939,15 +938,15 @@ endfunction
implicit none implicit none
integer(pInt) remainingChunks,unit,N integer(pInt) remainingChunks,unit,N
integer(pInt), parameter :: maxNchunks = 64 integer(pInt), parameter :: maxNchunks = 64_pInt
integer(pInt), dimension(1+2*maxNchunks) :: pos integer(pInt), dimension(1+2*maxNchunks) :: myPos
character(len=300) line character(len=300) line
remainingChunks = N remainingChunks = N
do while (remainingChunks > 0) do while (remainingChunks > 0)
read(unit,'(A300)',end=100) line read(unit,'(A300)',end=100) line
pos = IO_stringPos(line,maxNchunks) myPos = IO_stringPos(line,maxNchunks)
remainingChunks = remainingChunks - pos(1) remainingChunks = remainingChunks - myPos(1)
enddo enddo
100 endsubroutine 100 endsubroutine
@ -957,19 +956,18 @@ endfunction
!******************************************************************** !********************************************************************
pure function IO_extractValue (line,key) pure function IO_extractValue (line,key)
use prec, only: pReal,pInt
implicit none implicit none
character(len=*), intent(in) :: line,key character(len=*), intent(in) :: line,key
character(len=*), parameter :: sep = achar(61) ! '=' character(len=*), parameter :: sep = achar(61) ! '='
integer(pInt) pos integer :: myPos ! no pInt (scan returns default integer)
character(len=300) IO_extractValue character(len=300) IO_extractValue
IO_extractValue = '' IO_extractValue = ''
pos = scan(line,sep) myPos = scan(line,sep)
if (pos > 0 .and. line(:pos-1) == key(:pos-1)) & ! key matches expected key if (myPos > 0 .and. line(:myPos-1) == key(:myPos-1)) & ! key matches expected key
IO_extractValue = line(pos+1:) ! extract value IO_extractValue = line(myPos+1:) ! extract value
endfunction endfunction
@ -980,28 +978,28 @@ endfunction
! : is not changed back to the original version since *.inp_assembly does not ! : is not changed back to the original version since *.inp_assembly does not
! : contain any comment lines (12.07.2010) ! : contain any comment lines (12.07.2010)
!******************************************************************** !********************************************************************
function IO_countDataLines (unit) function IO_countDataLines (myUnit)
use prec, only: pReal,pInt use prec, only: pReal,pInt
implicit none implicit none
integer(pInt) IO_countDataLines,unit integer(pInt) :: IO_countDataLines,myUnit
integer(pInt), parameter :: maxNchunks = 1 integer(pInt), parameter :: maxNchunks = 1_pInt
integer(pInt), dimension(1+2*maxNchunks) :: pos integer(pInt), dimension(1+2*maxNchunks) :: myPos
character(len=300) line,tmp character(len=300) :: line,tmp
IO_countDataLines = 0 IO_countDataLines = 0_pInt
do do
read(unit,'(A300)',end=100) line read(myUnit,'(A300)',end=100) line
pos = IO_stringPos(line,maxNchunks) myPos = IO_stringPos(line,maxNchunks)
tmp = IO_lc(IO_stringValue(line,pos,1_pInt)) tmp = IO_lc(IO_stringValue(line,myPos,1_pInt))
if (tmp(1:1) == '*' .and. tmp(2:2) /= '*') then ! found keyword if (tmp(1:1) == '*' .and. tmp(2:2) /= '*') then ! found keyword
exit exit
else else
if (tmp(2:2) /= '*') IO_countDataLines = IO_countDataLines + 1_pInt if (tmp(2:2) /= '*') IO_countDataLines = IO_countDataLines + 1_pInt
endif endif
enddo enddo
100 backspace(unit) 100 backspace(myUnit)
endfunction endfunction
@ -1011,16 +1009,16 @@ endfunction
! Marc: ints concatenated by "c" as last char or range of values a "to" b ! Marc: ints concatenated by "c" as last char or range of values a "to" b
! Abaqus: triplet of start,stop,inc ! Abaqus: triplet of start,stop,inc
!******************************************************************** !********************************************************************
function IO_countContinousIntValues (unit) function IO_countContinousIntValues (myUnit)
use DAMASK_interface use DAMASK_interface
use prec, only: pReal,pInt use prec, only: pReal,pInt
implicit none implicit none
integer(pInt) unit,l,count integer(pInt) :: myUnit,l,c
integer(pInt) IO_countContinousIntValues integer(pInt) :: IO_countContinousIntValues
integer(pInt), parameter :: maxNchunks = 8192 integer(pInt), parameter :: maxNchunks = 8192_pInt
integer(pInt), dimension(1+2*maxNchunks) :: pos integer(pInt), dimension(1+2*maxNchunks) :: myPos
character(len=65536) line character(len=65536) line
IO_countContinousIntValues = 0_pInt IO_countContinousIntValues = 0_pInt
@ -1029,15 +1027,15 @@ endfunction
case ('Marc','Spectral') case ('Marc','Spectral')
do do
read(unit,'(A300)',end=100) line read(myUnit,'(A300)',end=100) line
pos = IO_stringPos(line,maxNchunks) myPos = IO_stringPos(line,maxNchunks)
if (IO_lc(IO_stringValue(line,pos,2_pInt)) == 'to' ) then ! found range indicator if (IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'to' ) then ! found range indicator
IO_countContinousIntValues = 1_pInt + IO_intValue(line,pos,3_pInt) - IO_intValue(line,pos,1_pInt) IO_countContinousIntValues = 1_pInt + IO_intValue(line,myPos,3_pInt) - IO_intValue(line,myPos,1_pInt)
exit ! only one single range indicator allowed exit ! only one single range indicator allowed
else else
IO_countContinousIntValues = IO_countContinousIntValues+pos(1)-1_pInt ! add line's count when assuming 'c' IO_countContinousIntValues = IO_countContinousIntValues+myPos(1)-1_pInt ! add line's count when assuming 'c'
if ( IO_lc(IO_stringValue(line,pos,pos(1))) /= 'c' ) then ! line finished, read last value if ( IO_lc(IO_stringValue(line,myPos,myPos(1))) /= 'c' ) then ! line finished, read last value
IO_countContinousIntValues = IO_countContinousIntValues+1 IO_countContinousIntValues = IO_countContinousIntValues+1_pInt
exit ! data ended exit ! data ended
endif endif
endif endif
@ -1045,17 +1043,17 @@ endfunction
case('Abaqus') case('Abaqus')
count = IO_countDataLines(unit) c = IO_countDataLines(myUnit)
do l = 1,count do l = 1_pInt,c
backspace(unit) backspace(myUnit)
enddo enddo
do l = 1,count do l = 1_pInt,c
read(unit,'(A300)',end=100) line read(myUnit,'(A300)',end=100) line
pos = IO_stringPos(line,maxNchunks) myPos = IO_stringPos(line,maxNchunks)
IO_countContinousIntValues = IO_countContinousIntValues + 1 + & ! assuming range generation IO_countContinousIntValues = IO_countContinousIntValues + 1_pInt + & ! assuming range generation
(IO_intValue(line,pos,2_pInt)-IO_intValue(line,pos,1_pInt))/& (IO_intValue(line,myPos,2_pInt)-IO_intValue(line,myPos,1_pInt))/&
max(1_pInt,IO_intValue(line,pos,3_pInt)) max(1_pInt,IO_intValue(line,myPos,3_pInt))
enddo enddo
endselect endselect
@ -1068,53 +1066,53 @@ endfunction
! Marc: ints concatenated by "c" as last char, range of a "to" b, or named set ! Marc: ints concatenated by "c" as last char, range of a "to" b, or named set
! Abaqus: triplet of start,stop,inc or named set ! Abaqus: triplet of start,stop,inc or named set
!******************************************************************** !********************************************************************
function IO_continousIntValues (unit,maxN,lookupName,lookupMap,lookupMaxN) function IO_continousIntValues (myUnit,maxN,lookupName,lookupMap,lookupMaxN)
use DAMASK_interface use DAMASK_interface
use prec, only: pReal,pInt use prec, only: pReal,pInt
implicit none implicit none
integer(pInt) unit,maxN,i,j,l,count,first,last integer(pInt) myUnit,maxN,i,j,l,c,first,last
integer(pInt), dimension(1+maxN) :: IO_continousIntValues integer(pInt), dimension(1+maxN) :: IO_continousIntValues
integer(pInt), parameter :: maxNchunks = 8192_pInt integer(pInt), parameter :: maxNchunks = 8192_pInt
integer(pInt), dimension(1+2*maxNchunks) :: pos integer(pInt), dimension(1+2*maxNchunks) :: myPos
character(len=64), dimension(:) :: lookupName character(len=64), dimension(:) :: lookupName
integer(pInt) :: lookupMaxN integer(pInt) :: lookupMaxN
integer(pInt), dimension(:,:) :: lookupMap integer(pInt), dimension(:,:) :: lookupMap
character(len=65536) line character(len=65536) line
logical rangeGeneration logical rangeGeneration
IO_continousIntValues = 0 IO_continousIntValues = 0_pInt
rangeGeneration = .false. rangeGeneration = .false.
select case (FEsolver) select case (FEsolver)
case ('Marc','Spectral') case ('Marc','Spectral')
do do
read(unit,'(A65536)',end=100) line read(myUnit,'(A65536)',end=100) line
pos = IO_stringPos(line,maxNchunks) myPos = IO_stringPos(line,maxNchunks)
if (verify(IO_stringValue(line,pos,1_pInt),'0123456789') > 0) then ! a non-int, i.e. set name if (verify(IO_stringValue(line,myPos,1_pInt),'0123456789') > 0) then ! a non-int, i.e. set name
do i = 1_pInt,lookupMaxN ! loop over known set names do i = 1_pInt, lookupMaxN ! loop over known set names
if (IO_stringValue(line,pos,1_pInt) == lookupName(i)) then ! found matching name if (IO_stringValue(line,myPos,1_pInt) == lookupName(i)) then ! found matching name
IO_continousIntValues = lookupMap(:,i) ! return resp. entity list IO_continousIntValues = lookupMap(:,i) ! return resp. entity list
exit exit
endif endif
enddo enddo
exit exit
else if (pos(1) > 2_pInt .and. IO_lc(IO_stringValue(line,pos,2_pInt)) == 'to' ) then ! found range indicator else if (myPos(1) > 2_pInt .and. IO_lc(IO_stringValue(line,myPos,2_pInt)) == 'to' ) then ! found range indicator
do i = IO_intValue(line,pos,1_pInt),IO_intValue(line,pos,3_pInt) do i = IO_intValue(line,myPos,1_pInt),IO_intValue(line,myPos,3_pInt)
IO_continousIntValues(1) = IO_continousIntValues(1) + 1_pInt IO_continousIntValues(1) = IO_continousIntValues(1) + 1_pInt
IO_continousIntValues(1+IO_continousIntValues(1)) = i IO_continousIntValues(1+IO_continousIntValues(1)) = i
enddo enddo
exit exit
else else
do i = 1,pos(1)-1 ! interpret up to second to last value do i = 1_pInt ,myPos(1)-1_pInt ! interpret up to second to last value
IO_continousIntValues(1) = IO_continousIntValues(1) + 1_pInt IO_continousIntValues(1) = IO_continousIntValues(1) + 1_pInt
IO_continousIntValues(1+IO_continousIntValues(1)) = IO_intValue(line,pos,i) IO_continousIntValues(1+IO_continousIntValues(1)) = IO_intValue(line,myPos,i)
enddo enddo
if ( IO_lc(IO_stringValue(line,pos,pos(1))) /= 'c' ) then ! line finished, read last value if ( IO_lc(IO_stringValue(line,myPos,myPos(1))) /= 'c' ) then ! line finished, read last value
IO_continousIntValues(1) = IO_continousIntValues(1) + 1_pInt IO_continousIntValues(1) = IO_continousIntValues(1) + 1_pInt
IO_continousIntValues(1+IO_continousIntValues(1)) = IO_intValue(line,pos,pos(1)) IO_continousIntValues(1+IO_continousIntValues(1)) = IO_intValue(line,myPos,myPos(1))
exit exit
endif endif
endif endif
@ -1122,43 +1120,43 @@ endfunction
case('Abaqus') case('Abaqus')
count = IO_countDataLines(unit) c = IO_countDataLines(myUnit)
do l = 1,count do l = 1_pInt,c
backspace(unit) backspace(myUnit)
enddo enddo
! check if the element values in the elset are auto generated ! check if the element values in the elset are auto generated
backspace(unit) backspace(myUnit)
read(unit,'(A65536)',end=100) line read(myUnit,'(A65536)',end=100) line
pos = IO_stringPos(line,maxNchunks) myPos = IO_stringPos(line,maxNchunks)
do i = 1,pos(1) do i = 1_pInt,myPos(1)
if (IO_lc(IO_stringValue(line,pos,i)) == 'generate') rangeGeneration = .true. if (IO_lc(IO_stringValue(line,myPos,i)) == 'generate') rangeGeneration = .true.
enddo enddo
do l = 1,count do l = 1_pInt,c
read(unit,'(A65536)',end=100) line read(myUnit,'(A65536)',end=100) line
pos = IO_stringPos(line,maxNchunks) myPos = IO_stringPos(line,maxNchunks)
if (verify(IO_stringValue(line,pos,1_pInt),'0123456789') > 0_pInt) then ! a non-int, i.e. set names follow on this line if (verify(IO_stringValue(line,myPos,1_pInt),'0123456789') > 0) then ! a non-int, i.e. set names follow on this line
do i = 1,pos(1) ! loop over set names in line do i = 1_pInt,myPos(1) ! loop over set names in line
do j = 1,lookupMaxN ! look thru known set names do j = 1_pInt,lookupMaxN ! look thru known set names
if (IO_stringValue(line,pos,i) == lookupName(j)) then ! found matching name if (IO_stringValue(line,myPos,i) == lookupName(j)) then ! found matching name
first = 2 + IO_continousIntValues(1) ! where to start appending data first = 2_pInt + IO_continousIntValues(1) ! where to start appending data
last = first + lookupMap(1,j) - 1 ! up to where to append data last = first + lookupMap(1,j) - 1_pInt ! up to where to append data
IO_continousIntValues(first:last) = lookupMap(2:1+lookupMap(1,j),j) ! add resp. entity list IO_continousIntValues(first:last) = lookupMap(2:1+lookupMap(1,j),j) ! add resp. entity list
IO_continousIntValues(1) = IO_continousIntValues(1) + lookupMap(1,j) ! count them IO_continousIntValues(1) = IO_continousIntValues(1) + lookupMap(1,j) ! count them
endif endif
enddo enddo
enddo enddo
else if (rangeGeneration) then ! range generation else if (rangeGeneration) then ! range generation
do i = IO_intValue(line,pos,1_pInt),IO_intValue(line,pos,2_pInt),max(1_pInt,IO_intValue(line,pos,3_pInt)) do i = IO_intValue(line,myPos,1_pInt),IO_intValue(line,myPos,2_pInt),max(1_pInt,IO_intValue(line,myPos,3_pInt))
IO_continousIntValues(1) = IO_continousIntValues(1) + 1 IO_continousIntValues(1) = IO_continousIntValues(1) + 1_pInt
IO_continousIntValues(1+IO_continousIntValues(1)) = i IO_continousIntValues(1+IO_continousIntValues(1)) = i
enddo enddo
else ! read individual elem nums else ! read individual elem nums
do i = 1,pos(1) do i = 1_pInt,myPos(1)
! write(*,*)'IO_CIV-int',IO_intValue(line,pos,i) ! write(*,*)'IO_CIV-int',IO_intValue(line,myPos,i)
IO_continousIntValues(1) = IO_continousIntValues(1) + 1 IO_continousIntValues(1) = IO_continousIntValues(1) + 1_pInt
IO_continousIntValues(1+IO_continousIntValues(1)) = IO_intValue(line,pos,i) IO_continousIntValues(1+IO_continousIntValues(1)) = IO_intValue(line,myPos,i)
enddo enddo
endif endif
enddo enddo

24
code/crystallite.f90
View File

@ -3178,7 +3178,7 @@ function crystallite_postResults(&
crystallite_postResults = 0.0_pReal crystallite_postResults = 0.0_pReal
c = 0_pInt c = 0_pInt
crystallite_postResults(c+1) = crystallite_sizePostResults(crystID) ! size of results from cryst crystallite_postResults(c+1) = real(crystallite_sizePostResults(crystID),pReal) ! size of results from cryst
c = c + 1_pInt c = c + 1_pInt
do o = 1,crystallite_Noutput(crystID) do o = 1,crystallite_Noutput(crystID)
@ -3186,10 +3186,10 @@ function crystallite_postResults(&
select case(crystallite_output(o,crystID)) select case(crystallite_output(o,crystID))
case ('phase') case ('phase')
mySize = 1_pInt mySize = 1_pInt
crystallite_postResults(c+1) = material_phase(g,i,e) ! phaseID of grain crystallite_postResults(c+1) = real(material_phase(g,i,e),pReal) ! phaseID of grain
case ('texture') case ('texture')
mySize = 1_pInt mySize = 1_pInt
crystallite_postResults(c+1) = material_texture(g,i,e) ! textureID of grain crystallite_postResults(c+1) = real(material_texture(g,i,e),pReal) ! textureID of grain
case ('volume') case ('volume')
mySize = 1_pInt mySize = 1_pInt
detF = math_det33(crystallite_partionedF(1:3,1:3,g,i,e)) ! V_current = det(F) * V_reference detF = math_det33(crystallite_partionedF(1:3,1:3,g,i,e)) ! V_current = det(F) * V_reference
@ -3210,36 +3210,36 @@ function crystallite_postResults(&
case ('defgrad','f') case ('defgrad','f')
mySize = 9_pInt mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_partionedF(1:3,1:3,g,i,e)),(/mySize/)) crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_partionedF(1:3,1:3,g,i,e)),[mySize])
case ('e') case ('e')
mySize = 9_pInt mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = 0.5_pReal * reshape((math_mul33x33( & crystallite_postResults(c+1:c+mySize) = 0.5_pReal * reshape((math_mul33x33( &
math_transpose33(crystallite_partionedF(1:3,1:3,g,i,e)), & math_transpose33(crystallite_partionedF(1:3,1:3,g,i,e)), &
crystallite_partionedF(1:3,1:3,g,i,e)) - math_I3),(/mySize/)) crystallite_partionedF(1:3,1:3,g,i,e)) - math_I3),[mySize])
case ('fe') case ('fe')
mySize = 9_pInt mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_Fe(1:3,1:3,g,i,e)),(/mySize/)) crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_Fe(1:3,1:3,g,i,e)),[mySize])
case ('ee') case ('ee')
Ee = 0.5_pReal * (math_mul33x33(math_transpose33(crystallite_Fe(1:3,1:3,g,i,e)), crystallite_Fe(1:3,1:3,g,i,e)) - math_I3) Ee = 0.5_pReal * (math_mul33x33(math_transpose33(crystallite_Fe(1:3,1:3,g,i,e)), crystallite_Fe(1:3,1:3,g,i,e)) - math_I3)
mySize = 9_pInt mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = reshape(Ee,(/mySize/)) crystallite_postResults(c+1:c+mySize) = reshape(Ee,[mySize])
case ('fp') case ('fp')
mySize = 9_pInt mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_Fp(1:3,1:3,g,i,e)),(/mySize/)) crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_Fp(1:3,1:3,g,i,e)),[mySize])
case ('lp') case ('lp')
mySize = 9_pInt mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_Lp(1:3,1:3,g,i,e)),(/mySize/)) crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_Lp(1:3,1:3,g,i,e)),[mySize])
case ('p','firstpiola','1stpiola') case ('p','firstpiola','1stpiola')
mySize = 9_pInt mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_P(1:3,1:3,g,i,e)),(/mySize/)) crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_P(1:3,1:3,g,i,e)),[mySize])
case ('s','tstar','secondpiola','2ndpiola') case ('s','tstar','secondpiola','2ndpiola')
mySize = 9_pInt mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = reshape(math_Mandel6to33(crystallite_Tstar_v(1:6,g,i,e)),(/mySize/)) crystallite_postResults(c+1:c+mySize) = reshape(math_Mandel6to33(crystallite_Tstar_v(1:6,g,i,e)),[mySize])
end select end select
c = c + mySize c = c + mySize
enddo enddo
crystallite_postResults(c+1) = constitutive_sizePostResults(g,i,e) ! size of constitutive results crystallite_postResults(c+1) = real(constitutive_sizePostResults(g,i,e),pReal) ! size of constitutive results
c = c + 1_pInt c = c + 1_pInt
crystallite_postResults(c+1:c+constitutive_sizePostResults(g,i,e)) = constitutive_postResults(crystallite_Tstar_v(1:6,g,i,e), & crystallite_postResults(c+1:c+constitutive_sizePostResults(g,i,e)) = constitutive_postResults(crystallite_Tstar_v(1:6,g,i,e), &
crystallite_Fe, & crystallite_Fe, &

4
code/homogenization.f90
View File

@ -562,7 +562,7 @@ subroutine materialpoint_postResults(dt)
thePos = 0_pInt thePos = 0_pInt
theSize = homogenization_sizePostResults(i,e) theSize = homogenization_sizePostResults(i,e)
materialpoint_results(thePos+1,i,e) = theSize ! tell size of homogenization results materialpoint_results(thePos+1,i,e) = real(theSize,pReal) ! tell size of homogenization results
thePos = thePos + 1_pInt thePos = thePos + 1_pInt
if (theSize > 0_pInt) then ! any homogenization results to mention? if (theSize > 0_pInt) then ! any homogenization results to mention?
@ -570,7 +570,7 @@ subroutine materialpoint_postResults(dt)
thePos = thePos + theSize thePos = thePos + theSize
endif endif
materialpoint_results(thePos+1,i,e) = myNgrains ! tell number of grains at materialpoint materialpoint_results(thePos+1,i,e) = real(myNgrains,pReal) ! tell number of grains at materialpoint
thePos = thePos + 1_pInt thePos = thePos + 1_pInt
do g = 1,myNgrains ! loop over all grains do g = 1,myNgrains ! loop over all grains

58
code/homogenization_RGC.f90
View File

@ -281,7 +281,7 @@ subroutine homogenization_RGC_partitionDeformation(&
write(6,'(1x,a,i3,a,i3,a)')'========== Increment: ',theInc,' Cycle: ',cycleCounter,' ==========' write(6,'(1x,a,i3,a,i3,a)')'========== Increment: ',theInc,' Cycle: ',cycleCounter,' =========='
write(6,'(1x,a32)')'Overall deformation gradient: ' write(6,'(1x,a32)')'Overall deformation gradient: '
do i = 1_pInt,3_pInt do i = 1_pInt,3_pInt
write(6,'(1x,3(e14.8,1x))')(avgF(i,j), j = 1_pInt,3_pInt) write(6,'(1x,3(e15.8,1x))')(avgF(i,j), j = 1_pInt,3_pInt)
enddo enddo
write(6,*)' ' write(6,*)' '
call flush(6) call flush(6)
@ -307,7 +307,7 @@ subroutine homogenization_RGC_partitionDeformation(&
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(1x,a32,1x,i3)')'Deformation gradient of grain: ',iGrain write(6,'(1x,a32,1x,i3)')'Deformation gradient of grain: ',iGrain
do i = 1_pInt,3_pInt do i = 1_pInt,3_pInt
write(6,'(1x,3(e14.8,1x))')(F(i,j,iGrain), j = 1_pInt,3_pInt) write(6,'(1x,3(e15.8,1x))')(F(i,j,iGrain), j = 1_pInt,3_pInt)
enddo enddo
write(6,*)' ' write(6,*)' '
call flush(6) call flush(6)
@ -394,7 +394,7 @@ function homogenization_RGC_updateState(&
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Obtained state: ' write(6,'(1x,a30)')'Obtained state: '
do i = 1,3*nIntFaceTot do i = 1,3*nIntFaceTot
write(6,'(1x,2(e14.8,1x))')state%p(i) write(6,'(1x,2(e15.8,1x))')state%p(i)
enddo enddo
write(6,*)' ' write(6,*)' '
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
@ -410,11 +410,11 @@ function homogenization_RGC_updateState(&
if (debug_verbosity == 4) then if (debug_verbosity == 4) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
do iGrain = 1,nGrain do iGrain = 1,nGrain
write(6,'(1x,a30,1x,i3,1x,a4,3(1x,e14.8))')'Mismatch magnitude of grain(',iGrain,') :',NN(1,iGrain),NN(2,iGrain),NN(3,iGrain) write(6,'(1x,a30,1x,i3,1x,a4,3(1x,e15.8))')'Mismatch magnitude of grain(',iGrain,') :',NN(1,iGrain),NN(2,iGrain),NN(3,iGrain)
write(6,*)' ' write(6,*)' '
write(6,'(1x,a30,1x,i3)')'Stress and penalties of grain: ',iGrain write(6,'(1x,a30,1x,i3)')'Stress and penalties of grain: ',iGrain
do i = 1,3 do i = 1,3
write(6,'(1x,3(e14.8,1x),1x,3(e14.8,1x),1x,3(e14.8,1x))')(P(i,j,iGrain), j = 1,3), & write(6,'(1x,3(e15.8,1x),1x,3(e15.8,1x),1x,3(e15.8,1x))')(P(i,j,iGrain), j = 1,3), &
(R(i,j,iGrain), j = 1,3), & (R(i,j,iGrain), j = 1,3), &
(D(i,j,iGrain), j = 1,3) (D(i,j,iGrain), j = 1,3)
enddo enddo
@ -459,7 +459,7 @@ function homogenization_RGC_updateState(&
if (debug_verbosity == 4) then if (debug_verbosity == 4) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(1x,a30,1x,i3)')'Traction at interface: ',iNum write(6,'(1x,a30,1x,i3)')'Traction at interface: ',iNum
write(6,'(1x,3(e14.8,1x))')(tract(iNum,j), j = 1,3) write(6,'(1x,3(e15.8,1x))')(tract(iNum,j), j = 1,3)
write(6,*)' ' write(6,*)' '
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
@ -478,9 +478,9 @@ function homogenization_RGC_updateState(&
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(1x,a)')' ' write(6,'(1x,a)')' '
write(6,'(1x,a,1x,i2,1x,i4)')'RGC residual check ...',ip,el write(6,'(1x,a,1x,i2,1x,i4)')'RGC residual check ...',ip,el
write(6,'(1x,a15,1x,e14.8,1x,a7,i3,1x,a12,i2,i2)')'Max stress: ',stresMax, & write(6,'(1x,a15,1x,e15.8,1x,a7,i3,1x,a12,i2,i2)')'Max stress: ',stresMax, &
'@ grain',stresLoc(3),'in component',stresLoc(1),stresLoc(2) '@ grain',stresLoc(3),'in component',stresLoc(1),stresLoc(2)
write(6,'(1x,a15,1x,e14.8,1x,a7,i3,1x,a12,i2)')'Max residual: ',residMax, & write(6,'(1x,a15,1x,e15.8,1x,a7,i3,1x,a12,i2)')'Max residual: ',residMax, &
'@ iface',residLoc(1),'in direction',residLoc(2) '@ iface',residLoc(1),'in direction',residLoc(2)
call flush(6) call flush(6)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
@ -523,15 +523,15 @@ function homogenization_RGC_updateState(&
if (debug_verbosity == 4 .and. debug_e == el .and. debug_i == ip) then if (debug_verbosity == 4 .and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(1x,a30,1x,e14.8)')'Constitutive work: ',constitutiveWork write(6,'(1x,a30,1x,e15.8)')'Constitutive work: ',constitutiveWork
write(6,'(1x,a30,3(1x,e14.8))')'Magnitude mismatch: ',sum(NN(1,:))/dble(nGrain), & write(6,'(1x,a30,3(1x,e15.8))')'Magnitude mismatch: ',sum(NN(1,:))/dble(nGrain), &
sum(NN(2,:))/dble(nGrain), & sum(NN(2,:))/dble(nGrain), &
sum(NN(3,:))/dble(nGrain) sum(NN(3,:))/dble(nGrain)
write(6,'(1x,a30,1x,e14.8)')'Penalty energy: ',penaltyEnergy write(6,'(1x,a30,1x,e15.8)')'Penalty energy: ',penaltyEnergy
write(6,'(1x,a30,1x,e14.8)')'Volume discrepancy: ',volDiscrep write(6,'(1x,a30,1x,e15.8)')'Volume discrepancy: ',volDiscrep
write(6,*)'' write(6,*)''
write(6,'(1x,a30,1x,e14.8)')'Maximum relaxation rate: ',maxval(abs(drelax))/dt write(6,'(1x,a30,1x,e15.8)')'Maximum relaxation rate: ',maxval(abs(drelax))/dt
write(6,'(1x,a30,1x,e14.8)')'Average relaxation rate: ',sum(abs(drelax))/dt/dble(3*nIntFaceTot) write(6,'(1x,a30,1x,e15.8)')'Average relaxation rate: ',sum(abs(drelax))/dt/dble(3*nIntFaceTot)
write(6,*)'' write(6,*)''
call flush(6) call flush(6)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
@ -619,7 +619,7 @@ function homogenization_RGC_updateState(&
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Jacobian matrix of stress' write(6,'(1x,a30)')'Jacobian matrix of stress'
do i = 1,3*nIntFaceTot do i = 1,3*nIntFaceTot
write(6,'(1x,100(e10.4,1x))')(smatrix(i,j), j = 1,3*nIntFaceTot) write(6,'(1x,100(e11.4,1x))')(smatrix(i,j), j = 1,3*nIntFaceTot)
enddo enddo
write(6,*)' ' write(6,*)' '
call flush(6) call flush(6)
@ -675,7 +675,7 @@ function homogenization_RGC_updateState(&
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Jacobian matrix of penalty' write(6,'(1x,a30)')'Jacobian matrix of penalty'
do i = 1,3*nIntFaceTot do i = 1,3*nIntFaceTot
write(6,'(1x,100(e10.4,1x))')(pmatrix(i,j), j = 1,3*nIntFaceTot) write(6,'(1x,100(e11.4,1x))')(pmatrix(i,j), j = 1,3*nIntFaceTot)
enddo enddo
write(6,*)' ' write(6,*)' '
call flush(6) call flush(6)
@ -695,7 +695,7 @@ function homogenization_RGC_updateState(&
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Jacobian matrix of penalty' write(6,'(1x,a30)')'Jacobian matrix of penalty'
do i = 1,3*nIntFaceTot do i = 1,3*nIntFaceTot
write(6,'(1x,100(e10.4,1x))')(rmatrix(i,j), j = 1,3*nIntFaceTot) write(6,'(1x,100(e11.4,1x))')(rmatrix(i,j), j = 1,3*nIntFaceTot)
enddo enddo
write(6,*)' ' write(6,*)' '
call flush(6) call flush(6)
@ -709,7 +709,7 @@ function homogenization_RGC_updateState(&
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Jacobian matrix (total)' write(6,'(1x,a30)')'Jacobian matrix (total)'
do i = 1,3*nIntFaceTot do i = 1,3*nIntFaceTot
write(6,'(1x,100(e10.4,1x))')(jmatrix(i,j), j = 1,3*nIntFaceTot) write(6,'(1x,100(e11.4,1x))')(jmatrix(i,j), j = 1,3*nIntFaceTot)
enddo enddo
write(6,*)' ' write(6,*)' '
call flush(6) call flush(6)
@ -728,7 +728,7 @@ function homogenization_RGC_updateState(&
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Jacobian inverse' write(6,'(1x,a30)')'Jacobian inverse'
do i = 1,3*nIntFaceTot do i = 1,3*nIntFaceTot
write(6,'(1x,100(e10.4,1x))')(jnverse(i,j), j = 1_pInt,3_pInt*nIntFaceTot) write(6,'(1x,100(e11.4,1x))')(jnverse(i,j), j = 1_pInt,3_pInt*nIntFaceTot)
enddo enddo
write(6,*)' ' write(6,*)' '
call flush(6) call flush(6)
@ -748,7 +748,7 @@ function homogenization_RGC_updateState(&
homogenization_RGC_updateState = (/.true.,.false./) homogenization_RGC_updateState = (/.true.,.false./)
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(1x,a,1x,i3,1x,a,1x,i3,1x,a)')'RGC_updateState: ip',ip,'| el',el,'enforces cutback' write(6,'(1x,a,1x,i3,1x,a,1x,i3,1x,a)')'RGC_updateState: ip',ip,'| el',el,'enforces cutback'
write(6,'(1x,a,1x,e14.8)')'due to large relaxation change =',maxval(abs(drelax)) write(6,'(1x,a,1x,e15.8)')'due to large relaxation change =',maxval(abs(drelax))
call flush(6) call flush(6)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
@ -758,7 +758,7 @@ function homogenization_RGC_updateState(&
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Returned state: ' write(6,'(1x,a30)')'Returned state: '
do i = 1,3*nIntFaceTot do i = 1,3*nIntFaceTot
write(6,'(1x,2(e14.8,1x))')state%p(i) write(6,'(1x,2(e15.8,1x))')state%p(i)
enddo enddo
write(6,*)' ' write(6,*)' '
call flush(6) call flush(6)
@ -810,7 +810,7 @@ subroutine homogenization_RGC_averageStressAndItsTangent(&
dPdF99 = math_Plain3333to99(dPdF(1:3,1:3,1:3,1:3,iGrain)) dPdF99 = math_Plain3333to99(dPdF(1:3,1:3,1:3,1:3,iGrain))
write(6,'(1x,a30,1x,i3)')'Stress tangent of grain: ',iGrain write(6,'(1x,a30,1x,i3)')'Stress tangent of grain: ',iGrain
do i = 1,9 do i = 1,9
write(6,'(1x,(e14.8,1x))') (dPdF99(i,j), j = 1,9) write(6,'(1x,(e15.8,1x))') (dPdF99(i,j), j = 1,9)
enddo enddo
write(6,*)' ' write(6,*)' '
enddo enddo
@ -955,7 +955,7 @@ subroutine homogenization_RGC_stressPenalty(&
!* Debugging the surface correction factor !* Debugging the surface correction factor
! if (ip == 1 .and. el == 1) then ! if (ip == 1 .and. el == 1) then
! write(6,'(1x,a20,2(1x,i3))')'Correction factor: ',ip,el ! write(6,'(1x,a20,2(1x,i3))')'Correction factor: ',ip,el
! write(6,'(1x,3(e10.4,1x))')(surfCorr(i), i = 1,3) ! write(6,'(1x,3(e11.4,1x))')(surfCorr(i), i = 1,3)
! endif ! endif
!* ------------------------------------------------------------------------------------------------------------- !* -------------------------------------------------------------------------------------------------------------
@ -1005,9 +1005,9 @@ subroutine homogenization_RGC_stressPenalty(&
! if (ip == 1 .and. el == 1) then ! if (ip == 1 .and. el == 1) then
! write(6,'(1x,a20,i2,1x,a20,1x,i3)')'Mismatch to face: ',intFace(1),'neighbor grain: ',iGNghb ! write(6,'(1x,a20,i2,1x,a20,1x,i3)')'Mismatch to face: ',intFace(1),'neighbor grain: ',iGNghb
! do i = 1,3 ! do i = 1,3
! write(6,'(1x,3(e10.4,1x))')(nDef(i,j), j = 1,3) ! write(6,'(1x,3(e11.4,1x))')(nDef(i,j), j = 1,3)
! enddo ! enddo
! write(6,'(1x,a20,e10.4))')'with magnitude: ',nDefNorm ! write(6,'(1x,a20,e11.4))')'with magnitude: ',nDefNorm
! endif ! endif
!* Compute the stress penalty of all interfaces !* Compute the stress penalty of all interfaces
@ -1030,7 +1030,7 @@ subroutine homogenization_RGC_stressPenalty(&
! if (ip == 1 .and. el == 1) then ! if (ip == 1 .and. el == 1) then
! write(6,'(1x,a20,i2)')'Penalty of grain: ',iGrain ! write(6,'(1x,a20,i2)')'Penalty of grain: ',iGrain
! do i = 1,3 ! do i = 1,3
! write(6,'(1x,3(e10.4,1x))')(rPen(i,j,iGrain), j = 1,3) ! write(6,'(1x,3(e11.4,1x))')(rPen(i,j,iGrain), j = 1,3)
! enddo ! enddo
! endif ! endif
@ -1090,7 +1090,7 @@ subroutine homogenization_RGC_volumePenalty(&
! if (ip == 1 .and. el == 1) then ! if (ip == 1 .and. el == 1) then
! write(6,'(1x,a30,i2)')'Volume penalty of grain: ',iGrain ! write(6,'(1x,a30,i2)')'Volume penalty of grain: ',iGrain
! do i = 1,3 ! do i = 1,3
! write(6,'(1x,3(e10.4,1x))')(vPen(i,j,iGrain), j = 1,3) ! write(6,'(1x,3(e11.4,1x))')(vPen(i,j,iGrain), j = 1,3)
! enddo ! enddo
! endif ! endif
@ -1231,7 +1231,7 @@ function homogenization_RGC_interfaceNormal(&
!* Get the normal of the interface, identified from the value of intFace(1) !* Get the normal of the interface, identified from the value of intFace(1)
homogenization_RGC_interfaceNormal = 0.0_pReal homogenization_RGC_interfaceNormal = 0.0_pReal
nPos = abs(intFace(1)) ! identify the position of the interface in global state array nPos = abs(intFace(1)) ! identify the position of the interface in global state array
homogenization_RGC_interfaceNormal(nPos) = intFace(1)/abs(intFace(1)) ! get the normal vector w.r.t. cluster axis homogenization_RGC_interfaceNormal(nPos) = real(intFace(1)/abs(intFace(1)),pReal) ! get the normal vector w.r.t. cluster axis
homogenization_RGC_interfaceNormal = & homogenization_RGC_interfaceNormal = &
math_mul33x3(homogenization_RGC_orientation(:,:,ip,el),homogenization_RGC_interfaceNormal) math_mul33x3(homogenization_RGC_orientation(:,:,ip,el),homogenization_RGC_interfaceNormal)
@ -1239,7 +1239,7 @@ function homogenization_RGC_interfaceNormal(&
! if (ip == 1 .and. el == 1) then ! if (ip == 1 .and. el == 1) then
! write(6,'(1x,a32,3(1x,i3))')'Interface normal: ',intFace(1) ! write(6,'(1x,a32,3(1x,i3))')'Interface normal: ',intFace(1)
! write(6,'(1x,3(e14.8,1x))')(nVect(i), i = 1,3) ! write(6,'(1x,3(e15.8,1x))')(nVect(i), i = 1,3)
! write(6,*)' ' ! write(6,*)' '
! call flush(6) ! call flush(6)
! endif ! endif

227
code/kdtree2.f90
View File

@ -64,7 +64,7 @@ module kdtree2_priority_queue_module
! There are heap_size active elements. ! There are heap_size active elements.
! Assumes the allocation is always sufficient. Will NOT increase it ! Assumes the allocation is always sufficient. Will NOT increase it
! to match. ! to match.
integer(pInt) :: heap_size = 0 integer(pInt) :: heap_size = 0_pInt
type(kdtree2_result), pointer :: elems(:) type(kdtree2_result), pointer :: elems(:)
end type pq end type pq
@ -97,14 +97,14 @@ contains
type(pq) :: res type(pq) :: res
! !
! !
integer(pInt) :: nalloc integer :: nalloc
nalloc = size(results_in,1) nalloc = size(results_in,1)
if (nalloc .lt. 1) then if (nalloc .lt. 1) then
write (*,*) 'PQ_CREATE: error, input arrays must be allocated.' write (*,*) 'PQ_CREATE: error, input arrays must be allocated.'
end if end if
res%elems => results_in res%elems => results_in
res%heap_size = 0 res%heap_size = 0_pInt
return return
end function pq_create end function pq_create
@ -160,8 +160,8 @@ contains
i = i_in i = i_in
bigloop: do bigloop: do
l = 2*i ! left(i) l = 2_pInt*i ! left(i)
r = l+1 ! right(i) r = l+1_pInt ! right(i)
! !
! set 'largest' to the index of either i, l, r ! set 'largest' to the index of either i, l, r
! depending on whose priority is largest. ! depending on whose priority is largest.
@ -296,11 +296,11 @@ bigloop: do
! write (*,*) 'PQ_INSERT: error, attempt made to insert element on full PQ' ! write (*,*) 'PQ_INSERT: error, attempt made to insert element on full PQ'
! stop ! stop
! else ! else
a%heap_size = a%heap_size + 1 a%heap_size = a%heap_size + 1_pInt
i = a%heap_size i = a%heap_size
do while (i .gt. 1) do while (i .gt. 1)
isparent = int(i/2) isparent = int(i/2_pInt, pInt) !needed casting?
parentdis = a%elems(isparent)%dis parentdis = a%elems(isparent)%dis
if (dis .gt. parentdis) then if (dis .gt. parentdis) then
! move what was in i's parent into i. ! move what was in i's parent into i.
@ -339,13 +339,13 @@ bigloop: do
e = a%elems(i) e = a%elems(i)
parent = i parent = i
child = 2*i child = 2_pInt*i
N = a%heap_size N = a%heap_size
do while (child .le. N) do while (child .le. N)
if (child .lt. N) then if (child .lt. N) then
if (a%elems(child)%dis .lt. a%elems(child+1)%dis) then if (a%elems(child)%dis .lt. a%elems(child+1)%dis) then
child = child+1 child = child+1_pInt
endif endif
endif endif
prichild = a%elems(child)%dis prichild = a%elems(child)%dis
@ -355,7 +355,7 @@ bigloop: do
! move child into parent. ! move child into parent.
a%elems(parent) = a%elems(child) a%elems(parent) = a%elems(child)
parent = child parent = child
child = 2*parent child = 2_pInt*parent
end if end if
end do end do
a%elems(parent) = e a%elems(parent) = e
@ -384,8 +384,8 @@ bigloop: do
if (.true.) then if (.true.) then
N=a%heap_size N=a%heap_size
if (N .ge. 1) then if (N .ge. 1) then
parent =1 parent =1_pInt
child=2 child=2_pInt
loop: do while (child .le. N) loop: do while (child .le. N)
prichild = a%elems(child)%dis prichild = a%elems(child)%dis
@ -396,9 +396,9 @@ bigloop: do
! !
if (child .lt. N) then if (child .lt. N) then
prichildp1 = a%elems(child+1)%dis prichildp1 = a%elems(child+1_pInt)%dis
if (prichild .lt. prichildp1) then if (prichild .lt. prichildp1) then
child = child+1 child = child+1_pInt
prichild = prichildp1 prichild = prichildp1
endif endif
endif endif
@ -411,7 +411,7 @@ bigloop: do
! move child into parent. ! move child into parent.
a%elems(parent) = a%elems(child) a%elems(parent) = a%elems(child)
parent = child parent = child
child = 2*parent child = 2_pInt*parent
end if end if
end do loop end do loop
a%elems(parent)%dis = dis a%elems(parent)%dis = dis
@ -449,7 +449,7 @@ bigloop: do
! swap the item to be deleted with the last element ! swap the item to be deleted with the last element
! and shorten heap by one. ! and shorten heap by one.
a%elems(i) = a%elems(a%heap_size) a%elems(i) = a%elems(a%heap_size)
a%heap_size = a%heap_size - 1 a%heap_size = a%heap_size - 1_pInt
call heapify(a,i) call heapify(a,i)
@ -469,10 +469,10 @@ module kdtree2_module
! This module is identical to 'kd_tree', except that the order ! This module is identical to 'kd_tree', except that the order
! of subscripts is reversed in the data file. ! of subscripts is reversed in the data file.
! In otherwords for an embedding of N D-dimensional vectors, the ! In otherwords for an embedding of N D-dimensional vectors, the
! data file is here, in natural Fortran order data(1:D, 1:N) ! data file is here, in natural Fortran order myData(1:D, 1:N)
! because Fortran lays out columns first, ! because Fortran lays out columns first,
! !
! whereas conventionally (C-style) it is data(1:N,1:D) ! whereas conventionally (C-style) it is myData(1:N,1:D)
! as in the original kd_tree module. ! as in the original kd_tree module.
! !
!-------------DATA TYPE, CREATION, DELETION--------------------- !-------------DATA TYPE, CREATION, DELETION---------------------
@ -498,7 +498,7 @@ module kdtree2_module
!---------------------------------------------------------------- !----------------------------------------------------------------
integer(pInt), parameter :: bucket_size = 12 integer(pInt), parameter :: bucket_size = 12_pInt
! The maximum number of points to keep in a terminal node. ! The maximum number of points to keep in a terminal node.
type interval type interval
@ -525,7 +525,7 @@ module kdtree2_module
type :: kdtree2 type :: kdtree2
! Global information about the tree, one per tree ! Global information about the tree, one per tree
integer(pInt) :: dimen=0, n=0 integer(pInt) :: dimen=0_pInt, n=0_pInt
! dimensionality and total # of points ! dimensionality and total # of points
real(pReal), pointer :: the_data(:,:) => null() real(pReal), pointer :: the_data(:,:) => null()
! pointer to the actual data array ! pointer to the actual data array
@ -570,7 +570,7 @@ module kdtree2_module
integer(pInt) :: dimen integer(pInt) :: dimen
integer(pInt) :: nn, nfound integer(pInt) :: nn, nfound
real(pReal) :: ballsize real(pReal) :: ballsize
integer(pInt) :: centeridx=999, correltime=9999 integer(pInt) :: centeridx=999_pInt, correltime=9999_pInt
! exclude points within 'correltime' of 'centeridx', iff centeridx >= 0 ! exclude points within 'correltime' of 'centeridx', iff centeridx >= 0
integer(pInt) :: nalloc ! how much allocated for results(:)? integer(pInt) :: nalloc ! how much allocated for results(:)?
logical :: rearrange ! are the data rearranged or original? logical :: rearrange ! are the data rearranged or original?
@ -579,7 +579,7 @@ module kdtree2_module
real(pReal), pointer :: qv(:) ! query vector real(pReal), pointer :: qv(:) ! query vector
type(kdtree2_result), pointer :: results(:) ! results type(kdtree2_result), pointer :: results(:) ! results
type(pq) :: pq type(pq) :: pq
real(pReal), pointer :: data(:,:) ! temp pointer to data real(pReal), pointer :: myData(:,:) ! temp pointer to data
integer(pInt), pointer :: ind(:) ! temp pointer to indexes integer(pInt), pointer :: ind(:) ! temp pointer to indexes
end type tree_search_record end type tree_search_record
@ -590,7 +590,7 @@ module kdtree2_module
contains contains
function kdtree2_create(input_data,dim,sort,rearrange) result (mr) function kdtree2_create(input_data,myDim,sort,rearrange) result (mr)
! !
! create the actual tree structure, given an input array of data. ! create the actual tree structure, given an input array of data.
! !
@ -598,9 +598,9 @@ contains
! THIS IS THE REVERSE OF THE PREVIOUS VERSION OF THIS MODULE. ! THIS IS THE REVERSE OF THE PREVIOUS VERSION OF THIS MODULE.
! The reason for it is cache friendliness, improving performance. ! The reason for it is cache friendliness, improving performance.
! !
! Optional arguments: If 'dim' is specified, then the tree ! Optional arguments: If 'myDim' is specified, then the tree
! will only search the first 'dim' components ! will only search the first 'myDim' components
! of input_data, otherwise, dim is inferred ! of input_data, otherwise, myDim is inferred
! from SIZE(input_data,1). ! from SIZE(input_data,1).
! !
! if sort .eqv. .true. then output results ! if sort .eqv. .true. then output results
@ -615,7 +615,7 @@ contains
! !
! .. Function Return Cut_value .. ! .. Function Return Cut_value ..
type (kdtree2), pointer :: mr type (kdtree2), pointer :: mr
integer(pInt), intent(in), optional :: dim integer(pInt), intent(in), optional :: myDim
logical, intent(in), optional :: sort logical, intent(in), optional :: sort
logical, intent(in), optional :: rearrange logical, intent(in), optional :: rearrange
! .. ! ..
@ -628,19 +628,19 @@ contains
mr%the_data => input_data mr%the_data => input_data
! pointer assignment ! pointer assignment
if (present(dim)) then if (present(myDim)) then
mr%dimen = dim mr%dimen = myDim
else else
mr%dimen = size(input_data,1) mr%dimen = int(size(input_data,1), pInt) ! size returns default integer
end if end if
mr%n = size(input_data,2) mr%n = int(size(input_data,2), pInt) ! size returns default integer
if (mr%dimen > mr%n) then if (mr%dimen > mr%n) then
! unlikely to be correct ! unlikely to be correct
write (*,*) 'KD_TREE_TRANS: likely user error.' write (*,*) 'KD_TREE_TRANS: likely user error.'
write (*,*) 'KD_TREE_TRANS: You passed in matrix with D=',mr%dimen write (*,*) 'KD_TREE_TRANS: You passed in matrix with D=',mr%dimen
write (*,*) 'KD_TREE_TRANS: and N=',mr%n write (*,*) 'KD_TREE_TRANS: and N=',mr%n
write (*,*) 'KD_TREE_TRANS: note, that new format is data(1:D,1:N)' write (*,*) 'KD_TREE_TRANS: note, that new format is myData(1:D,1:N)'
write (*,*) 'KD_TREE_TRANS: with usually N >> D. If N =approx= D, then a k-d tree' write (*,*) 'KD_TREE_TRANS: with usually N >> D. If N =approx= D, then a k-d tree'
write (*,*) 'KD_TREE_TRANS: is not an appropriate data structure.' write (*,*) 'KD_TREE_TRANS: is not an appropriate data structure.'
stop stop
@ -662,7 +662,7 @@ contains
if (mr%rearrange) then if (mr%rearrange) then
allocate(mr%rearranged_data(mr%dimen,mr%n)) allocate(mr%rearranged_data(mr%dimen,mr%n))
do i=1,mr%n do i=1_pInt,mr%n
mr%rearranged_data(:,i) = mr%the_data(:, & mr%rearranged_data(:,i) = mr%the_data(:, &
mr%ind(i)) mr%ind(i))
enddo enddo
@ -679,7 +679,7 @@ contains
type(tree_node), pointer :: dummy => null() type(tree_node), pointer :: dummy => null()
! .. ! ..
allocate (tp%ind(tp%n)) allocate (tp%ind(tp%n))
forall (j=1:tp%n) forall (j=1_pInt:tp%n)
tp%ind(j) = j tp%ind(j) = j
end forall end forall
tp%root => build_tree_for_range(tp,1_pInt,tp%n, dummy) tp%root => build_tree_for_range(tp,1_pInt,tp%n, dummy)
@ -729,10 +729,10 @@ contains
! !
! always compute true bounding box for terminal nodes. ! always compute true bounding box for terminal nodes.
! !
do i=1,dimen do i=1_pInt,dimen
call spread_in_coordinate(tp,i,l,u,res%box(i)) call spread_in_coordinate(tp,i,l,u,res%box(i))
end do end do
res%cut_dim = 0 res%cut_dim = 0_pInt
res%cut_val = 0.0_pReal res%cut_val = 0.0_pReal
res%l = l res%l = l
res%u = u res%u = u
@ -764,7 +764,7 @@ contains
end do end do
c = maxloc(res%box(1:dimen)%upper-res%box(1:dimen)%lower,1) c = int(maxloc(res%box(1:dimen)%upper-res%box(1:dimen)%lower,1), pInt)
! !
! c is the identity of which coordinate has the greatest spread. ! c is the identity of which coordinate has the greatest spread.
! !
@ -867,11 +867,11 @@ contains
do while (lb < rb) do while (lb < rb)
if ( v(c,ind(lb)) <= alpha ) then if ( v(c,ind(lb)) <= alpha ) then
! it is good where it is. ! it is good where it is.
lb = lb+1 lb = lb+1_pInt
else else
! swap it with rb. ! swap it with rb.
tmp = ind(lb); ind(lb) = ind(rb); ind(rb) = tmp tmp = ind(lb); ind(lb) = ind(rb); ind(rb) = tmp
rb = rb-1 rb = rb-1_pInt
endif endif
end do end do
@ -879,7 +879,7 @@ contains
if (v(c,ind(lb)) <= alpha) then if (v(c,ind(lb)) <= alpha) then
res = lb res = lb
else else
res = lb-1 res = lb-1_pInt
endif endif
end function select_on_coordinate_value end function select_on_coordinate_value
@ -901,9 +901,9 @@ contains
do while (l<u) do while (l<u)
t = ind(l) t = ind(l)
m = l m = l
do i = l + 1, u do i = l + 1_pInt, u
if (v(c,ind(i))<v(c,t)) then if (v(c,ind(i))<v(c,t)) then
m = m + 1 m = m + 1_pInt
s = ind(m) s = ind(m)
ind(m) = ind(i) ind(m) = ind(i)
ind(i) = s ind(i) = s
@ -912,8 +912,8 @@ contains
s = ind(l) s = ind(l)
ind(l) = ind(m) ind(l) = ind(m)
ind(m) = s ind(m) = s
if (m<=k) l = m + 1 if (m<=k) l = m + 1_pInt
if (m>=k) u = m - 1 if (m>=k) u = m - 1_pInt
end do end do
end subroutine select_on_coordinate end subroutine select_on_coordinate
@ -944,8 +944,8 @@ contains
ulocal = u ulocal = u
do i = l + 2, ulocal, 2 do i = l + 2_pInt, ulocal, 2_pInt
lmin = v(c,ind(i-1)) lmin = v(c,ind(i-1_pInt))
lmax = v(c,ind(i)) lmax = v(c,ind(i))
if (lmin>lmax) then if (lmin>lmax) then
t = lmin t = lmin
@ -1022,9 +1022,9 @@ contains
sr%ballsize = huge(1.0_pReal) sr%ballsize = huge(1.0_pReal)
sr%qv => qv sr%qv => qv
sr%nn = nn sr%nn = nn
sr%nfound = 0 sr%nfound = 0_pInt
sr%centeridx = -1 sr%centeridx = -1_pInt
sr%correltime = 0 sr%correltime = 0_pInt
sr%overflow = .false. sr%overflow = .false.
sr%results => results sr%results => results
@ -1034,9 +1034,9 @@ contains
sr%ind => tp%ind sr%ind => tp%ind
sr%rearrange = tp%rearrange sr%rearrange = tp%rearrange
if (tp%rearrange) then if (tp%rearrange) then
sr%Data => tp%rearranged_data sr%myData => tp%rearranged_data
else else
sr%Data => tp%the_data sr%myData => tp%the_data
endif endif
sr%dimen = tp%dimen sr%dimen = tp%dimen
@ -1067,7 +1067,7 @@ contains
sr%correltime = correltime sr%correltime = correltime
sr%nn = nn sr%nn = nn
sr%nfound = 0 sr%nfound = 0_pInt
sr%dimen = tp%dimen sr%dimen = tp%dimen
sr%nalloc = nn sr%nalloc = nn
@ -1078,9 +1078,9 @@ contains
sr%rearrange = tp%rearrange sr%rearrange = tp%rearrange
if (sr%rearrange) then if (sr%rearrange) then
sr%Data => tp%rearranged_data sr%myData => tp%rearranged_data
else else
sr%Data => tp%the_data sr%myData => tp%the_data
endif endif
call validate_query_storage(nn) call validate_query_storage(nn)
@ -1116,10 +1116,10 @@ contains
! !
sr%qv => qv sr%qv => qv
sr%ballsize = r2 sr%ballsize = r2
sr%nn = 0 ! flag for fixed ball search sr%nn = 0_pInt ! flag for fixed ball search
sr%nfound = 0 sr%nfound = 0_pInt
sr%centeridx = -1 sr%centeridx = -1_pInt
sr%correltime = 0 sr%correltime = 0_pInt
sr%results => results sr%results => results
@ -1130,9 +1130,9 @@ contains
sr%rearrange= tp%rearrange sr%rearrange= tp%rearrange
if (tp%rearrange) then if (tp%rearrange) then
sr%Data => tp%rearranged_data sr%myData => tp%rearranged_data
else else
sr%Data => tp%the_data sr%myData => tp%the_data
endif endif
sr%dimen = tp%dimen sr%dimen = tp%dimen
@ -1176,8 +1176,8 @@ contains
allocate (sr%qv(tp%dimen)) allocate (sr%qv(tp%dimen))
sr%qv = tp%the_data(:,idxin) ! copy the vector sr%qv = tp%the_data(:,idxin) ! copy the vector
sr%ballsize = r2 sr%ballsize = r2
sr%nn = 0 ! flag for fixed r search sr%nn = 0_pInt ! flag for fixed r search
sr%nfound = 0 sr%nfound = 0_pInt
sr%centeridx = idxin sr%centeridx = idxin
sr%correltime = correltime sr%correltime = correltime
@ -1195,9 +1195,9 @@ contains
sr%rearrange = tp%rearrange sr%rearrange = tp%rearrange
if (tp%rearrange) then if (tp%rearrange) then
sr%Data => tp%rearranged_data sr%myData => tp%rearranged_data
else else
sr%Data => tp%the_data sr%myData => tp%the_data
endif endif
sr%rearrange = tp%rearrange sr%rearrange = tp%rearrange
sr%dimen = tp%dimen sr%dimen = tp%dimen
@ -1236,21 +1236,21 @@ contains
sr%qv => qv sr%qv => qv
sr%ballsize = r2 sr%ballsize = r2
sr%nn = 0 ! flag for fixed r search sr%nn = 0_pInt ! flag for fixed r search
sr%nfound = 0 sr%nfound = 0_pInt
sr%centeridx = -1 sr%centeridx = -1_pInt
sr%correltime = 0 sr%correltime = 0_pInt
nullify(sr%results) ! for some reason, FTN 95 chokes on '=> null()' nullify(sr%results) ! for some reason, FTN 95 chokes on '=> null()'
sr%nalloc = 0 ! we do not allocate any storage but that's OK sr%nalloc = 0_pInt ! we do not allocate any storage but that's OK
! for counting. ! for counting.
sr%ind => tp%ind sr%ind => tp%ind
sr%rearrange = tp%rearrange sr%rearrange = tp%rearrange
if (tp%rearrange) then if (tp%rearrange) then
sr%Data => tp%rearranged_data sr%myData => tp%rearranged_data
else else
sr%Data => tp%the_data sr%myData => tp%the_data
endif endif
sr%dimen = tp%dimen sr%dimen = tp%dimen
@ -1285,22 +1285,22 @@ contains
sr%qv = tp%the_data(:,idxin) sr%qv = tp%the_data(:,idxin)
sr%ballsize = r2 sr%ballsize = r2
sr%nn = 0 ! flag for fixed r search sr%nn = 0_pInt ! flag for fixed r search
sr%nfound = 0 sr%nfound = 0_pInt
sr%centeridx = idxin sr%centeridx = idxin
sr%correltime = correltime sr%correltime = correltime
nullify(sr%results) nullify(sr%results)
sr%nalloc = 0 ! we do not allocate any storage but that's OK sr%nalloc = 0_pInt ! we do not allocate any storage but that's OK
! for counting. ! for counting.
sr%ind => tp%ind sr%ind => tp%ind
sr%rearrange = tp%rearrange sr%rearrange = tp%rearrange
if (sr%rearrange) then if (sr%rearrange) then
sr%Data => tp%rearranged_data sr%myData => tp%rearranged_data
else else
sr%Data => tp%the_data sr%myData => tp%the_data
endif endif
sr%dimen = tp%dimen sr%dimen = tp%dimen
@ -1324,7 +1324,7 @@ contains
! !
integer(pInt), intent(in) :: n integer(pInt), intent(in) :: n
if (size(sr%results,1) .lt. n) then if (int(size(sr%results,1),pInt) .lt. n) then
write (*,*) 'KD_TREE_TRANS: you did not provide enough storage for results(1:n)' write (*,*) 'KD_TREE_TRANS: you did not provide enough storage for results(1:n)'
stop stop
return return
@ -1408,7 +1408,7 @@ contains
! check will also be false. ! check will also be false.
! !
box => node%box(1:) box => node%box(1:)
do i=1,sr%dimen do i=1_pInt,sr%dimen
if (i .ne. cut_dim) then if (i .ne. cut_dim) then
dis = dis + dis2_from_bnd(qv(i),box(i)%lower,box(i)%upper) dis = dis + dis2_from_bnd(qv(i),box(i)%lower,box(i)%upper)
if (dis > ballsize) then if (dis > ballsize) then
@ -1486,7 +1486,7 @@ contains
! !
real(pReal), pointer :: qv(:) real(pReal), pointer :: qv(:)
integer(pInt), pointer :: ind(:) integer(pInt), pointer :: ind(:)
real(pReal), pointer :: data(:,:) real(pReal), pointer :: myData(:,:)
! !
integer(pInt) :: dimen, i, indexofi, k, centeridx, correltime integer(pInt) :: dimen, i, indexofi, k, centeridx, correltime
real(pReal) :: ballsize, sd, newpri real(pReal) :: ballsize, sd, newpri
@ -1505,7 +1505,7 @@ contains
ballsize = sr%ballsize ballsize = sr%ballsize
rearrange = sr%rearrange rearrange = sr%rearrange
ind => sr%ind(1:) ind => sr%ind(1:)
data => sr%Data(1:,1:) myData => sr%myData(1:,1:)
centeridx = sr%centeridx centeridx = sr%centeridx
correltime = sr%correltime correltime = sr%correltime
@ -1517,7 +1517,7 @@ contains
if (rearrange) then if (rearrange) then
sd = 0.0_pReal sd = 0.0_pReal
do k = 1_pInt,dimen do k = 1_pInt,dimen
sd = sd + (data(k,i) - qv(k))**2.0_pReal sd = sd + (myData(k,i) - qv(k))**2.0_pReal
if (sd>ballsize) cycle mainloop if (sd>ballsize) cycle mainloop
end do end do
indexofi = ind(i) ! only read it if we have not broken out indexofi = ind(i) ! only read it if we have not broken out
@ -1525,7 +1525,7 @@ contains
indexofi = ind(i) indexofi = ind(i)
sd = 0.0_pReal sd = 0.0_pReal
do k = 1_pInt,dimen do k = 1_pInt,dimen
sd = sd + (data(k,indexofi) - qv(k))**2.0_pReal sd = sd + (myData(k,indexofi) - qv(k))**2.0_pReal
if (sd>ballsize) cycle mainloop if (sd>ballsize) cycle mainloop
end do end do
endif endif
@ -1557,7 +1557,7 @@ contains
! !
! add this point unconditionally to fill list. ! add this point unconditionally to fill list.
! !
sr%nfound = sr%nfound +1 sr%nfound = sr%nfound +1_pInt
newpri = pq_insert(pqp,sd,indexofi) newpri = pq_insert(pqp,sd,indexofi)
if (sr%nfound .eq. sr%nn) ballsize = newpri if (sr%nfound .eq. sr%nn) ballsize = newpri
! we have just filled the working list. ! we have just filled the working list.
@ -1592,7 +1592,7 @@ contains
! !
real(pReal), pointer :: qv(:) real(pReal), pointer :: qv(:)
integer(pInt), pointer :: ind(:) integer(pInt), pointer :: ind(:)
real(pReal), pointer :: data(:,:) real(pReal), pointer :: myData(:,:)
! !
integer(pInt) :: nfound integer(pInt) :: nfound
integer(pInt) :: dimen, i, indexofi, k integer(pInt) :: dimen, i, indexofi, k
@ -1608,7 +1608,7 @@ contains
ballsize = sr%ballsize ballsize = sr%ballsize
rearrange = sr%rearrange rearrange = sr%rearrange
ind => sr%ind(1:) ind => sr%ind(1:)
data => sr%Data(1:,1:) myData => sr%myData(1:,1:)
centeridx = sr%centeridx centeridx = sr%centeridx
correltime = sr%correltime correltime = sr%correltime
nn = sr%nn ! number to search for nn = sr%nn ! number to search for
@ -1640,7 +1640,7 @@ contains
if (rearrange) then if (rearrange) then
sd = 0.0_pReal sd = 0.0_pReal
do k = 1_pInt,dimen do k = 1_pInt,dimen
sd = sd + (data(k,i) - qv(k))**2.0_pReal sd = sd + (myData(k,i) - qv(k))**2.0_pReal
if (sd>ballsize) cycle mainloop if (sd>ballsize) cycle mainloop
end do end do
indexofi = ind(i) ! only read it if we have not broken out indexofi = ind(i) ! only read it if we have not broken out
@ -1648,7 +1648,7 @@ contains
indexofi = ind(i) indexofi = ind(i)
sd = 0.0_pReal sd = 0.0_pReal
do k = 1_pInt,dimen do k = 1_pInt,dimen
sd = sd + (data(k,indexofi) - qv(k))**2.0_pReal sd = sd + (myData(k,indexofi) - qv(k))**2.0_pReal
if (sd>ballsize) cycle mainloop if (sd>ballsize) cycle mainloop
end do end do
endif endif
@ -1698,12 +1698,12 @@ contains
do i = 1_pInt, tp%n do i = 1_pInt, tp%n
if (all_distances(i)<results(nn)%dis) then if (all_distances(i)<results(nn)%dis) then
! insert it somewhere on the list ! insert it somewhere on the list
do j = 1, nn do j = 1_pInt, nn
if (all_distances(i)<results(j)%dis) exit if (all_distances(i)<results(j)%dis) exit
end do end do
! now we know 'j' ! now we know 'j'
do k = nn - 1_pInt, j, -1_pInt do k = nn - 1_pInt, j, -1_pInt
results(k+1) = results(k) results(k+1_pInt) = results(k)
end do end do
results(j)%dis = all_distances(i) results(j)%dis = all_distances(i)
results(j)%idx = i results(j)%idx = i
@ -1724,22 +1724,23 @@ contains
integer(pInt), intent(out) :: nfound integer(pInt), intent(out) :: nfound
type(kdtree2_result) :: results(:) type(kdtree2_result) :: results(:)
integer(pInt) :: i, nalloc integer(pInt) :: i
integer :: nalloc
real(pReal), allocatable :: all_distances(:) real(pReal), allocatable :: all_distances(:)
! .. ! ..
allocate (all_distances(tp%n)) allocate (all_distances(tp%n))
do i = 1, tp%n do i = 1_pInt, tp%n
all_distances(i) = square_distance(tp%dimen,qv,tp%the_data(:,i)) all_distances(i) = square_distance(tp%dimen,qv,tp%the_data(:,i))
end do end do
nfound = 0 nfound = 0_pInt
nalloc = size(results,1) nalloc = size(results,1)
do i = 1, tp%n do i = 1_pInt, tp%n
if (all_distances(i)< r2) then if (all_distances(i)< r2) then
! insert it somewhere on the list ! insert it somewhere on the list
if (nfound .lt. nalloc) then if (nfound .lt. nalloc) then
nfound = nfound+1 nfound = nfound+1_pInt
results(nfound)%dis = all_distances(i) results(nfound)%dis = all_distances(i)
results(nfound)%idx = i results(nfound)%idx = i
endif endif
@ -1763,7 +1764,7 @@ contains
!THIS IS BUGGY WITH INTEL FORTRAN !THIS IS BUGGY WITH INTEL FORTRAN
! If (nfound .Gt. 1) Call heapsort(results(1:nfound)%dis,results(1:nfound)%ind,nfound) ! If (nfound .Gt. 1) Call heapsort(results(1:nfound)%dis,results(1:nfound)%ind,nfound)
! !
if (nfound .gt. 1) call heapsort_struct(results,nfound) if (nfound .gt. 1_pInt) call heapsort_struct(results,nfound)
return return
end subroutine kdtree2_sort_results end subroutine kdtree2_sort_results
@ -1786,7 +1787,7 @@ contains
integer(pInt) :: i,j integer(pInt) :: i,j
integer(pInt) :: ileft,iright integer(pInt) :: ileft,iright
ileft=n/2+1 ileft=n/2_pInt+1_pInt
iright=n iright=n
! do i=1,n ! do i=1,n
@ -1794,33 +1795,33 @@ contains
! Generate initial idum array ! Generate initial idum array
! end do ! end do
if(n.eq.1) return if(n.eq.1_pInt) return
do do
if(ileft > 1)then if(ileft > 1_pInt)then
ileft=ileft-1 ileft=ileft-1_pInt
value=a(ileft); ivalue=ind(ileft) value=a(ileft); ivalue=ind(ileft)
else else
value=a(iright); ivalue=ind(iright) value=a(iright); ivalue=ind(iright)
a(iright)=a(1); ind(iright)=ind(1) a(iright)=a(1); ind(iright)=ind(1)
iright=iright-1 iright=iright-1_pInt
if (iright == 1) then if (iright == 1_pInt) then
a(1)=value;ind(1)=ivalue a(1)=value;ind(1)=ivalue
return return
endif endif
endif endif
i=ileft i=ileft
j=2*ileft j=2_pInt*ileft
do while (j <= iright) do while (j <= iright)
if(j < iright) then if(j < iright) then
if(a(j) < a(j+1)) j=j+1 if(a(j) < a(j+1_pInt)) j=j+1_pInt
endif endif
if(value < a(j)) then if(value < a(j)) then
a(i)=a(j); ind(i)=ind(j) a(i)=a(j); ind(i)=ind(j)
i=j i=j
j=j+j j=j+j
else else
j=iright+1 j=iright+1_pInt
endif endif
end do end do
a(i)=value; ind(i)=ivalue a(i)=value; ind(i)=ivalue
@ -1842,7 +1843,7 @@ contains
integer(pInt) :: i,j integer(pInt) :: i,j
integer(pInt) :: ileft,iright integer(pInt) :: ileft,iright
ileft=n/2+1 ileft=n/2_pInt+1_pInt
iright=n iright=n
! do i=1,n ! do i=1,n
@ -1850,33 +1851,33 @@ contains
! Generate initial idum array ! Generate initial idum array
! end do ! end do
if(n.eq.1) return if(n.eq.1_pInt) return
do do
if(ileft > 1)then if(ileft > 1_pInt)then
ileft=ileft-1 ileft=ileft-1_pInt
value=a(ileft) value=a(ileft)
else else
value=a(iright) value=a(iright)
a(iright)=a(1) a(iright)=a(1)
iright=iright-1 iright=iright-1_pInt
if (iright == 1) then if (iright == 1_pInt) then
a(1) = value a(1) = value
return return
endif endif
endif endif
i=ileft i=ileft
j=2*ileft j=2_pInt*ileft
do while (j <= iright) do while (j <= iright)
if(j < iright) then if(j < iright) then
if(a(j)%dis < a(j+1)%dis) j=j+1 if(a(j)%dis < a(j+1_pInt)%dis) j=j+1_pInt
endif endif
if(value%dis < a(j)%dis) then if(value%dis < a(j)%dis) then
a(i)=a(j); a(i)=a(j);
i=j i=j
j=j+j j=j+j
else else
j=iright+1 j=iright+1_pInt
endif endif
end do end do
a(i)=value a(i)=value

42
code/lattice.f90
View File

@ -39,11 +39,11 @@ implicit none
integer(pInt) lattice_Nhexagonal, & ! # of hexagonal lattice structure (from tag CoverA_ratio) integer(pInt) lattice_Nhexagonal, & ! # of hexagonal lattice structure (from tag CoverA_ratio)
lattice_Nstructure ! # of lattice structures (1: fcc,2: bcc,3+: hexagonal) lattice_Nstructure ! # of lattice structures (1: fcc,2: bcc,3+: hexagonal)
integer(pInt), parameter :: lattice_maxNslipFamily = 5 ! max # of slip system families over lattice structures integer(pInt), parameter :: lattice_maxNslipFamily = 5_pInt ! max # of slip system families over lattice structures
integer(pInt), parameter :: lattice_maxNtwinFamily = 4 ! max # of twin system families over lattice structures integer(pInt), parameter :: lattice_maxNtwinFamily = 4_pInt ! max # of twin system families over lattice structures
integer(pInt), parameter :: lattice_maxNslip = 54 ! max # of slip systems over lattice structures integer(pInt), parameter :: lattice_maxNslip = 54_pInt ! max # of slip systems over lattice structures
integer(pInt), parameter :: lattice_maxNtwin = 24 ! max # of twin systems over lattice structures integer(pInt), parameter :: lattice_maxNtwin = 24_pInt ! max # of twin systems over lattice structures
integer(pInt), parameter :: lattice_maxNinteraction = 30 ! max # of interaction types (in hardening matrix part) integer(pInt), parameter :: lattice_maxNinteraction = 30_pInt ! max # of interaction types (in hardening matrix part)
integer(pInt), pointer, dimension(:,:) :: interactionSlipSlip, & integer(pInt), pointer, dimension(:,:) :: interactionSlipSlip, &
interactionSlipTwin, & interactionSlipTwin, &
@ -81,10 +81,10 @@ integer(pInt), allocatable, dimension(:,:,:) :: lattice_interactionSlipSlip, &
!============================== fcc (1) ================================= !============================== fcc (1) =================================
integer(pInt), parameter, dimension(lattice_maxNslipFamily) :: lattice_fcc_NslipSystem = (/12, 0, 0, 0, 0/) integer(pInt), parameter, dimension(lattice_maxNslipFamily) :: lattice_fcc_NslipSystem = int([12, 0, 0, 0, 0],pInt)
integer(pInt), parameter, dimension(lattice_maxNtwinFamily) :: lattice_fcc_NtwinSystem = (/12, 0, 0, 0/) integer(pInt), parameter, dimension(lattice_maxNtwinFamily) :: lattice_fcc_NtwinSystem = int([12, 0, 0, 0],pInt)
integer(pInt), parameter :: lattice_fcc_Nslip = 12 ! sum(lattice_fcc_NslipSystem) integer(pInt), parameter :: lattice_fcc_Nslip = 12_pInt ! sum(lattice_fcc_NslipSystem)
integer(pInt), parameter :: lattice_fcc_Ntwin = 12 ! sum(lattice_fcc_NtwinSystem) integer(pInt), parameter :: lattice_fcc_Ntwin = 12_pInt ! sum(lattice_fcc_NtwinSystem)
integer(pInt) :: lattice_fcc_Nstructure = 0_pInt integer(pInt) :: lattice_fcc_Nstructure = 0_pInt
real(pReal), dimension(3+3,lattice_fcc_Nslip), parameter :: lattice_fcc_systemSlip = & real(pReal), dimension(3+3,lattice_fcc_Nslip), parameter :: lattice_fcc_systemSlip = &
@ -442,10 +442,10 @@ integer(pInt), allocatable, dimension(:,:,:) :: lattice_interactionSlipSlip, &
!============================== hex (3+) ================================= !============================== hex (3+) =================================
integer(pInt), parameter, dimension(lattice_maxNslipFamily) :: lattice_hex_NslipSystem = (/ 3, 3, 6,12, 6/) integer(pInt), parameter, dimension(lattice_maxNslipFamily) :: lattice_hex_NslipSystem = int([ 3, 3, 6,12, 6],pInt)
integer(pInt), parameter, dimension(lattice_maxNtwinFamily) :: lattice_hex_NtwinSystem = (/ 6, 6, 6, 6/) integer(pInt), parameter, dimension(lattice_maxNtwinFamily) :: lattice_hex_NtwinSystem = int([ 6, 6, 6, 6],pInt)
integer(pInt), parameter :: lattice_hex_Nslip = 30 ! sum(lattice_hex_NslipSystem) integer(pInt), parameter :: lattice_hex_Nslip = 30_pInt ! sum(lattice_hex_NslipSystem)
integer(pInt), parameter :: lattice_hex_Ntwin = 24 ! sum(lattice_hex_NtwinSystem) integer(pInt), parameter :: lattice_hex_Ntwin = 24_pInt ! sum(lattice_hex_NtwinSystem)
integer(pInt) :: lattice_hex_Nstructure = 0_pInt integer(pInt) :: lattice_hex_Nstructure = 0_pInt
!* sorted by A. Alankar & P. Eisenlohr !* sorted by A. Alankar & P. Eisenlohr
@ -824,12 +824,12 @@ function lattice_initializeStructure(struct,CoverA)
lattice_fcc_Nstructure = lattice_fcc_Nstructure + 1_pInt ! count fcc instances lattice_fcc_Nstructure = lattice_fcc_Nstructure + 1_pInt ! count fcc instances
if (lattice_fcc_Nstructure == 1_pInt) then ! me is first fcc structure if (lattice_fcc_Nstructure == 1_pInt) then ! me is first fcc structure
processMe = .true. processMe = .true.
do i = 1,myNslip ! calculate slip system vectors do i = 1_pInt,myNslip ! calculate slip system vectors
sd(1:3,i) = lattice_fcc_systemSlip(1:3,i)/sqrt(math_mul3x3(lattice_fcc_systemSlip(1:3,i),lattice_fcc_systemSlip(1:3,i))) sd(1:3,i) = lattice_fcc_systemSlip(1:3,i)/sqrt(math_mul3x3(lattice_fcc_systemSlip(1:3,i),lattice_fcc_systemSlip(1:3,i)))
sn(1:3,i) = lattice_fcc_systemSlip(4:6,i)/sqrt(math_mul3x3(lattice_fcc_systemSlip(4:6,i),lattice_fcc_systemSlip(4:6,i))) sn(1:3,i) = lattice_fcc_systemSlip(4:6,i)/sqrt(math_mul3x3(lattice_fcc_systemSlip(4:6,i),lattice_fcc_systemSlip(4:6,i)))
st(1:3,i) = math_vectorproduct(sd(1:3,i),sn(1:3,i)) st(1:3,i) = math_vectorproduct(sd(1:3,i),sn(1:3,i))
enddo enddo
do i = 1,myNtwin ! calculate twin system vectors and (assign) shears do i = 1_pInt,myNtwin ! calculate twin system vectors and (assign) shears
td(1:3,i) = lattice_fcc_systemTwin(1:3,i)/sqrt(math_mul3x3(lattice_fcc_systemTwin(1:3,i),lattice_fcc_systemTwin(1:3,i))) td(1:3,i) = lattice_fcc_systemTwin(1:3,i)/sqrt(math_mul3x3(lattice_fcc_systemTwin(1:3,i),lattice_fcc_systemTwin(1:3,i)))
tn(1:3,i) = lattice_fcc_systemTwin(4:6,i)/sqrt(math_mul3x3(lattice_fcc_systemTwin(4:6,i),lattice_fcc_systemTwin(4:6,i))) tn(1:3,i) = lattice_fcc_systemTwin(4:6,i)/sqrt(math_mul3x3(lattice_fcc_systemTwin(4:6,i),lattice_fcc_systemTwin(4:6,i)))
tt(1:3,i) = math_vectorproduct(td(1:3,i),tn(1:3,i)) tt(1:3,i) = math_vectorproduct(td(1:3,i),tn(1:3,i))
@ -850,12 +850,12 @@ function lattice_initializeStructure(struct,CoverA)
lattice_bcc_Nstructure = lattice_bcc_Nstructure + 1_pInt ! count bcc instances lattice_bcc_Nstructure = lattice_bcc_Nstructure + 1_pInt ! count bcc instances
if (lattice_bcc_Nstructure == 1_pInt) then ! me is first bcc structure if (lattice_bcc_Nstructure == 1_pInt) then ! me is first bcc structure
processMe = .true. processMe = .true.
do i = 1,myNslip ! calculate slip system vectors do i = 1_pInt,myNslip ! calculate slip system vectors
sd(1:3,i) = lattice_bcc_systemSlip(1:3,i)/sqrt(math_mul3x3(lattice_bcc_systemSlip(1:3,i),lattice_bcc_systemSlip(1:3,i))) sd(1:3,i) = lattice_bcc_systemSlip(1:3,i)/sqrt(math_mul3x3(lattice_bcc_systemSlip(1:3,i),lattice_bcc_systemSlip(1:3,i)))
sn(1:3,i) = lattice_bcc_systemSlip(4:6,i)/sqrt(math_mul3x3(lattice_bcc_systemSlip(4:6,i),lattice_bcc_systemSlip(4:6,i))) sn(1:3,i) = lattice_bcc_systemSlip(4:6,i)/sqrt(math_mul3x3(lattice_bcc_systemSlip(4:6,i),lattice_bcc_systemSlip(4:6,i)))
st(1:3,i) = math_vectorproduct(sd(1:3,i),sn(1:3,i)) st(1:3,i) = math_vectorproduct(sd(1:3,i),sn(1:3,i))
enddo enddo
do i = 1,myNtwin ! calculate twin system vectors and (assign) shears do i = 1_pInt,myNtwin ! calculate twin system vectors and (assign) shears
td(1:3,i) = lattice_bcc_systemTwin(1:3,i)/sqrt(math_mul3x3(lattice_bcc_systemTwin(1:3,i),lattice_bcc_systemTwin(1:3,i))) td(1:3,i) = lattice_bcc_systemTwin(1:3,i)/sqrt(math_mul3x3(lattice_bcc_systemTwin(1:3,i),lattice_bcc_systemTwin(1:3,i)))
tn(1:3,i) = lattice_bcc_systemTwin(4:6,i)/sqrt(math_mul3x3(lattice_bcc_systemTwin(4:6,i),lattice_bcc_systemTwin(4:6,i))) tn(1:3,i) = lattice_bcc_systemTwin(4:6,i)/sqrt(math_mul3x3(lattice_bcc_systemTwin(4:6,i),lattice_bcc_systemTwin(4:6,i)))
tt(1:3,i) = math_vectorproduct(td(1:3,i),tn(1:3,i)) tt(1:3,i) = math_vectorproduct(td(1:3,i),tn(1:3,i))
@ -877,7 +877,7 @@ function lattice_initializeStructure(struct,CoverA)
myNtwin = lattice_hex_Ntwin ! overall number of twin systems myNtwin = lattice_hex_Ntwin ! overall number of twin systems
processMe = .true. processMe = .true.
! converting from 4 axes coordinate system (a1=a2=a3=c) to ortho-hexgonal system (a, b, c) ! converting from 4 axes coordinate system (a1=a2=a3=c) to ortho-hexgonal system (a, b, c)
do i = 1,myNslip do i = 1_pInt,myNslip
hex_d(1) = lattice_hex_systemSlip(1,i)*1.5_pReal ! direction [uvtw]->[3u/2 (u+2v)*sqrt(3)/2 w*(c/a)] hex_d(1) = lattice_hex_systemSlip(1,i)*1.5_pReal ! direction [uvtw]->[3u/2 (u+2v)*sqrt(3)/2 w*(c/a)]
hex_d(2) = (lattice_hex_systemSlip(1,i)+2.0_pReal*lattice_hex_systemSlip(2,i))*(0.5_pReal*sqrt(3.0_pReal)) hex_d(2) = (lattice_hex_systemSlip(1,i)+2.0_pReal*lattice_hex_systemSlip(2,i))*(0.5_pReal*sqrt(3.0_pReal))
hex_d(3) = lattice_hex_systemSlip(4,i)*CoverA hex_d(3) = lattice_hex_systemSlip(4,i)*CoverA
@ -889,7 +889,7 @@ function lattice_initializeStructure(struct,CoverA)
sn(1:3,i) = hex_n/sqrt(math_mul3x3(hex_n,hex_n)) sn(1:3,i) = hex_n/sqrt(math_mul3x3(hex_n,hex_n))
st(1:3,i) = math_vectorproduct(sd(1:3,i),sn(1:3,i)) st(1:3,i) = math_vectorproduct(sd(1:3,i),sn(1:3,i))
enddo enddo
do i = 1,myNtwin do i = 1_pInt,myNtwin
hex_d(1) = lattice_hex_systemTwin(1,i)*1.5_pReal hex_d(1) = lattice_hex_systemTwin(1,i)*1.5_pReal
hex_d(2) = (lattice_hex_systemTwin(1,i)+2.0_pReal*lattice_hex_systemTwin(2,i))*(0.5_pReal*sqrt(3.0_pReal)) hex_d(2) = (lattice_hex_systemTwin(1,i)+2.0_pReal*lattice_hex_systemTwin(2,i))*(0.5_pReal*sqrt(3.0_pReal))
hex_d(3) = lattice_hex_systemTwin(4,i)*CoverA hex_d(3) = lattice_hex_systemTwin(4,i)*CoverA
@ -923,14 +923,14 @@ function lattice_initializeStructure(struct,CoverA)
if (processMe) then if (processMe) then
if (myStructure > lattice_Nstructure) & if (myStructure > lattice_Nstructure) &
call IO_error(666_pInt,0_pInt,0_pInt,0_pInt,'structure index too large') ! check for memory leakage call IO_error(666_pInt,0_pInt,0_pInt,0_pInt,'structure index too large') ! check for memory leakage
do i = 1,myNslip ! store slip system vectors and Schmid matrix for my structure do i = 1_pInt,myNslip ! store slip system vectors and Schmid matrix for my structure
lattice_sd(1:3,i,myStructure) = sd(1:3,i) lattice_sd(1:3,i,myStructure) = sd(1:3,i)
lattice_st(1:3,i,myStructure) = st(1:3,i) lattice_st(1:3,i,myStructure) = st(1:3,i)
lattice_sn(1:3,i,myStructure) = sn(1:3,i) lattice_sn(1:3,i,myStructure) = sn(1:3,i)
lattice_Sslip(1:3,1:3,i,myStructure) = math_tensorproduct(sd(1:3,i),sn(1:3,i)) lattice_Sslip(1:3,1:3,i,myStructure) = math_tensorproduct(sd(1:3,i),sn(1:3,i))
lattice_Sslip_v(1:6,i,myStructure) = math_Mandel33to6(math_symmetric33(lattice_Sslip(1:3,1:3,i,myStructure))) lattice_Sslip_v(1:6,i,myStructure) = math_Mandel33to6(math_symmetric33(lattice_Sslip(1:3,1:3,i,myStructure)))
enddo enddo
do i = 1,myNtwin ! store twin system vectors and Schmid plus rotation matrix for my structure do i = 1_pInt,myNtwin ! store twin system vectors and Schmid plus rotation matrix for my structure
lattice_td(1:3,i,myStructure) = td(1:3,i) lattice_td(1:3,i,myStructure) = td(1:3,i)
lattice_tt(1:3,i,myStructure) = tt(1:3,i) lattice_tt(1:3,i,myStructure) = tt(1:3,i)
lattice_tn(1:3,i,myStructure) = tn(1:3,i) lattice_tn(1:3,i,myStructure) = tn(1:3,i)

151
code/makefile
View File

@ -32,42 +32,13 @@
# PREFIX = arbitrary prefix # PREFIX = arbitrary prefix
# SUFFIX = arbitrary suffix # SUFFIX = arbitrary suffix
# STANDARD_CHECK = checking for Fortran 2008, compiler dependend # STANDARD_CHECK = checking for Fortran 2008, compiler dependend
########################################################################################
# Here are some useful debugging switches for ifort. Switch on by uncommenting the #SUFFIX line at the end of this section:
# 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!
DEBUG1 =-check bounds -g
#will cause a lot of warnings because we create a bunch of temporary arrays
DEBUG2 =-check arg_temp_created
#check from time to time
DEBUG3 =-fp-stack-check -g -traceback -gen-interfaces -warn interfaces
#should not be done for OpenMP, but set "ulimit -s unlimited" on shell. Probably it helps also to unlimit other limits
DEBUG4 =-heap-arrays
#additional warnings
DEBUG5 =-warn all
# or one or more of those: alignments, declarations,general, ignore_loc, uncalled, unuses, usage
#set precision (check for missing _pInt and _pReal)
DEBUG6= -real-size 32 -integer-size 16
#or one of those 16/32/64/128 (= 2,4,8,16 bytes)
#SUFFIX =$(DEBUG1) $(DEBUG2) $(DEBUG3) $(DEBUG4) $(DEBUG5) $(DEBUG6)
# Here are some useful debugging switches for gfortran
# fcheck-bounds: eqv to DEBUG1 of ifort
######################################################################################## ########################################################################################
#auto values will be set by setup_code.py #auto values will be set by setup_code.py
FFTWROOT := $(DAMASK_ROOT)/lib/fftw FFTWROOT :=/$(DAMASK_ROOT)/lib/fftw
IKMLROOT := IKMLROOT :=
ACMLROOT := /opt/acml4.4.0 ACMLROOT :=/opt/acml4.4.0
LAPACKROOT := #LAPACKROOT := /usr
F90 ?= ifort F90 ?= ifort
COMPILERNAME ?= $(F90) COMPILERNAME ?= $(F90)
@ -144,38 +115,71 @@ OPTIMIZATION_OFF_ifort :=-O0
OPTIMIZATION_OFF_gfortran :=-O0 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) -ipo -static -no-prec-div -fp-model fast=2
OPTIMIZATION_AGGRESSIVE_gfortran :=-O3 $(PORTABLE_SWITCH) -ffast-math -funroll-loops -ftree-vectorize OPTIMIZATION_AGGRESSIVE_gfortran :=-O3 $(PORTABLE_SWITCH) -ffast-math -funroll-loops -ftree-vectorize
COMPILE_OPTIONS_ifort :=-fpp\ COMPILE_OPTIONS_ifort :=-fpp\
-implicitnone\
-diag-enable sc3\ -diag-enable sc3\
-diag-disable 8291,8290,5268\ -diag-disable 5268\
-warn declarations\ -warn declarations\
-warn general\ -warn general\
-warn usage -warn usage\
-warn interfaces\
#alignments: Determines whether warnings occur for data that is not naturally aligned. -warn ignore_loc\
#declarations: Determines whether warnings occur for any undeclared names. -warn alignments\
#errors: Determines whether warnings are changed to errors. -warn unused\
#general: Determines whether warning messages and informational messages are issued by the compiler. -warn errors\
#ignore_loc: Determines whether warnings occur when %LOC is stripped from an actual argument. -warn stderrors
#interfaces: Determines whether the compiler checks the interfaces of all SUBROUTINEs called and FUNCTIONs invoked in your compilation against an external set of interface blocks.
#stderrors: Determines whether warnings about Fortran standard violations are changed to errors.
#truncated_source: Determines whether warnings occur when source exceeds the maximum column width in fixed-format files.
#uncalled: Determines whether warnings occur when a statement function is never called
#unused: Determines whether warnings occur for declared variables that are never used.
#usage: Determines whether warnings occur for questionable programming practices.
#-fpp: preprocessor #-fpp: preprocessor
#-fimplicit-none: assume "implicit-none" even if not present in source
#-diag-disable: disables warnings, where #-diag-disable: disables warnings, where
# warning ID 9291: # warning ID 5268: the text exceeds right hand column allowed on the line (we have only comments there)
# warning ID 8290: #-warn: enables warnings, where
# warning ID 5268: The text exceeds right hand column allowed on the line (we have only comments there) # declarations: any undeclared names
# general: warning messages and informational messages are issued by the compiler
# usage: questionable programming practices
# interfaces: checks the interfaces of all SUBROUTINEs called and FUNCTIONs invoked in your compilation against an external set of interface blocks
# ignore_loc: %LOC is stripped from an actual argument
# alignments: data that is not naturally aligned
# unused: declared variables that are never used
# errors: warnings are changed to errors
# stderrors: warnings about Fortran standard violations are changed to errors
#
###################################################################################################
#MORE OPTIONS FOR DEBUGGING DURING COMPILING
#-warn: enables warnings, where
# truncated_source: Determines whether warnings occur when source exceeds the maximum column width in fixed-format files. (too many warnings because we have comments beyond character 132)
# uncalled: Determines whether warnings occur when a statement function is never called
# all:
#
#OPTIONS FOR DEGUBBING DURING RUNTIME
# information on http://software.intel.com/en-us/articles/determining-root-cause-of-sigsegv-or-sigbus-errors/
#-g: Generate symbolic debugging information in the object file
#-traceback: Generate extra information in the object file to provide source file traceback information when a severe error occurs at run time.
#-gen-interfaces: Generate an interface block for each routine. http://software.intel.com/en-us/blogs/2012/01/05/doctor-fortran-gets-explicit-again/
#-fp-stack-check: Generate extra code after every function call to ensure that the floating-point (FP) stack is in the expected state.
#-check: checks at runtime, where
# bounds: check if an array index is too small (<1) or too large!
# arg_temp_created: will cause a lot of warnings because we create a bunch of temporary arrays
# format: Checking for the data type of an item being formatted for output.
# output_conversion: Checking for the fit of data items within a designated format descriptor field.
# pointers: Checking for certain disassociated or uninitialized pointers or unallocated allocatable objects.
# uninit: Checking for uninitialized variables.
#-heap-arrays: should not be done for OpenMP, but set "ulimit -s unlimited" on shell. Probably it helps also to unlimit other limits
#
#OPTIONS FOR TYPE DEBUGGING
#-real-size 32: set precision to one of those 32/64/128 (= 4/8/16 bytes) for standard real (=8 for pReal)
#-integer-size 16: set precision to one of those 16/32/64 (= 2/4/8 bytes) for standard integer (=4 for pInt)
###################################################################################################
COMPILE_OPTIONS_gfortran :=-xf95-cpp-input\ COMPILE_OPTIONS_gfortran :=-xf95-cpp-input\
-ffree-line-length-132\ -ffree-line-length-132\
-fno-range-check\ -fno-range-check\
-fimplicit-none\ -fimplicit-none\
-fall-intrinsics\
-pedantic\ -pedantic\
-Warray-bounds\ -Warray-bounds\
-Wunused-parameter\ -Wunused-parameter\
@ -193,6 +197,7 @@ COMPILE_OPTIONS_gfortran :=-xf95-cpp-input\
#-ffree-line-length-132: restrict line length to the standard 132 characters #-ffree-line-length-132: restrict line length to the standard 132 characters
#-fno-range-check: disables checking if result can be represented by variable. Needs to be set to enable DAMASK_NaN #-fno-range-check: disables checking if result can be represented by variable. Needs to be set to enable DAMASK_NaN
#-fimplicit-none: assume "implicit-none" even if not present in source #-fimplicit-none: assume "implicit-none" even if not present in source
#-fall-intrinsics:
#-pedantic: more strict on standard, enables some of the warnings below #-pedantic: more strict on standard, enables some of the warnings below
#-Warray-bounds: checks if array reference is out of bounds at compile time. use -fcheck-bounds to also check during runtime #-Warray-bounds: checks if array reference is out of bounds at compile time. use -fcheck-bounds to also check during runtime
#-Wunused-parameter: find usused variables with "parameter" attribute #-Wunused-parameter: find usused variables with "parameter" attribute
@ -205,31 +210,33 @@ COMPILE_OPTIONS_gfortran :=-xf95-cpp-input\
#-Wsurprising: warn when "suspicious" code constructs are encountered. While technically legal these usually indicate that an error has been made. #-Wsurprising: warn when "suspicious" code constructs are encountered. While technically legal these usually indicate that an error has been made.
#-Wunused-value: #-Wunused-value:
#-Wunderflow: produce a warning when numerical constant expressions are encountered, which yield an UNDERFLOW during compilation #-Wunderflow: produce a warning when numerical constant expressions are encountered, which yield an UNDERFLOW during compilation
#
###################################################################################################
#OPTIONS FOR GFORTRAN 4.6
#MORE OPTIONS #-Wsuggest-attribute=const:
#-Wsuggest-attribute=noreturn:
# only for gfortran 4.6: #-Wsuggest-attribute=pure:
#-Wsuggest-attribute=const #
#-Wsuggest-attribute=noreturn #MORE OPTIONS FOR DEBUGGING DURING COMPILING
#-Wsuggest-attribute=pure #-Wline-truncation: too many warnings because we have comments beyond character 132
#-Wintrinsic-std: warnings because of "flush" is not longer in the standard, but still an intrinsic fuction of the compilers:
# too many warnings because we have comments beyond character 132: #-Warray-temporarieswarnings:
#-Wline-truncation # because we have many temporary arrays (performance issue?):
# warnings because of "flush" is not longer in the standard, but still an intrinsic fuction of the compilers: #-Wimplicit-interface
#-Wintrinsic-std #-pedantic-errors
# warnings because we have many temporary arrays (performance issue?): #-fmodule-private
#-Warray-temporaries #
#OPTIONS FOR DEGUBBING DURING RUNTIME
# -Wimplicit-interface #-fcheck-bounds: check if an array index is too small (<1) or too large!
# -pedantic-errors #
# -fmodule-private #OPTIONS FOR TYPE DEBUGGING
#-fdefault-real-8: set precision to 8 bytes for standard real (=8 for pReal). Will set size of double to 16 bytes as long as -fdefault-double-8 is not set
#-fdefault-integer-8: set precision to 8 bytes for standard integer (=4 for pInt)
##################################################################################################
COMPILE =$(OPENMP_FLAG_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(OPTI)_$(F90)) -c COMPILE =$(OPENMP_FLAG_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(OPTI)_$(F90)) -c
COMPILE_MAXOPTI =$(OPENMP_FLAG_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) -c COMPILE_MAXOPTI =$(OPENMP_FLAG_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) -c
###################################################################################################
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 \ $(PREFIX) $(COMPILERNAME) ${OPENMP_FLAG_${F90}} -o DAMASK_spectral.exe DAMASK_spectral.o CPFEM.a \

242
code/material.f90
View File

@ -114,7 +114,7 @@ subroutine material_init()
implicit none implicit none
!* Definition of variables !* Definition of variables
integer(pInt), parameter :: fileunit = 200 integer(pInt), parameter :: fileunit = 200_pInt
integer(pInt) i,j integer(pInt) i,j
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
@ -128,46 +128,46 @@ subroutine material_init()
call IO_open_file(fileunit,material_configFile) ! ...open material.config file call IO_open_file(fileunit,material_configFile) ! ...open material.config file
endif endif
call material_parseHomogenization(fileunit,material_partHomogenization) call material_parseHomogenization(fileunit,material_partHomogenization)
if (debug_verbosity > 0) then if (debug_verbosity > 0_pInt) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write (6,*) 'Homogenization parsed' write (6,*) 'Homogenization parsed'
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
call material_parseMicrostructure(fileunit,material_partMicrostructure) call material_parseMicrostructure(fileunit,material_partMicrostructure)
if (debug_verbosity > 0) then if (debug_verbosity > 0_pInt) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write (6,*) 'Microstructure parsed' write (6,*) 'Microstructure parsed'
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
call material_parseCrystallite(fileunit,material_partCrystallite) call material_parseCrystallite(fileunit,material_partCrystallite)
if (debug_verbosity > 0) then if (debug_verbosity > 0_pInt) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write (6,*) 'Crystallite parsed' write (6,*) 'Crystallite parsed'
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
call material_parseTexture(fileunit,material_partTexture) call material_parseTexture(fileunit,material_partTexture)
if (debug_verbosity > 0) then if (debug_verbosity > 0_pInt) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write (6,*) 'Texture parsed' write (6,*) 'Texture parsed'
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
call material_parsePhase(fileunit,material_partPhase) call material_parsePhase(fileunit,material_partPhase)
if (debug_verbosity > 0) then if (debug_verbosity > 0_pInt) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write (6,*) 'Phase parsed' write (6,*) 'Phase parsed'
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
close(fileunit) close(fileunit)
do i = 1,material_Nmicrostructure do i = 1_pInt,material_Nmicrostructure
if (microstructure_crystallite(i) < 1 .or. & if (microstructure_crystallite(i) < 1_pInt .or. &
microstructure_crystallite(i) > material_Ncrystallite) call IO_error(150_pInt,i) microstructure_crystallite(i) > material_Ncrystallite) call IO_error(150_pInt,i)
if (minval(microstructure_phase(1:microstructure_Nconstituents(i),i)) < 1 .or. & if (minval(microstructure_phase(1:microstructure_Nconstituents(i),i)) < 1_pInt .or. &
maxval(microstructure_phase(1:microstructure_Nconstituents(i),i)) > material_Nphase) call IO_error(151_pInt,i) maxval(microstructure_phase(1:microstructure_Nconstituents(i),i)) > material_Nphase) call IO_error(151_pInt,i)
if (minval(microstructure_texture(1:microstructure_Nconstituents(i),i)) < 1 .or. & if (minval(microstructure_texture(1:microstructure_Nconstituents(i),i)) < 1_pInt .or. &
maxval(microstructure_texture(1:microstructure_Nconstituents(i),i)) > material_Ntexture) call IO_error(152_pInt,i) maxval(microstructure_texture(1:microstructure_Nconstituents(i),i)) > material_Ntexture) call IO_error(152_pInt,i)
if (abs(sum(microstructure_fraction(:,i)) - 1.0_pReal) >= 1.0e-10_pReal) then if (abs(sum(microstructure_fraction(:,i)) - 1.0_pReal) >= 1.0e-10_pReal) then
if (debug_verbosity > 0) then if (debug_verbosity > 0_pInt) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*)'sum of microstructure fraction = ',sum(microstructure_fraction(:,i)) write(6,*)'sum of microstructure fraction = ',sum(microstructure_fraction(:,i))
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
@ -175,25 +175,25 @@ subroutine material_init()
call IO_error(153_pInt,i) call IO_error(153_pInt,i)
endif endif
enddo enddo
if (debug_verbosity > 0) then if (debug_verbosity > 0_pInt) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write (6,*) write (6,*)
write (6,*) 'MATERIAL configuration' write (6,*) 'MATERIAL configuration'
write (6,*) write (6,*)
write (6,'(a32,1x,a16,1x,a6)') 'homogenization ','type ','grains' write (6,'(a32,1x,a16,1x,a6)') 'homogenization ','type ','grains'
do i = 1,material_Nhomogenization do i = 1_pInt,material_Nhomogenization
write (6,'(1x,a32,1x,a16,1x,i4)') homogenization_name(i),homogenization_type(i),homogenization_Ngrains(i) write (6,'(1x,a32,1x,a16,1x,i4)') homogenization_name(i),homogenization_type(i),homogenization_Ngrains(i)
enddo enddo
write (6,*) write (6,*)
write (6,'(a32,1x,a11,1x,a12,1x,a13)') 'microstructure ','crystallite','constituents','homogeneous' write (6,'(a32,1x,a11,1x,a12,1x,a13)') 'microstructure ','crystallite','constituents','homogeneous'
do i = 1,material_Nmicrostructure do i = 1_pInt,material_Nmicrostructure
write (6,'(a32,4x,i4,8x,i4,8x,l1)') microstructure_name(i), & write (6,'(a32,4x,i4,8x,i4,8x,l1)') microstructure_name(i), &
microstructure_crystallite(i), & microstructure_crystallite(i), &
microstructure_Nconstituents(i), & microstructure_Nconstituents(i), &
microstructure_elemhomo(i) microstructure_elemhomo(i)
if (microstructure_Nconstituents(i) > 0_pInt) then if (microstructure_Nconstituents(i) > 0_pInt) then
do j = 1,microstructure_Nconstituents(i) do j = 1_pInt,microstructure_Nconstituents(i)
write (6,'(a1,1x,a32,1x,a32,1x,f6.4)') '>',phase_name(microstructure_phase(j,i)),& write (6,'(a1,1x,a32,1x,a32,1x,f7.4)') '>',phase_name(microstructure_phase(j,i)),&
texture_name(microstructure_texture(j,i)),& texture_name(microstructure_texture(j,i)),&
microstructure_fraction(j,i) microstructure_fraction(j,i)
enddo enddo
@ -209,7 +209,7 @@ endsubroutine
!********************************************************************* !*********************************************************************
subroutine material_parseHomogenization(file,myPart) subroutine material_parseHomogenization(myFile,myPart)
!********************************************************************* !*********************************************************************
use prec, only: pInt use prec, only: pInt
@ -218,14 +218,14 @@ subroutine material_parseHomogenization(file,myPart)
implicit none implicit none
character(len=*), intent(in) :: myPart character(len=*), intent(in) :: myPart
integer(pInt), intent(in) :: file integer(pInt), intent(in) :: myFile
integer(pInt), parameter :: maxNchunks = 2_pInt integer(pInt), parameter :: maxNchunks = 2_pInt
integer(pInt), dimension(1+2*maxNchunks) :: positions integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt) Nsections, section, s integer(pInt) Nsections, section, s
character(len=64) tag character(len=64) tag
character(len=1024) line character(len=1024) line
Nsections = IO_countSections(file,myPart) Nsections = IO_countSections(myFile,myPart)
material_Nhomogenization = Nsections material_Nhomogenization = Nsections
if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart) if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart)
@ -236,23 +236,23 @@ subroutine material_parseHomogenization(file,myPart)
allocate(homogenization_Noutput(Nsections)); homogenization_Noutput = 0_pInt allocate(homogenization_Noutput(Nsections)); homogenization_Noutput = 0_pInt
allocate(homogenization_active(Nsections)); homogenization_active = .false. allocate(homogenization_active(Nsections)); homogenization_active = .false.
forall (s = 1:Nsections) homogenization_active(s) = any(mesh_element(3,:) == s) ! current homogenization used in model? Homogenization view, maximum operations depend on maximum number of homog schemes forall (s = 1_pInt:Nsections) homogenization_active(s) = any(mesh_element(3,:) == s) ! current homogenization used in model? Homogenization view, maximum operations depend on maximum number of homog schemes
homogenization_Noutput = IO_countTagInPart(file,myPart,'(output)',Nsections) homogenization_Noutput = IO_countTagInPart(myFile,myPart,'(output)',Nsections)
rewind(file) rewind(myFile)
line = '' line = ''
section = 0 section = 0_pInt
do while (IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to myPart do while (IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to myPart
read(file,'(a1024)',END=100) line read(myFile,'(a1024)',END=100) line
enddo enddo
do do
read(file,'(a1024)',END=100) line read(myFile,'(a1024)',END=100) line
if (IO_isBlank(line)) cycle ! skip empty lines if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') then ! next section if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1 section = section + 1_pInt
homogenization_name(section) = IO_getTag(line,'[',']') homogenization_name(section) = IO_getTag(line,'[',']')
endif endif
if (section > 0_pInt) then if (section > 0_pInt) then
@ -277,7 +277,7 @@ subroutine material_parseHomogenization(file,myPart)
!********************************************************************* !*********************************************************************
subroutine material_parseMicrostructure(file,myPart) subroutine material_parseMicrostructure(myFile,myPart)
!********************************************************************* !*********************************************************************
use prec, only: pInt use prec, only: pInt
@ -286,14 +286,14 @@ subroutine material_parseMicrostructure(file,myPart)
implicit none implicit none
character(len=*), intent(in) :: myPart character(len=*), intent(in) :: myPart
integer(pInt), intent(in) :: file integer(pInt), intent(in) :: myFile
integer(pInt), parameter :: maxNchunks = 7_pInt integer(pInt), parameter :: maxNchunks = 7_pInt
integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions
integer(pInt) Nsections, section, constituent, e, i integer(pInt) Nsections, section, constituent, e, i
character(len=64) tag character(len=64) tag
character(len=1024) line character(len=1024) line
Nsections = IO_countSections(file,myPart) Nsections = IO_countSections(myFile,myPart)
material_Nmicrostructure = Nsections material_Nmicrostructure = Nsections
if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart) if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart)
@ -303,42 +303,42 @@ subroutine material_parseMicrostructure(file,myPart)
allocate(microstructure_active(Nsections)) allocate(microstructure_active(Nsections))
allocate(microstructure_elemhomo(Nsections)) allocate(microstructure_elemhomo(Nsections))
forall (e = 1:mesh_NcpElems) microstructure_active(mesh_element(4,e)) = .true. ! current microstructure used in model? Elementwise view, maximum N operations for N elements forall (e = 1_pInt:mesh_NcpElems) microstructure_active(mesh_element(4,e)) = .true. ! current microstructure used in model? Elementwise view, maximum N operations for N elements
microstructure_Nconstituents = IO_countTagInPart(file,myPart,'(constituent)',Nsections) microstructure_Nconstituents = IO_countTagInPart(myFile,myPart,'(constituent)',Nsections)
microstructure_maxNconstituents = maxval(microstructure_Nconstituents) microstructure_maxNconstituents = maxval(microstructure_Nconstituents)
microstructure_elemhomo = IO_spotTagInPart(file,myPart,'/elementhomogeneous/',Nsections) microstructure_elemhomo = IO_spotTagInPart(myFile,myPart,'/elementhomogeneous/',Nsections)
allocate(microstructure_phase (microstructure_maxNconstituents,Nsections)); microstructure_phase = 0_pInt allocate(microstructure_phase (microstructure_maxNconstituents,Nsections)); microstructure_phase = 0_pInt
allocate(microstructure_texture (microstructure_maxNconstituents,Nsections)); microstructure_texture = 0_pInt allocate(microstructure_texture (microstructure_maxNconstituents,Nsections)); microstructure_texture = 0_pInt
allocate(microstructure_fraction(microstructure_maxNconstituents,Nsections)); microstructure_fraction = 0.0_pReal allocate(microstructure_fraction(microstructure_maxNconstituents,Nsections)); microstructure_fraction = 0.0_pReal
rewind(file) rewind(myFile)
line = '' line = ''
section = 0 section = 0_pInt
do while (IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to myPart do while (IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to myPart
read(file,'(a1024)',END=100) line read(myFile,'(a1024)',END=100) line
enddo enddo
do do
read(file,'(a1024)',END=100) line read(myFile,'(a1024)',END=100) line
if (IO_isBlank(line)) cycle ! skip empty lines if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') then ! next section if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1 section = section + 1_pInt
constituent = 0 constituent = 0_pInt
microstructure_name(section) = IO_getTag(line,'[',']') microstructure_name(section) = IO_getTag(line,'[',']')
endif endif
if (section > 0) then if (section > 0_pInt) then
positions = IO_stringPos(line,maxNchunks) positions = IO_stringPos(line,maxNchunks)
tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key
select case(tag) select case(tag)
case ('crystallite') case ('crystallite')
microstructure_crystallite(section) = IO_intValue(line,positions,2_pInt) microstructure_crystallite(section) = IO_intValue(line,positions,2_pInt)
case ('(constituent)') case ('(constituent)')
constituent = constituent + 1 constituent = constituent + 1_pInt
do i=2,6,2 do i=2_pInt,6_pInt,2_pInt
tag = IO_lc(IO_stringValue(line,positions,i)) tag = IO_lc(IO_stringValue(line,positions,i))
select case (tag) select case (tag)
case('phase') case('phase')
@ -357,7 +357,7 @@ subroutine material_parseMicrostructure(file,myPart)
!********************************************************************* !*********************************************************************
subroutine material_parseCrystallite(file,myPart) subroutine material_parseCrystallite(myFile,myPart)
!********************************************************************* !*********************************************************************
use prec, only: pInt use prec, only: pInt
@ -366,33 +366,33 @@ subroutine material_parseCrystallite(file,myPart)
implicit none implicit none
character(len=*), intent(in) :: myPart character(len=*), intent(in) :: myPart
integer(pInt), intent(in) :: file integer(pInt), intent(in) :: myFile
integer(pInt) Nsections, section integer(pInt) Nsections, section
character(len=1024) line character(len=1024) line
Nsections = IO_countSections(file,myPart) Nsections = IO_countSections(myFile,myPart)
material_Ncrystallite = Nsections material_Ncrystallite = Nsections
if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart) if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart)
allocate(crystallite_name(Nsections)); crystallite_name = '' allocate(crystallite_name(Nsections)); crystallite_name = ''
allocate(crystallite_Noutput(Nsections)); crystallite_Noutput = 0_pInt allocate(crystallite_Noutput(Nsections)); crystallite_Noutput = 0_pInt
crystallite_Noutput = IO_countTagInPart(file,myPart,'(output)',Nsections) crystallite_Noutput = IO_countTagInPart(myFile,myPart,'(output)',Nsections)
rewind(file) rewind(myFile)
line = '' line = ''
section = 0 section = 0_pInt
do while (IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to myPart do while (IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to myPart
read(file,'(a1024)',END=100) line read(myFile,'(a1024)',END=100) line
enddo enddo
do do
read(file,'(a1024)',END=100) line read(myFile,'(a1024)',END=100) line
if (IO_isBlank(line)) cycle ! skip empty lines if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') then ! next section if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1 section = section + 1_pInt
crystallite_name(section) = IO_getTag(line,'[',']') crystallite_name(section) = IO_getTag(line,'[',']')
endif endif
enddo enddo
@ -401,7 +401,7 @@ subroutine material_parseCrystallite(file,myPart)
!********************************************************************* !*********************************************************************
subroutine material_parsePhase(file,myPart) subroutine material_parsePhase(myFile,myPart)
!********************************************************************* !*********************************************************************
use prec, only: pInt use prec, only: pInt
@ -409,14 +409,14 @@ subroutine material_parsePhase(file,myPart)
implicit none implicit none
character(len=*), intent(in) :: myPart character(len=*), intent(in) :: myPart
integer(pInt), intent(in) :: file integer(pInt), intent(in) :: myFile
integer(pInt), parameter :: maxNchunks = 2 integer(pInt), parameter :: maxNchunks = 2_pInt
integer(pInt), dimension(1+2*maxNchunks) :: positions integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt) Nsections, section, s integer(pInt) Nsections, section, s
character(len=64) tag character(len=64) tag
character(len=1024) line character(len=1024) line
Nsections = IO_countSections(file,myPart) Nsections = IO_countSections(myFile,myPart)
material_Nphase = Nsections material_Nphase = Nsections
if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart) if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart)
@ -426,23 +426,23 @@ subroutine material_parsePhase(file,myPart)
allocate(phase_Noutput(Nsections)) allocate(phase_Noutput(Nsections))
allocate(phase_localConstitution(Nsections)) allocate(phase_localConstitution(Nsections))
phase_Noutput = IO_countTagInPart(file,myPart,'(output)',Nsections) phase_Noutput = IO_countTagInPart(myFile,myPart,'(output)',Nsections)
phase_localConstitution = .not. IO_spotTagInPart(file,myPart,'/nonlocal/',Nsections) phase_localConstitution = .not. IO_spotTagInPart(myFile,myPart,'/nonlocal/',Nsections)
rewind(file) rewind(myFile)
line = '' line = ''
section = 0 section = 0_pInt
do while (IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to myPart do while (IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to myPart
read(file,'(a1024)',END=100) line read(myFile,'(a1024)',END=100) line
enddo enddo
do do
read(file,'(a1024)',END=100) line read(myFile,'(a1024)',END=100) line
if (IO_isBlank(line)) cycle ! skip empty lines if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') then ! next section if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1 section = section + 1_pInt
phase_name(section) = IO_getTag(line,'[',']') phase_name(section) = IO_getTag(line,'[',']')
endif endif
if (section > 0_pInt) then if (section > 0_pInt) then
@ -451,9 +451,9 @@ subroutine material_parsePhase(file,myPart)
select case(tag) select case(tag)
case ('constitution') case ('constitution')
phase_constitution(section) = IO_lc(IO_stringValue(line,positions,2_pInt)) phase_constitution(section) = IO_lc(IO_stringValue(line,positions,2_pInt))
do s = 1,section do s = 1_pInt,section
if (phase_constitution(s) == phase_constitution(section)) & if (phase_constitution(s) == phase_constitution(section)) &
phase_constitutionInstance(section) = phase_constitutionInstance(section) + 1 ! count instances phase_constitutionInstance(section) = phase_constitutionInstance(section) + 1_pInt ! count instances
enddo enddo
end select end select
endif endif
@ -463,7 +463,7 @@ subroutine material_parsePhase(file,myPart)
!********************************************************************* !*********************************************************************
subroutine material_parseTexture(file,myPart) subroutine material_parseTexture(myFile,myPart)
!********************************************************************* !*********************************************************************
use prec, only: pInt, pReal use prec, only: pInt, pReal
@ -472,15 +472,15 @@ subroutine material_parseTexture(file,myPart)
implicit none implicit none
character(len=*), intent(in) :: myPart character(len=*), intent(in) :: myPart
integer(pInt), intent(in) :: file integer(pInt), intent(in) :: myFile
integer(pInt), parameter :: maxNchunks = 13 integer(pInt), parameter :: maxNchunks = 13_pInt
integer(pInt), dimension(1+2*maxNchunks) :: positions integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt) Nsections, section, gauss, fiber, i integer(pInt) Nsections, section, gauss, fiber, i
character(len=64) tag character(len=64) tag
character(len=1024) line character(len=1024) line
Nsections = IO_countSections(file,myPart) Nsections = IO_countSections(myFile,myPart)
material_Ntexture = Nsections material_Ntexture = Nsections
if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart) if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart)
@ -490,33 +490,33 @@ subroutine material_parseTexture(file,myPart)
allocate(texture_Ngauss(Nsections)); texture_Ngauss = 0_pInt allocate(texture_Ngauss(Nsections)); texture_Ngauss = 0_pInt
allocate(texture_Nfiber(Nsections)); texture_Nfiber = 0_pInt allocate(texture_Nfiber(Nsections)); texture_Nfiber = 0_pInt
texture_Ngauss = IO_countTagInPart(file,myPart,'(gauss)', Nsections) + & texture_Ngauss = IO_countTagInPart(myFile,myPart,'(gauss)', Nsections) + &
IO_countTagInPart(file,myPart,'(random)',Nsections) IO_countTagInPart(myFile,myPart,'(random)',Nsections)
texture_Nfiber = IO_countTagInPart(file,myPart,'(fiber)', Nsections) texture_Nfiber = IO_countTagInPart(myFile,myPart,'(fiber)', Nsections)
texture_maxNgauss = maxval(texture_Ngauss) texture_maxNgauss = maxval(texture_Ngauss)
texture_maxNfiber = maxval(texture_Nfiber) texture_maxNfiber = maxval(texture_Nfiber)
allocate(texture_Gauss (5,texture_maxNgauss,Nsections)); texture_Gauss = 0.0_pReal allocate(texture_Gauss (5,texture_maxNgauss,Nsections)); texture_Gauss = 0.0_pReal
allocate(texture_Fiber (6,texture_maxNfiber,Nsections)); texture_Fiber = 0.0_pReal allocate(texture_Fiber (6,texture_maxNfiber,Nsections)); texture_Fiber = 0.0_pReal
rewind(file) rewind(myFile)
line = '' line = ''
section = 0 section = 0_pInt
do while (IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to myPart do while (IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to myPart
read(file,'(a1024)',END=100) line read(myFile,'(a1024)',END=100) line
enddo enddo
do do
read(file,'(a1024)',END=100) line read(myFile,'(a1024)',END=100) line
if (IO_isBlank(line)) cycle ! skip empty lines if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') then ! next section if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1 section = section + 1_pInt
gauss = 0 gauss = 0_pInt
fiber = 0 fiber = 0_pInt
texture_name(section) = IO_getTag(line,'[',']') texture_name(section) = IO_getTag(line,'[',']')
endif endif
if (section > 0) then if (section > 0_pInt) then
positions = IO_stringPos(line,maxNchunks) positions = IO_stringPos(line,maxNchunks)
tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key
select case(tag) select case(tag)
@ -528,17 +528,17 @@ subroutine material_parseTexture(file,myPart)
tag = IO_lc(IO_stringValue(line,positions,2_pInt)) tag = IO_lc(IO_stringValue(line,positions,2_pInt))
select case (tag) select case (tag)
case('orthotropic') case('orthotropic')
texture_symmetry(section) = 4 texture_symmetry(section) = 4_pInt
case('monoclinic') case('monoclinic')
texture_symmetry(section) = 2 texture_symmetry(section) = 2_pInt
case default case default
texture_symmetry(section) = 1 texture_symmetry(section) = 1_pInt
end select end select
case ('(random)') case ('(random)')
gauss = gauss + 1 gauss = gauss + 1_pInt
texture_Gauss(1:3,gauss,section) = math_sampleRandomOri() texture_Gauss(1:3,gauss,section) = math_sampleRandomOri()
do i = 2,4,2 do i = 2_pInt,4_pInt,2_pInt
tag = IO_lc(IO_stringValue(line,positions,i)) tag = IO_lc(IO_stringValue(line,positions,i))
select case (tag) select case (tag)
case('scatter') case('scatter')
@ -549,8 +549,8 @@ subroutine material_parseTexture(file,myPart)
enddo enddo
case ('(gauss)') case ('(gauss)')
gauss = gauss + 1 gauss = gauss + 1_pInt
do i = 2,10,2 do i = 2_pInt,10_pInt,2_pInt
tag = IO_lc(IO_stringValue(line,positions,i)) tag = IO_lc(IO_stringValue(line,positions,i))
select case (tag) select case (tag)
case('phi1') case('phi1')
@ -567,8 +567,8 @@ subroutine material_parseTexture(file,myPart)
enddo enddo
case ('(fiber)') case ('(fiber)')
fiber = fiber + 1 fiber = fiber + 1_pInt
do i = 2,12,2 do i = 2_pInt,12_pInt,2_pInt
tag = IO_lc(IO_stringValue(line,positions,i)) tag = IO_lc(IO_stringValue(line,positions,i))
select case (tag) select case (tag)
case('alpha1') case('alpha1')
@ -629,7 +629,7 @@ subroutine material_populateGrains()
allocate(Nelems(material_Nhomogenization,material_Nmicrostructure)); Nelems = 0_pInt allocate(Nelems(material_Nhomogenization,material_Nmicrostructure)); Nelems = 0_pInt
! precounting of elements for each homog/micro pair ! precounting of elements for each homog/micro pair
do e = 1, mesh_NcpElems do e = 1_pInt, mesh_NcpElems
homog = mesh_element(3,e) homog = mesh_element(3,e)
micro = mesh_element(4,e) micro = mesh_element(4,e)
Nelems(homog,micro) = Nelems(homog,micro) + 1_pInt Nelems(homog,micro) = Nelems(homog,micro) + 1_pInt
@ -641,12 +641,12 @@ subroutine material_populateGrains()
Nelems = 0_pInt ! reuse as counter Nelems = 0_pInt ! reuse as counter
! identify maximum grain count per IP (from element) and find grains per homog/micro pair ! identify maximum grain count per IP (from element) and find grains per homog/micro pair
do e = 1,mesh_NcpElems do e = 1_pInt,mesh_NcpElems
homog = mesh_element(3,e) homog = mesh_element(3,e)
micro = mesh_element(4,e) micro = mesh_element(4,e)
if (homog < 1 .or. homog > material_Nhomogenization) & ! out of bounds if (homog < 1_pInt .or. homog > material_Nhomogenization) & ! out of bounds
call IO_error(154_pInt,e,0_pInt,0_pInt) call IO_error(154_pInt,e,0_pInt,0_pInt)
if (micro < 1 .or. micro > material_Nmicrostructure) & ! out of bounds if (micro < 1_pInt .or. micro > material_Nmicrostructure) & ! out of bounds
call IO_error(155_pInt,e,0_pInt,0_pInt) call IO_error(155_pInt,e,0_pInt,0_pInt)
if (microstructure_elemhomo(micro)) then if (microstructure_elemhomo(micro)) then
dGrains = homogenization_Ngrains(homog) dGrains = homogenization_Ngrains(homog)
@ -664,7 +664,7 @@ subroutine material_populateGrains()
allocate(textureOfGrain(maxval(Ngrains))) ! reserve memory for maximum case allocate(textureOfGrain(maxval(Ngrains))) ! reserve memory for maximum case
allocate(orientationOfGrain(3,maxval(Ngrains))) ! reserve memory for maximum case allocate(orientationOfGrain(3,maxval(Ngrains))) ! reserve memory for maximum case
if (debug_verbosity > 0) then if (debug_verbosity > 0_pInt) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write (6,*) write (6,*)
write (6,*) 'MATERIAL grain population' write (6,*) 'MATERIAL grain population'
@ -672,12 +672,12 @@ subroutine material_populateGrains()
write (6,'(a32,1x,a32,1x,a6)') 'homogenization_name','microstructure_name','grain#' write (6,'(a32,1x,a32,1x,a6)') 'homogenization_name','microstructure_name','grain#'
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
do homog = 1,material_Nhomogenization ! loop over homogenizations do homog = 1_pInt,material_Nhomogenization ! loop over homogenizations
dGrains = homogenization_Ngrains(homog) ! grain number per material point dGrains = homogenization_Ngrains(homog) ! grain number per material point
do micro = 1,material_Nmicrostructure ! all pairs of homog and micro do micro = 1_pInt,material_Nmicrostructure ! all pairs of homog and micro
if (Ngrains(homog,micro) > 0) then ! an active pair of homog and micro if (Ngrains(homog,micro) > 0_pInt) then ! an active pair of homog and micro
myNgrains = Ngrains(homog,micro) ! assign short name for total number of grains to populate myNgrains = Ngrains(homog,micro) ! assign short name for total number of grains to populate
if (debug_verbosity > 0) then if (debug_verbosity > 0_pInt) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write (6,*) write (6,*)
write (6,'(a32,1x,a32,1x,i6)') homogenization_name(homog),microstructure_name(micro),myNgrains write (6,'(a32,1x,a32,1x,i6)') homogenization_name(homog),microstructure_name(micro),myNgrains
@ -690,12 +690,12 @@ subroutine material_populateGrains()
do hme = 1_pInt, Nelems(homog,micro) do hme = 1_pInt, Nelems(homog,micro)
e = elemsOfHomogMicro(hme,homog,micro) ! my combination of homog and micro, only perform calculations for elements with homog, micro combinations which is indexed in cpElemsindex e = elemsOfHomogMicro(hme,homog,micro) ! my combination of homog and micro, only perform calculations for elements with homog, micro combinations which is indexed in cpElemsindex
if (microstructure_elemhomo(micro)) then ! homogeneous distribution of grains over each element's IPs if (microstructure_elemhomo(micro)) then ! homogeneous distribution of grains over each element's IPs
volumeOfGrain(grain+1:grain+dGrains) = sum(mesh_ipVolume(1:FE_Nips(mesh_element(2,e)),e))/& volumeOfGrain(grain+1_pInt:grain+dGrains) = sum(mesh_ipVolume(1:FE_Nips(mesh_element(2,e)),e))/&
real(dGrains,pReal) real(dGrains,pReal)
grain = grain + dGrains ! wind forward by NgrainsPerIP grain = grain + dGrains ! wind forward by NgrainsPerIP
else else
forall (i = 1:FE_Nips(mesh_element(2,e))) & ! loop over IPs forall (i = 1_pInt:FE_Nips(mesh_element(2,e))) & ! loop over IPs
volumeOfGrain(grain+(i-1)*dGrains+1:grain+i*dGrains) = & volumeOfGrain(grain+(i-1)*dGrains+1_pInt:grain+i*dGrains) = &
mesh_ipVolume(i,e)/dGrains ! assign IPvolume/Ngrains to all grains of IP mesh_ipVolume(i,e)/dGrains ! assign IPvolume/Ngrains to all grains of IP
grain = grain + FE_Nips(mesh_element(2,e)) * dGrains ! wind forward by Nips*NgrainsPerIP grain = grain + FE_Nips(mesh_element(2,e)) * dGrains ! wind forward by Nips*NgrainsPerIP
endif endif
@ -703,13 +703,13 @@ subroutine material_populateGrains()
! ---------------------------------------------------------------------------- divide myNgrains as best over constituents ! ---------------------------------------------------------------------------- divide myNgrains as best over constituents
NgrainsOfConstituent = 0_pInt NgrainsOfConstituent = 0_pInt
forall (i = 1:microstructure_Nconstituents(micro)) & forall (i = 1_pInt:microstructure_Nconstituents(micro)) &
NgrainsOfConstituent(i) = nint(microstructure_fraction(i,micro) * myNgrains, pInt) ! do rounding integer conversion NgrainsOfConstituent(i) = nint(microstructure_fraction(i,micro) * myNgrains, pInt) ! do rounding integer conversion
do while (sum(NgrainsOfConstituent) /= myNgrains) ! total grain count over constituents wrong? do while (sum(NgrainsOfConstituent) /= myNgrains) ! total grain count over constituents wrong?
sgn = sign(1_pInt, myNgrains - sum(NgrainsOfConstituent)) ! direction of required change sgn = sign(1_pInt, myNgrains - sum(NgrainsOfConstituent)) ! direction of required change
extreme = 0.0_pReal extreme = 0.0_pReal
t = 0_pInt t = 0_pInt
do i = 1,microstructure_Nconstituents(micro) ! find largest deviator do i = 1_pInt,microstructure_Nconstituents(micro) ! find largest deviator
if (real(sgn,pReal)*log(NgrainsOfConstituent(i)/myNgrains/microstructure_fraction(i,micro)) > extreme) then if (real(sgn,pReal)*log(NgrainsOfConstituent(i)/myNgrains/microstructure_fraction(i,micro)) > extreme) then
extreme = real(sgn,pReal)*log(NgrainsOfConstituent(i)/myNgrains/microstructure_fraction(i,micro)) extreme = real(sgn,pReal)*log(NgrainsOfConstituent(i)/myNgrains/microstructure_fraction(i,micro))
t = i t = i
@ -723,21 +723,21 @@ subroutine material_populateGrains()
orientationOfGrain = 0.0_pReal orientationOfGrain = 0.0_pReal
grain = 0_pInt ! reset microstructure grain index grain = 0_pInt ! reset microstructure grain index
do i = 1,microstructure_Nconstituents(micro) ! loop over constituents do i = 1_pInt,microstructure_Nconstituents(micro) ! loop over constituents
phaseID = microstructure_phase(i,micro) phaseID = microstructure_phase(i,micro)
textureID = microstructure_texture(i,micro) textureID = microstructure_texture(i,micro)
phaseOfGrain(grain+1:grain+NgrainsOfConstituent(i)) = phaseID ! assign resp. phase phaseOfGrain(grain+1_pInt:grain+NgrainsOfConstituent(i)) = phaseID ! assign resp. phase
textureOfGrain(grain+1:grain+NgrainsOfConstituent(i)) = textureID ! assign resp. texture textureOfGrain(grain+1_pInt:grain+NgrainsOfConstituent(i)) = textureID ! assign resp. texture
myNorientations = ceiling(real(NgrainsOfConstituent(i),pReal)/& myNorientations = ceiling(real(NgrainsOfConstituent(i),pReal)/&
real(texture_symmetry(textureID),pReal)) ! max number of unique orientations (excl. symmetry) real(texture_symmetry(textureID),pReal),pInt) ! max number of unique orientations (excl. symmetry)
constituentGrain = 0_pInt ! constituent grain index constituentGrain = 0_pInt ! constituent grain index
! --------- ! ---------
if (texture_ODFfile(textureID) == '') then ! dealing with texture components if (texture_ODFfile(textureID) == '') then ! dealing with texture components
! --------- ! ---------
do t = 1,texture_Ngauss(textureID) ! loop over Gauss components do t = 1_pInt,texture_Ngauss(textureID) ! loop over Gauss components
do g = 1,int(myNorientations*texture_Gauss(5,t,textureID)) ! loop over required grain count do g = 1_pInt,int(myNorientations*texture_Gauss(5,t,textureID),pInt) ! loop over required grain count
orientationOfGrain(:,grain+constituentGrain+g) = & orientationOfGrain(:,grain+constituentGrain+g) = &
math_sampleGaussOri(texture_Gauss(1:3,t,textureID),& math_sampleGaussOri(texture_Gauss(1:3,t,textureID),&
texture_Gauss( 4,t,textureID)) texture_Gauss( 4,t,textureID))
@ -745,17 +745,17 @@ subroutine material_populateGrains()
constituentGrain = constituentGrain + int(myNorientations*texture_Gauss(5,t,textureID)) constituentGrain = constituentGrain + int(myNorientations*texture_Gauss(5,t,textureID))
enddo enddo
do t = 1,texture_Nfiber(textureID) ! loop over fiber components do t = 1_pInt,texture_Nfiber(textureID) ! loop over fiber components
do g = 1,int(myNorientations*texture_Fiber(6,t,textureID)) ! loop over required grain count do g = 1_pInt,int(myNorientations*texture_Fiber(6,t,textureID),pInt) ! loop over required grain count
orientationOfGrain(:,grain+constituentGrain+g) = & orientationOfGrain(:,grain+constituentGrain+g) = &
math_sampleFiberOri(texture_Fiber(1:2,t,textureID),& math_sampleFiberOri(texture_Fiber(1:2,t,textureID),&
texture_Fiber(3:4,t,textureID),& texture_Fiber(3:4,t,textureID),&
texture_Fiber( 5,t,textureID)) texture_Fiber( 5,t,textureID))
enddo enddo
constituentGrain = constituentGrain + int(myNorientations*texture_fiber(6,t,textureID)) constituentGrain = constituentGrain + int(myNorientations*texture_fiber(6,t,textureID),pInt)
enddo enddo
do j = constituentGrain+1,myNorientations ! fill remainder with random do j = constituentGrain+1_pInt,myNorientations ! fill remainder with random
orientationOfGrain(:,grain+j) = math_sampleRandomOri() orientationOfGrain(:,grain+j) = math_sampleRandomOri()
enddo enddo
! --------- ! ---------
@ -770,12 +770,12 @@ subroutine material_populateGrains()
symExtension = texture_symmetry(textureID) - 1_pInt symExtension = texture_symmetry(textureID) - 1_pInt
if (symExtension > 0_pInt) then ! sample symmetry if (symExtension > 0_pInt) then ! sample symmetry
constituentGrain = NgrainsOfConstituent(i)-myNorientations ! calc remainder of array constituentGrain = NgrainsOfConstituent(i)-myNorientations ! calc remainder of array
do j = 1,myNorientations ! loop over each "real" orientation do j = 1_pInt,myNorientations ! loop over each "real" orientation
symOrientation = math_symmetricEulers(texture_symmetry(textureID),orientationOfGrain(:,j)) ! get symmetric equivalents symOrientation = math_symmetricEulers(texture_symmetry(textureID),orientationOfGrain(:,j)) ! get symmetric equivalents
e = min(symExtension,constituentGrain) ! are we at end of constituent grain array? e = min(symExtension,constituentGrain) ! are we at end of constituent grain array?
if (e > 0_pInt) then if (e > 0_pInt) then
orientationOfGrain(:,grain+myNorientations+1+(j-1)*symExtension:& orientationOfGrain(:,grain+myNorientations+1+(j-1_pInt)*symExtension:&
grain+myNorientations+e+(j-1)*symExtension) = & grain+myNorientations+e+(j-1_pInt)*symExtension) = &
symOrientation(:,1:e) symOrientation(:,1:e)
constituentGrain = constituentGrain - e ! remainder shrinks by e constituentGrain = constituentGrain - e ! remainder shrinks by e
endif endif
@ -787,7 +787,7 @@ subroutine material_populateGrains()
! ---------------------------------------------------------------------------- ! ----------------------------------------------------------------------------
if (.not. microstructure_elemhomo(micro)) then ! unless element homogeneous, reshuffle grains if (.not. microstructure_elemhomo(micro)) then ! unless element homogeneous, reshuffle grains
do i=1,myNgrains-1 ! walk thru grains do i=1_pInt,myNgrains-1_pInt ! walk thru grains
call random_number(rnd) call random_number(rnd)
t = nint(rnd*(myNgrains-i)+i+0.5_pReal,pInt) ! select a grain in remaining list t = nint(rnd*(myNgrains-i)+i+0.5_pReal,pInt) ! select a grain in remaining list
m = phaseOfGrain(t) ! exchange current with random m = phaseOfGrain(t) ! exchange current with random
@ -809,7 +809,7 @@ subroutine material_populateGrains()
do hme = 1_pInt, Nelems(homog,micro) do hme = 1_pInt, Nelems(homog,micro)
e = elemsOfHomogMicro(hme,homog,micro) ! only perform calculations for elements with homog, micro combinations which is indexed in cpElemsindex e = elemsOfHomogMicro(hme,homog,micro) ! only perform calculations for elements with homog, micro combinations which is indexed in cpElemsindex
if (microstructure_elemhomo(micro)) then ! homogeneous distribution of grains over each element's IPs if (microstructure_elemhomo(micro)) then ! homogeneous distribution of grains over each element's IPs
forall (i = 1:FE_Nips(mesh_element(2,e)), g = 1:dGrains) ! loop over IPs and grains forall (i = 1_pInt:FE_Nips(mesh_element(2,e)), g = 1_pInt:dGrains) ! loop over IPs and grains
material_volume(g,i,e) = volumeOfGrain(grain+g) material_volume(g,i,e) = volumeOfGrain(grain+g)
material_phase(g,i,e) = phaseOfGrain(grain+g) material_phase(g,i,e) = phaseOfGrain(grain+g)
material_texture(g,i,e) = textureOfGrain(grain+g) material_texture(g,i,e) = textureOfGrain(grain+g)
@ -818,11 +818,11 @@ subroutine material_populateGrains()
FEsolving_execIP(2,e) = 1_pInt ! restrict calculation to first IP only, since all other results are to be copied from this FEsolving_execIP(2,e) = 1_pInt ! restrict calculation to first IP only, since all other results are to be copied from this
grain = grain + dGrains ! wind forward by NgrainsPerIP grain = grain + dGrains ! wind forward by NgrainsPerIP
else else
forall (i = 1:FE_Nips(mesh_element(2,e)), g = 1:dGrains) ! loop over IPs and grains forall (i = 1_pInt:FE_Nips(mesh_element(2,e)), g = 1_pInt:dGrains) ! loop over IPs and grains
material_volume(g,i,e) = volumeOfGrain(grain+(i-1)*dGrains+g) material_volume(g,i,e) = volumeOfGrain(grain+(i-1_pInt)*dGrains+g)
material_phase(g,i,e) = phaseOfGrain(grain+(i-1)*dGrains+g) material_phase(g,i,e) = phaseOfGrain(grain+(i-1_pInt)*dGrains+g)
material_texture(g,i,e) = textureOfGrain(grain+(i-1)*dGrains+g) material_texture(g,i,e) = textureOfGrain(grain+(i-1_pInt)*dGrains+g)
material_EulerAngles(:,g,i,e) = orientationOfGrain(:,grain+(i-1)*dGrains+g) material_EulerAngles(:,g,i,e) = orientationOfGrain(:,grain+(i-1_pInt)*dGrains+g)
end forall end forall
grain = grain + FE_Nips(mesh_element(2,e)) * dGrains ! wind forward by Nips*NgrainsPerIP grain = grain + FE_Nips(mesh_element(2,e)) * dGrains ! wind forward by Nips*NgrainsPerIP
endif endif

6
code/math.f90
View File

@ -33,7 +33,7 @@
real(pReal), parameter :: pi = 3.14159265358979323846264338327950288419716939937510_pReal real(pReal), parameter :: pi = 3.14159265358979323846264338327950288419716939937510_pReal
real(pReal), parameter :: inDeg = 180.0_pReal/pi real(pReal), parameter :: inDeg = 180.0_pReal/pi
real(pReal), parameter :: inRad = pi/180.0_pReal real(pReal), parameter :: inRad = pi/180.0_pReal
complex(pReal), parameter :: two_pi_img = cmplx(0.0_pReal,2.0_pReal* pi, pReal) complex(pReal), parameter :: two_pi_img = (0.0_pReal,2.0_pReal)* pi
! *** 3x3 Identity *** ! *** 3x3 Identity ***
real(pReal), dimension(3,3), parameter :: math_I3 = & real(pReal), dimension(3,3), parameter :: math_I3 = &
@ -265,7 +265,7 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = &
integer(pInt), intent(in) :: istart,iend integer(pInt), intent(in) :: istart,iend
integer(pInt) :: d,i,j,k,x,tmp integer(pInt) :: d,i,j,k,x,tmp
d = size(a,1_pInt) ! number of linked data d = int(size(a,1_pInt), pInt) ! number of linked data
! set the starting and ending points, and the pivot point ! set the starting and ending points, and the pivot point
i = istart i = istart
@ -3433,7 +3433,7 @@ subroutine curl_fft(res,geomdim,vec_tens,field,curl)
xi(i,j,k,1:3) = real(k_s, pReal)/geomdim xi(i,j,k,1:3) = real(k_s, pReal)/geomdim
enddo; enddo; enddo enddo; enddo; enddo
do k = 1_pInt, res(3); do j = 1_pInt, res(2);do i = 1_pInt, res1_red do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red
do l = 1_pInt, vec_tens do l = 1_pInt, vec_tens
curl_fourier(i,j,k,l,1) = ( field_fourier(i,j,k,l,3)*xi(i,j,k,2)& curl_fourier(i,j,k,l,1) = ( field_fourier(i,j,k,l,3)*xi(i,j,k,2)&
-field_fourier(i,j,k,l,2)*xi(i,j,k,3) )*two_pi_img -field_fourier(i,j,k,l,2)*xi(i,j,k,3) )*two_pi_img

653
code/mesh.f90
View File

File diff suppressed because it is too large Load Diff

37
code/numerics.f90
View File

@ -69,18 +69,18 @@ real(pReal) :: relevantStrain = 1.0e-7_pReal, &
err_stress_tolrel = 0.01_pReal , & ! relative tolerance for fullfillment of stress BC, Default: 0.01 allowing deviation of 1% of maximum stress err_stress_tolrel = 0.01_pReal , & ! relative tolerance for fullfillment of stress BC, Default: 0.01 allowing deviation of 1% of maximum stress
fftw_timelimit = -1.0_pReal, & ! sets the timelimit of plan creation for FFTW, see manual on www.fftw.org, Default -1.0: disable timelimit fftw_timelimit = -1.0_pReal, & ! sets the timelimit of plan creation for FFTW, see manual on www.fftw.org, Default -1.0: disable timelimit
rotation_tol = 1.0e-12_pReal ! tolerance of rotation specified in loadcase, Default 1.0e-12: first guess rotation_tol = 1.0e-12_pReal ! tolerance of rotation specified in loadcase, Default 1.0e-12: first guess
character(len=64) :: fftw_planner_string = 'FFTW_PATIENT' ! reads the planing-rigor flag, see manual on www.fftw.org, Default FFTW_PATIENT: use patiant planner flag character(len=64) :: fftw_plan_mode = 'FFTW_PATIENT' ! reads the planing-rigor flag, see manual on www.fftw.org, Default FFTW_PATIENT: use patiant planner flag
integer(pInt) :: fftw_planner_flag = -1_pInt ! conversion of fftw_planner_string to integer, basically what is usually done in the include file of fftw integer(pInt) :: fftw_planner_flag = -1_pInt, & ! conversion of fftw_plan_mode to integer, basically what is usually done in the include file of fftw
logical :: memory_efficient = .true. ,& ! for fast execution (pre calculation of gamma_hat), Default .true.: do not precalculate itmax = 20_pInt ! maximum number of iterations
divergence_correction = .false. ,& ! correct divergence calculation in fourier space, Default .false.: no correction logical :: memory_efficient = .true., & ! for fast execution (pre calculation of gamma_hat), Default .true.: do not precalculate
update_gamma = .false.,& ! update gamma operator with current stiffness, Default .false.: use initial stiffness divergence_correction = .false., & ! correct divergence calculation in fourier space, Default .false.: no correction
update_gamma = .false., & ! update gamma operator with current stiffness, Default .false.: use initial stiffness
simplified_algorithm = .true. ! use short algorithm without fluctuation field, Default .true.: use simplified algorithm simplified_algorithm = .true. ! use short algorithm without fluctuation field, Default .true.: use simplified algorithm
real(pReal) :: cut_off_value = 0.0_pReal ! percentage of frequencies to cut away, Default 0.0: use all frequencies
integer(pInt) :: itmax = 20_pInt , & ! maximum number of iterations
!* Random seeding parameters !* Random seeding parameters
fixedSeed = 0_pInt ! fixed seeding for pseudo-random number generator, Default 0: use random seed integer(pInt) :: fixedSeed = 0_pInt ! fixed seeding for pseudo-random number generator, Default 0: use random seed
!* OpenMP variable !* OpenMP variable
integer(pInt) :: DAMASK_NumThreadsInt = 0_pInt ! value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive integer(pInt) :: DAMASK_NumThreadsInt = 0_pInt ! value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive
@ -111,7 +111,7 @@ subroutine numerics_init()
!*** local variables ***! !*** local variables ***!
integer(pInt), parameter :: fileunit = 300_pInt integer(pInt), parameter :: fileunit = 300_pInt
integer(pInt), parameter :: maxNchunks = 2_pInt integer(pInt), parameter :: maxNchunks = 2_pInt
integer(pInt) :: gotDAMASK_NUM_THREADS = 1_pInt integer :: gotDAMASK_NUM_THREADS = 1
integer(pInt), dimension(1+2*maxNchunks) :: positions integer(pInt), dimension(1+2*maxNchunks) :: positions
character(len=64) :: tag character(len=64) :: tag
character(len=1024) :: line character(len=1024) :: line
@ -127,9 +127,9 @@ subroutine numerics_init()
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
!$ call GET_ENVIRONMENT_VARIABLE(NAME='DAMASK_NUM_THREADS',VALUE=DAMASK_NumThreadsString,STATUS=gotDAMASK_NUM_THREADS) ! get environment variable DAMASK_NUM_THREADS... !$ call GET_ENVIRONMENT_VARIABLE(NAME='DAMASK_NUM_THREADS',VALUE=DAMASK_NumThreadsString,STATUS=gotDAMASK_NUM_THREADS) ! get environment variable DAMASK_NUM_THREADS...
!$ if(gotDAMASK_NUM_THREADS /= 0_pInt) call IO_warning(47_pInt,ext_msg=DAMASK_NumThreadsString) !$ if(gotDAMASK_NUM_THREADS /= 0) call IO_warning(47_pInt,ext_msg=DAMASK_NumThreadsString)
!$ read(DAMASK_NumThreadsString,'(i6)') DAMASK_NumThreadsInt ! ...convert it to integer... !$ read(DAMASK_NumThreadsString,'(i6)') DAMASK_NumThreadsInt ! ...convert it to integer...
!$ if (DAMASK_NumThreadsInt < 1) DAMASK_NumThreadsInt = 1 ! ...ensure that its at least one... !$ if (DAMASK_NumThreadsInt < 1_pInt) DAMASK_NumThreadsInt = 1_pInt ! ...ensure that its at least one...
!$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! ...and use it as number of threads for parallel execution !$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! ...and use it as number of threads for parallel execution
! try to open the config file ! try to open the config file
@ -238,8 +238,8 @@ subroutine numerics_init()
memory_efficient = IO_intValue(line,positions,2_pInt) > 0_pInt memory_efficient = IO_intValue(line,positions,2_pInt) > 0_pInt
case ('fftw_timelimit') case ('fftw_timelimit')
fftw_timelimit = IO_floatValue(line,positions,2_pInt) fftw_timelimit = IO_floatValue(line,positions,2_pInt)
case ('fftw_planner_string') case ('fftw_plan_mode')
fftw_planner_string = IO_stringValue(line,positions,2_pInt) fftw_plan_mode = IO_stringValue(line,positions,2_pInt)
case ('rotation_tol') case ('rotation_tol')
rotation_tol = IO_floatValue(line,positions,2_pInt) rotation_tol = IO_floatValue(line,positions,2_pInt)
case ('divergence_correction') case ('divergence_correction')
@ -248,8 +248,6 @@ subroutine numerics_init()
update_gamma = IO_intValue(line,positions,2_pInt) > 0_pInt update_gamma = IO_intValue(line,positions,2_pInt) > 0_pInt
case ('simplified_algorithm') case ('simplified_algorithm')
simplified_algorithm = IO_intValue(line,positions,2_pInt) > 0_pInt simplified_algorithm = IO_intValue(line,positions,2_pInt) > 0_pInt
case ('cut_off_value')
cut_off_value = IO_floatValue(line,positions,2_pInt)
!* Random seeding parameters !* Random seeding parameters
@ -272,7 +270,7 @@ subroutine numerics_init()
endif endif
select case(IO_lc(fftw_planner_string)) ! setting parameters for the plan creation of FFTW. Basically a translation from fftw3.f select case(IO_lc(fftw_plan_mode)) ! setting parameters for the plan creation of FFTW. Basically a translation from fftw3.f
case('estimate','fftw_estimate') ! ordered from slow execution (but fast plan creation) to fast execution case('estimate','fftw_estimate') ! ordered from slow execution (but fast plan creation) to fast execution
fftw_planner_flag = 64_pInt fftw_planner_flag = 64_pInt
case('measure','fftw_measure') case('measure','fftw_measure')
@ -282,7 +280,7 @@ subroutine numerics_init()
case('exhaustive','fftw_exhaustive') case('exhaustive','fftw_exhaustive')
fftw_planner_flag = 8_pInt fftw_planner_flag = 8_pInt
case default case default
call IO_warning(warning_ID=47_pInt,ext_msg=trim(IO_lc(fftw_planner_string))) call IO_warning(warning_ID=47_pInt,ext_msg=trim(IO_lc(fftw_plan_mode)))
fftw_planner_flag = 32_pInt fftw_planner_flag = 32_pInt
end select end select
@ -341,13 +339,12 @@ subroutine numerics_init()
else else
write(6,'(a24,1x,e8.1)') ' fftw_timelimit: ',fftw_timelimit write(6,'(a24,1x,e8.1)') ' fftw_timelimit: ',fftw_timelimit
endif endif
write(6,'(a24,1x,a)') ' fftw_planner_string: ',trim(fftw_planner_string) write(6,'(a24,1x,a)') ' fftw_plan_mode: ',trim(fftw_plan_mode)
write(6,'(a24,1x,i8)') ' fftw_planner_flag: ',fftw_planner_flag write(6,'(a24,1x,i8)') ' fftw_planner_flag: ',fftw_planner_flag
write(6,'(a24,1x,e8.1)') ' rotation_tol: ',rotation_tol write(6,'(a24,1x,e8.1)') ' rotation_tol: ',rotation_tol
write(6,'(a24,1x,L8,/)') ' divergence_correction: ',divergence_correction write(6,'(a24,1x,L8,/)') ' divergence_correction: ',divergence_correction
write(6,'(a24,1x,L8,/)') ' update_gamma: ',update_gamma write(6,'(a24,1x,L8,/)') ' update_gamma: ',update_gamma
write(6,'(a24,1x,L8,/)') ' simplified_algorithm: ',simplified_algorithm write(6,'(a24,1x,L8,/)') ' simplified_algorithm: ',simplified_algorithm
write(6,'(a24,1x,e8.1)') ' cut_off_value: ',cut_off_value
!* Random seeding parameters !* Random seeding parameters
@ -411,7 +408,7 @@ subroutine numerics_init()
if (fixedSeed <= 0_pInt) then if (fixedSeed <= 0_pInt) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(a)') 'Random is random!' write(6,'(a)') ' Random is random!'
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
endsubroutine endsubroutine

4
code/prec.f90
View File

@ -35,13 +35,13 @@ real(pReal), parameter, public :: tol_gravityNodePos = 1.0e-100_pReal
! from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html ! from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html
! copy can be found in documentation/Code/Fortran ! copy can be found in documentation/Code/Fortran
#ifdef __INTEL_COMPILER #ifdef __INTEL_COMPILER
#if __INTEL_COMPILER<12000 #if __INTEL_COMPILER<1200
real(pReal), parameter, public :: DAMASK_NaN = Z'7FF0000000000001' real(pReal), parameter, public :: DAMASK_NaN = Z'7FF0000000000001'
#else #else
real(pReal), parameter, public :: DAMASK_NaN = real(Z'7FF0000000000001', pReal) real(pReal), parameter, public :: DAMASK_NaN = real(Z'7FF0000000000001', pReal)
#endif #endif
#else #else
real(pReal), parameter, public :: DAMASK_NaN = real(Z'7FF0000000000001', pReal) real(pReal), parameter, public :: DAMASK_NaN = real(Z'7FF0000000000001', pReal)
#endif #endif
type :: p_vec type :: p_vec
real(pReal), dimension(:), pointer :: p real(pReal), dimension(:), pointer :: p

4
code/prec_single.f90
View File

@ -35,13 +35,13 @@ real(pReal), parameter, public :: tol_gravityNodePos = 1.0e-36_pReal
! from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html ! from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html
! copy can be found in documentation/Code/Fortran ! copy can be found in documentation/Code/Fortran
#ifdef __INTEL_COMPILER #ifdef __INTEL_COMPILER
#if __INTEL_COMPILER<12000 #if __INTEL_COMPILER<1200
real(pReal), parameter, public :: DAMASK_NaN = Z'Z'7F800001', pReal' real(pReal), parameter, public :: DAMASK_NaN = Z'Z'7F800001', pReal'
#else #else
real(pReal), parameter, public :: DAMASK_NaN = real(Z'7F800001', pReal) real(pReal), parameter, public :: DAMASK_NaN = real(Z'7F800001', pReal)
#endif #endif
#else #else
real(pReal), parameter, public :: DAMASK_NaN = real(Z'7F800001', pReal) real(pReal), parameter, public :: DAMASK_NaN = real(Z'7F800001', pReal)
#endif #endif
type :: p_vec type :: p_vec
real(pReal), dimension(:), pointer :: p real(pReal), dimension(:), pointer :: p