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

made code standard conform to Fortran 2008 (ignoring warning concerning comments beyond character 132). Basically, changing "x" format specifier to "Nx" ("1x") plus removing $ format specifier 

added compiler switches for gfortran and ifort to check for standard conformity
old gnu compilers <4.4 are not longer supported because they don't provide the c binding for fftw
This commit is contained in:
Martin Diehl 2012-01-31 19:18:55 +00:00
parent 93f311f2b9
commit 800e291240
25 changed files with 472 additions and 470 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
code/CPFEM.f90
View File

@ -207,9 +207,9 @@ subroutine CPFEM_init()
write(6,*) '$Id$'
#include "compilation_info.f90"
if (debug_verbosity > 0) then
write(6,'(a32,x,6(i8,x))') 'CPFEM_cs: ', shape(CPFEM_cs)
write(6,'(a32,x,6(i8,x))') 'CPFEM_dcsdE: ', shape(CPFEM_dcsdE)
write(6,'(a32,x,6(i8,x))') 'CPFEM_dcsdE_knownGood: ', shape(CPFEM_dcsdE_knownGood)
write(6,'(a32,1x,6(i8,1x))') 'CPFEM_cs: ', shape(CPFEM_cs)
write(6,'(a32,1x,6(i8,1x))') 'CPFEM_dcsdE: ', shape(CPFEM_dcsdE)
write(6,'(a32,1x,6(i8,1x))') 'CPFEM_dcsdE_knownGood: ', shape(CPFEM_dcsdE_knownGood)
write(6,*)
write(6,*) 'parallelExecution: ', parallelExecution
write(6,*) 'symmetricSolver: ', symmetricSolver
@ -367,11 +367,11 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
!$OMP CRITICAL (write2out)
write(6,*)
write(6,'(a)') '#############################################'
write(6,'(a1,a22,x,f15.7,a6)') '#','theTime',theTime,'#'
write(6,'(a1,a22,x,f15.7,a6)') '#','theDelta',theDelta,'#'
write(6,'(a1,a22,x,i8,a13)') '#','theInc',theInc,'#'
write(6,'(a1,a22,x,i8,a13)') '#','cycleCounter',cycleCounter,'#'
write(6,'(a1,a22,x,i8,a13)') '#','computationMode',mode,'#'
write(6,'(a1,a22,1x,f15.7,a6)') '#','theTime',theTime,'#'
write(6,'(a1,a22,1x,f15.7,a6)') '#','theDelta',theDelta,'#'
write(6,'(a1,a22,1x,i8,a13)') '#','theInc',theInc,'#'
write(6,'(a1,a22,1x,i8,a13)') '#','cycleCounter',cycleCounter,'#'
write(6,'(a1,a22,1x,i8,a13)') '#','computationMode',mode,'#'
write(6,'(a)') '#############################################'
write(6,*)
call flush (6)
@ -404,7 +404,7 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
!$OMP CRITICAL (write2out)
write(6,'(a)') '<< CPFEM >> Aging states'
if (debug_e == cp_en .and. debug_i == IP) then
write(6,'(a,x,i8,x,i2,x,i4,/,(12(x),6(e20.8,x)))') '<< CPFEM >> AGED state of element ip grain',&
write(6,'(a,1x,i8,1x,i2,1x,i4,/,(12x,6(e20.8,1x)))') '<< CPFEM >> AGED state of element ip grain',&
cp_en, IP, 1, constitutive_state(1,IP,cp_en)%p
write(6,*)
endif
@ -492,9 +492,9 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
if (.not. terminallyIll .and. .not. outdatedFFN1) then
if (debug_verbosity > 0) then
!$OMP CRITICAL (write2out)
write(6,'(a,x,i8,x,i2)') '<< CPFEM >> OUTDATED at element ip',cp_en,IP
write(6,'(a,/,3(12(x),3(f10.6,x),/))') '<< CPFEM >> FFN1 old:',math_transpose33(materialpoint_F(1:3,1:3,IP,cp_en))
write(6,'(a,/,3(12(x),3(f10.6,x),/))') '<< CPFEM >> FFN1 now:',math_transpose33(ffn1)
write(6,'(a,1x,i8,1x,i2)') '<< CPFEM >> OUTDATED at element ip',cp_en,IP
write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 old:',math_transpose33(materialpoint_F(1:3,1:3,IP,cp_en))
write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 now:',math_transpose33(ffn1)
!$OMP END CRITICAL (write2out)
endif
outdatedFFN1 = .true.
@ -519,7 +519,7 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
FEsolving_execIP(2,cp_en) = IP
if (debug_verbosity > 0) then
!$OMP CRITICAL (write2out)
write(6,'(a,i8,x,i2)') '<< CPFEM >> Calculation for element ip ',cp_en,IP
write(6,'(a,i8,1x,i2)') '<< CPFEM >> Calculation for element ip ',cp_en,IP
!$OMP END CRITICAL (write2out)
endif
call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent
@ -645,8 +645,8 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
if (mode < 3 .and. debug_verbosity > 0 .and. ((debug_e == cp_en .and. debug_i == IP) .or. .not. debug_selectiveDebugger)) then
!$OMP CRITICAL (write2out)
write(6,'(a,i8,x,i2,/,12(x),6(f10.3,x)/)') '<< CPFEM >> stress/MPa at el ip ', cp_en, IP, cauchyStress/1e6
write(6,'(a,i8,x,i2,/,6(12(x),6(f10.3,x)/))') '<< CPFEM >> jacobian/GPa at el ip ', cp_en, IP, transpose(jacobian)/1e9
write(6,'(a,i8,1x,i2,/,12x,6(f10.3,1x)/)') '<< CPFEM >> stress/MPa at el ip ', cp_en, IP, cauchyStress/1e6
write(6,'(a,i8,1x,i2,/,6(12x,6(f10.3,1x)/))') '<< CPFEM >> jacobian/GPa at el ip ', cp_en, IP, transpose(jacobian)/1e9
call flush(6)
!$OMP END CRITICAL (write2out)
endif

89
code/DAMASK_spectral.f90
View File

@ -148,8 +148,6 @@ program DAMASK_spectral
! loop variables, convergence etc.
real(pReal) :: time = 0.0_pReal, time0 = 0.0_pReal, timeinc ! elapsed time, begin of interval, time interval
real(pReal) :: guessmode, err_div, err_stress, err_stress_tol
complex(pReal), parameter :: differentationFactor = cmplx(0.0_pReal,1.0_pReal)*& ! cmplx(0.0_pReal,2.0_pReal*pi) will give rounding error (wrong type?)
2.0_pReal * pi
real(pReal), dimension(3,3), parameter :: ones = 1.0_pReal, zeroes = 0.0_pReal
complex(pReal), dimension(3) :: temp3_Complex
complex(pReal), dimension(3,3) :: temp33_Complex
@ -413,13 +411,14 @@ program DAMASK_spectral
else
print '(a)','deformation gradient rate:'
endif
print '(3(3(f12.6,x)/)$)', 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))
print '(a,/,3(3(f12.6,x)/)$)','stress / GPa:',1e-9*merge(math_transpose33(bc(loadcase)%stress)&
,reshape(spread(DAMASK_NaN,1,9),(/3,3/))&
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/))&
,transpose(bc(loadcase)%maskStress))
if (any(bc(loadcase)%rotation /= math_I3)) &
print '(a,3(3(f12.6,x)/)$)','rotation of loadframe:',math_transpose33(bc(loadcase)%rotation)
write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') ' rotation of loadframe:',&
math_transpose33(bc(loadcase)%rotation)
print '(a,f12.6)','temperature:',bc(loadcase)%temperature
print '(a,f12.6)','time: ',bc(loadcase)%time
print '(a,i5)' ,'increments: ',bc(loadcase)%incs
@ -667,25 +666,25 @@ program DAMASK_spectral
! write header of output file
open(538,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())&
//'.spectralOut',form='UNFORMATTED',status='REPLACE')
write(538), 'load', trim(getLoadcaseName())
write(538), 'workingdir', trim(getSolverWorkingDirectoryName())
write(538), 'geometry', trim(getSolverJobName())//InputFileExtension
write(538), 'resolution', res
write(538), 'dimension', geomdim
write(538), 'materialpoint_sizeResults', materialpoint_sizeResults
write(538), 'loadcases', N_Loadcases
write(538), 'frequencies', bc(1:N_Loadcases)%outputfrequency ! one entry per loadcase
write(538), 'times', bc(1:N_Loadcases)%time ! one entry per loadcase
write(538), 'logscales', bc(1:N_Loadcases)%logscale
write(538) 'load', trim(getLoadcaseName())
write(538) 'workingdir', trim(getSolverWorkingDirectoryName())
write(538) 'geometry', trim(getSolverJobName())//InputFileExtension
write(538) 'resolution', res
write(538) 'dimension', geomdim
write(538) 'materialpoint_sizeResults', materialpoint_sizeResults
write(538) 'loadcases', N_Loadcases
write(538) 'frequencies', bc(1:N_Loadcases)%outputfrequency ! one entry per loadcase
write(538) 'times', bc(1:N_Loadcases)%time ! one entry per loadcase
write(538) 'logscales', bc(1:N_Loadcases)%logscale
bc(1)%incs = bc(1)%incs + 1_pInt ! additional for zero deformation
write(538), 'increments', bc(1:N_Loadcases)%incs ! one entry per loadcase
write(538) 'increments', bc(1:N_Loadcases)%incs ! one entry per loadcase
bc(1)%incs = bc(1)%incs - 1_pInt
write(538), 'startingIncrement', restartReadInc ! start with writing out the previous inc
write(538) 'startingIncrement', restartReadInc ! start with writing out the previous inc
write(538), 'eoh' ! end of header
write(538), materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! initial (non-deformed or read-in) results
endif
write(538) 'eoh' ! end of header
write(538) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! initial (non-deformed or read-in) results
if (debugGeneral) print '(a)' , 'Header of result file written out'
endif
if(totalIncsCounter > restartReadInc) then ! Do calculations (otherwise just forwarding)
if(bc(loadcase)%restartFrequency>0_pInt) &
@ -793,7 +792,7 @@ program DAMASK_spectral
do n = 1_pInt,3_pInt; do m = 1_pInt,3_pInt
defgrad_av_lab(m,n) = sum(defgrad(1:res(1),1:res(2),1:res(3),m,n)) * wgt
enddo; enddo
print '(a,/,3(3(f12.7,x)/)$)', 'deformation gradient:',&
write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient:',&
math_transpose33(math_rotate_forward33(defgrad_av_lab,bc(loadcase)%rotation))
print '(a)', ''
print '(a)', '... update stress field P(F) ................................'
@ -868,8 +867,8 @@ program DAMASK_spectral
! comparing 1 and 3x3 FT results
if (debugFFTW) then
call fftw_execute_dft_r2c(plan_scalarField_forth,scalarField_real,scalarField_complex)
print '(a,i1,x,i1)', 'checking FT results of compontent ', row, column
print '(a,2(es10.4,x))', 'max FT relative error ',&
print '(a,i1,1x,i1)', 'checking FT results of compontent ', row, column
print '(a,2(es10.4,1x))', 'max FT relative error ',&
maxval( real((scalarField_complex(1:res1_red,1:res(2),1:res(3))-&
tensorField_complex(1:res1_red,1:res(2),1:res(3),row,column))/&
scalarField_complex(1:res1_red,1:res(2),1:res(3)))), &
@ -887,19 +886,19 @@ program DAMASK_spectral
do i = 2_pInt, res1_red -1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice.
err_div_RMS = err_div_RMS &
+ 2.0_pReal*(sum (real(math_mul33x3_complex(tensorField_complex(i,j,k,1:3,1:3),& ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2. do not take square root and square again
xi(1:3,i,j,k))*differentationFactor)**2.0_pReal)&! --> sum squared L_2 norm of vector
xi(1:3,i,j,k))*two_pi_img)**2.0_pReal)& ! --> sum squared L_2 norm of vector
+sum(aimag(math_mul33x3_complex(tensorField_complex(i,j,k,1:3,1:3),&
xi(1:3,i,j,k))*differentationFactor)**2.0_pReal))
xi(1:3,i,j,k))*two_pi_img)**2.0_pReal))
enddo
err_div_RMS = err_div_RMS & ! Those two layers do not have a conjugate complex counterpart
+ sum(real(math_mul33x3_complex(tensorField_complex(1 ,j,k,1:3,1:3),&
xi(1:3,1 ,j,k))*differentationFactor)**2.0_pReal)&
xi(1:3,1 ,j,k))*two_pi_img)**2.0_pReal)&
+ sum(aimag(math_mul33x3_complex(tensorField_complex(1 ,j,k,1:3,1:3),&
xi(1:3,1 ,j,k))*differentationFactor)**2.0_pReal)&
xi(1:3,1 ,j,k))*two_pi_img)**2.0_pReal)&
+ sum(real(math_mul33x3_complex(tensorField_complex(res1_red,j,k,1:3,1:3),&
xi(1:3,res1_red,j,k))*differentationFactor)**2.0_pReal)&
xi(1:3,res1_red,j,k))*two_pi_img)**2.0_pReal)&
+ sum(aimag(math_mul33x3_complex(tensorField_complex(res1_red,j,k,1:3,1:3),&
xi(1:3,res1_red,j,k))*differentationFactor)**2.0_pReal)
xi(1:3,res1_red,j,k))*two_pi_img)**2.0_pReal)
enddo; enddo
err_div_RMS = sqrt(err_div_RMS)*wgt ! RMS in real space calculated with Parsevals theorem from Fourier space
err_div = err_div_RMS/p_hat_avg/sqrt(wgt) * correctionFactor ! criterion to stop iterations
@ -910,7 +909,7 @@ program DAMASK_spectral
err_div_max = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red
temp3_Complex = math_mul33x3_complex(tensorField_complex(i,j,k,1:3,1:3),&
xi(1:3,i,j,k))*differentationFactor
xi(1:3,i,j,k))*two_pi_img
err_div_max = max(err_div_max,sqrt(sum(abs(temp3_Complex)**2.0_pReal)))
divergence_complex(i,j,k,1:3) = temp3_Complex ! need divergence NOT squared
enddo; enddo; enddo
@ -991,7 +990,7 @@ program DAMASK_spectral
!--------------------------------------------------------------------------------------------------
! comparing 1 and 3x3 inverse FT results
if (debugFFTW) then
print '(a,i1,x,i1)', 'checking iFT results of compontent ', row, column
print '(a,i1,1x,i1)', 'checking iFT results of compontent ', row, column
call fftw_execute_dft_c2r(plan_scalarField_back,scalarField_complex,scalarField_real)
print '(a,es10.4)', 'max iFT relative error ',&
maxval((scalarField_real(1:res(1),1:res(2),1:res(3))-&
@ -1012,11 +1011,11 @@ program DAMASK_spectral
maxval(math_skew33(tensorField_real(i,j,k,1:3,1:3))))
temp33_Real = temp33_Real + tensorField_real(i,j,k,1:3,1:3)
enddo; enddo; enddo
print '(a,x,es10.4)' , 'max symmetrix correction of deformation:',&
print '(a,1x,es10.4)' , 'max symmetrix correction of deformation:',&
maxCorrectionSym*wgt
print '(a,x,es10.4)' , 'max skew correction of deformation:',&
print '(a,1x,es10.4)' , 'max skew correction of deformation:',&
maxCorrectionSkew*wgt
print '(a,x,es10.4)' , 'max sym/skew of avg correction: ',&
print '(a,1x,es10.4)' , 'max sym/skew of avg correction: ',&
maxval(math_symmetric33(temp33_real))/&
maxval(math_skew33(temp33_real))
endif
@ -1040,7 +1039,8 @@ program DAMASK_spectral
!--------------------------------------------------------------------------------------------------
! stress BC handling
pstress_av = math_rotate_forward33(pstress_av_lab,bc(loadcase)%rotation)
print '(a,/,3(3(f12.7,x)/)$)', 'Piola-Kirchhoff stress / MPa:',math_transpose33(pstress_av)/1.e6
write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') ' Piola-Kirchhoff stress / MPa:',&
math_transpose33(pstress_av)/1.e6
if(size_reduced > 0_pInt) then ! calculate stress BC if applied
err_stress = maxval(abs(mask_stress * (pstress_av - bc(loadcase)%stress))) ! maximum deviaton (tensor norm not applicable)
@ -1050,8 +1050,9 @@ program DAMASK_spectral
print '(a,es10.4,a,f6.2)', 'error stress = ',err_stress, ', rel. error = ', err_stress/err_stress_tol
defgradAimCorr = - math_mul3333xx33(s_prev, ((pstress_av - bc(loadcase)%stress))) ! residual on given stress components
defgradAim = defgradAim + defgradAimCorr
print '(a,/,3(3(f12.7,x)/)$)', 'new deformation aim: ', math_transpose33(defgradAim)
print '(a,x,es10.4)' , 'with determinant: ', math_det33(defgradAim)
write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') ' new deformation aim: ',&
math_transpose33(defgradAim)
print '(a,1x,es10.4)' , 'with determinant: ', math_det33(defgradAim)
else
err_stress_tol = 0.0_pReal
endif
@ -1062,7 +1063,7 @@ program DAMASK_spectral
defgrad(1:res(1),1:res(2),1:res(3),m,n) + (defgradAim_lab(m,n) - defgrad_av_lab(m,n)) ! anticipated target minus current state
enddo; enddo
if(debugGeneral) then
print '(a,/,3(3(f12.7,x)/)$)', 'average deformation gradient correction:',&
write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') ' average deformation gradient correction:',&
math_transpose33(defgradAim_lab- defgrad_av_lab)
!--------------------------------------------------------------------------------------------------
@ -1076,13 +1077,12 @@ program DAMASK_spectral
defgradDetMin = min(defgradDetMin,defgradDet)
enddo; enddo; enddo
print '(a,x,es10.4)' , 'max determinant of deformation:', defgradDetMax
print '(a,x,es10.4)' , 'min determinant of deformation:', defgradDetMin
print '(a,1x,es10.4)' , 'max determinant of deformation:', defgradDetMax
print '(a,1x,es10.4)' , 'min determinant of deformation:', defgradDetMin
endif
enddo ! end looping when convergency is achieved
!$OMP CRITICAL (write2out)
print '(a)', ''
print '(a)', '============================================================='
if(err_div > err_div_tol .or. err_stress > err_stress_tol) then
@ -1094,25 +1094,22 @@ program DAMASK_spectral
if (mod(totalIncsCounter -1_pInt,bc(loadcase)%outputfrequency) == 0_pInt) then ! at output frequency
print '(a)', ''
print '(a)', '... writing results to file .................................'
write(538), materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! write result to file
write(538) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! write result to file
endif
if (update_gamma) then
print*, 'update c0_reference '
c0_reference = c_current*wgt
endif
!$OMP END CRITICAL (write2out)
endif ! end calculation/forwarding
enddo ! end looping over incs in current loadcase
deallocate(c_reduced)
deallocate(s_reduced)
enddo ! end looping over loadcases
!$OMP CRITICAL (write2out)
print '(a)', ''
print '(a)', '#############################################################'
print '(i6.6,a,i6.6,a)', notConvergedCounter, ' out of ', &
totalIncsCounter - restartReadInc, ' increments did not converge!'
!$OMP END CRITICAL (write2out)
close(538)
call fftw_destroy_plan(plan_stress); call fftw_destroy_plan(plan_correction)
if (debugDivergence) call fftw_destroy_plan(plan_divergence)

4
code/DAMASK_spectral_interface.f90
View File

@ -132,10 +132,10 @@ subroutine DAMASK_interface_init()
write(6,*) '<<<+- DAMASK_spectral_interface init -+>>>'
write(6,*) '$Id$'
#include "compilation_info.f90"
write(6,'(a,2(i2.2,a),i4.4)'), ' Date: ',date_and_time_values(3),'/',&
write(6,'(a,2(i2.2,a),i4.4)') ' Date: ',date_and_time_values(3),'/',&
date_and_time_values(2),'/',&
date_and_time_values(1)
write(6,'(a,2(i2.2,a),i2.2)'), ' Time: ',date_and_time_values(5),':',&
write(6,'(a,2(i2.2,a),i2.2)') ' Time: ',date_and_time_values(5),':',&
date_and_time_values(6),':',&
date_and_time_values(7)
write(6,*) 'Host Name: ', trim(hostName)

10
code/IO.f90
View File

@ -1042,7 +1042,7 @@ endfunction
do
read(unit,'(A65536)',end=100) line
pos = IO_stringPos(line,maxNchunks)
if (verify(IO_stringValue(line,pos,1),"0123456789") > 0) then ! a non-int, i.e. set name
if (verify(IO_stringValue(line,pos,1),'0123456789') > 0) then ! a non-int, i.e. set name
do i = 1,lookupMaxN ! loop over known set names
if (IO_stringValue(line,pos,1) == lookupName(i)) then ! found matching name
IO_continousIntValues = lookupMap(:,i) ! return resp. entity list
@ -1087,7 +1087,7 @@ endfunction
do l = 1,count
read(unit,'(A65536)',end=100) line
pos = IO_stringPos(line,maxNchunks)
if (verify(IO_stringValue(line,pos,1),"0123456789") > 0) then ! a non-int, i.e. set names follow on this line
if (verify(IO_stringValue(line,pos,1),'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 j = 1,lookupMaxN ! look thru known set names
if (IO_stringValue(line,pos,i) == lookupName(j)) then ! found matching name
@ -1447,12 +1447,12 @@ endfunction
if (present(e)) then
if (present(i)) then
if (present(g)) then
write(6,'(a12,x,i6,x,a2,x,i2,x,a5,x,i4,a2)') '+ at element',e,'IP',i,'grain',g,' +'
write(6,'(a12,1x,i6,1x,a2,1x,i2,1x,a5,1x,i4,a2)') '+ at element',e,'IP',i,'grain',g,' +'
else
write(6,'(a12,x,i6,x,a2,x,i2,a13)') '+ at element',e,'IP',i,' +'
write(6,'(a12,1x,i6,1x,a2,1x,i2,a13)') '+ at element',e,'IP',i,' +'
endif
else
write(6,'(a12,x,i6,a19)') '+ at element',e,' +'
write(6,'(a12,1x,i6,a19)') '+ at element',e,' +'
endif
endif
write(6,'(a38)') '+------------------------------------+'

2
code/compilation_info.f90
View File

@ -9,7 +9,7 @@
#endif
#endif
#ifdef __INTEL_COMPILER
write(6,'(a,i4.4,a,i8.8)'), ' Compiled with Intel fortran version ', __INTEL_COMPILER,&
write(6,'(a,i4.4,a,i8.8)') ' Compiled with Intel fortran version ', __INTEL_COMPILER,&
', build date ', __INTEL_COMPILER_BUILD_DATE
#endif
write(6,*) 'Compiled on ', __DATE__,' at ',__TIME__

24
code/constitutive.f90
View File

@ -343,19 +343,19 @@ constitutive_maxSizePostResults = maxval(constitutive_sizePostResults)
write(6,*) '$Id$'
#include "compilation_info.f90"
if (debug_verbosity > 0) then
write(6,'(a32,x,7(i8,x))') 'constitutive_state0: ', shape(constitutive_state0)
write(6,'(a32,x,7(i8,x))') 'constitutive_partionedState0: ', shape(constitutive_partionedState0)
write(6,'(a32,x,7(i8,x))') 'constitutive_subState0: ', shape(constitutive_subState0)
write(6,'(a32,x,7(i8,x))') 'constitutive_state: ', shape(constitutive_state)
write(6,'(a32,x,7(i8,x))') 'constitutive_aTolState: ', shape(constitutive_aTolState)
write(6,'(a32,x,7(i8,x))') 'constitutive_dotState: ', shape(constitutive_dotState)
write(6,'(a32,x,7(i8,x))') 'constitutive_sizeState: ', shape(constitutive_sizeState)
write(6,'(a32,x,7(i8,x))') 'constitutive_sizeDotState: ', shape(constitutive_sizeDotState)
write(6,'(a32,x,7(i8,x))') 'constitutive_sizePostResults: ', shape(constitutive_sizePostResults)
write(6,'(a32,1x,7(i8,1x))') 'constitutive_state0: ', shape(constitutive_state0)
write(6,'(a32,1x,7(i8,1x))') 'constitutive_partionedState0: ', shape(constitutive_partionedState0)
write(6,'(a32,1x,7(i8,1x))') 'constitutive_subState0: ', shape(constitutive_subState0)
write(6,'(a32,1x,7(i8,1x))') 'constitutive_state: ', shape(constitutive_state)
write(6,'(a32,1x,7(i8,1x))') 'constitutive_aTolState: ', shape(constitutive_aTolState)
write(6,'(a32,1x,7(i8,1x))') 'constitutive_dotState: ', shape(constitutive_dotState)
write(6,'(a32,1x,7(i8,1x))') 'constitutive_sizeState: ', shape(constitutive_sizeState)
write(6,'(a32,1x,7(i8,1x))') 'constitutive_sizeDotState: ', shape(constitutive_sizeDotState)
write(6,'(a32,1x,7(i8,1x))') 'constitutive_sizePostResults: ', shape(constitutive_sizePostResults)
write(6,*)
write(6,'(a32,x,7(i8,x))') 'maxSizeState: ', constitutive_maxSizeState
write(6,'(a32,x,7(i8,x))') 'maxSizeDotState: ', constitutive_maxSizeDotState
write(6,'(a32,x,7(i8,x))') 'maxSizePostResults: ', constitutive_maxSizePostResults
write(6,'(a32,1x,7(i8,1x))') 'maxSizeState: ', constitutive_maxSizeState
write(6,'(a32,1x,7(i8,1x))') 'maxSizeDotState: ', constitutive_maxSizeDotState
write(6,'(a32,1x,7(i8,1x))') 'maxSizePostResults: ', constitutive_maxSizePostResults
endif
call flush(6)
!$OMP END CRITICAL (write2out)

2
code/constitutive_j2.f90
View File

@ -107,7 +107,7 @@ subroutine constitutive_j2_init(file)
if (debug_verbosity > 0) then
!$OMP CRITICAL (write2out)
write(6,'(a16,x,i5)') '# instances:',maxNinstance
write(6,'(a16,1x,i5)') '# instances:',maxNinstance
write(6,*)
!$OMP END CRITICAL (write2out)
endif

2
code/constitutive_nonlocal.f90
View File

@ -216,7 +216,7 @@ if (maxNinstance == 0) return
if (debug_verbosity > 0) then
!$OMP CRITICAL (write2out)
write(6,'(a16,x,i5)') '# instances:',maxNinstance
write(6,'(a16,1x,i5)') '# instances:',maxNinstance
!$OMP END CRITICAL (write2out)
endif

2
code/constitutive_phenopowerlaw.f90
View File

@ -182,7 +182,7 @@ subroutine constitutive_phenopowerlaw_init(file)
if (debug_verbosity > 0) then
!$OMP CRITICAL (write2out)
write(6,'(a16,x,i5)') '# instances:',maxNinstance
write(6,'(a16,1x,i5)') '# instances:',maxNinstance
write(6,*)
!$OMP END CRITICAL (write2out)
endif

200
code/crystallite.f90
View File

@ -386,50 +386,50 @@ crystallite_fallbackdPdF = crystallite_dPdF ! use initial ela
! *** Output to MARC output file ***
if (debug_verbosity > 0) then
!$OMP CRITICAL (write2out)
write(6,'(a35,x,7(i8,x))') 'crystallite_Temperature: ', shape(crystallite_Temperature)
write(6,'(a35,x,7(i8,x))') 'crystallite_dotTemperature: ', shape(crystallite_dotTemperature)
write(6,'(a35,x,7(i8,x))') 'crystallite_Fe: ', shape(crystallite_Fe)
write(6,'(a35,x,7(i8,x))') 'crystallite_Fp: ', shape(crystallite_Fp)
write(6,'(a35,x,7(i8,x))') 'crystallite_Lp: ', shape(crystallite_Lp)
write(6,'(a35,x,7(i8,x))') 'crystallite_F0: ', shape(crystallite_F0)
write(6,'(a35,x,7(i8,x))') 'crystallite_Fp0: ', shape(crystallite_Fp0)
write(6,'(a35,x,7(i8,x))') 'crystallite_Lp0: ', shape(crystallite_Lp0)
write(6,'(a35,x,7(i8,x))') 'crystallite_partionedF: ', shape(crystallite_partionedF)
write(6,'(a35,x,7(i8,x))') 'crystallite_partionedTemp0: ', shape(crystallite_partionedTemperature0)
write(6,'(a35,x,7(i8,x))') 'crystallite_partionedF0: ', shape(crystallite_partionedF0)
write(6,'(a35,x,7(i8,x))') 'crystallite_partionedFp0: ', shape(crystallite_partionedFp0)
write(6,'(a35,x,7(i8,x))') 'crystallite_partionedLp0: ', shape(crystallite_partionedLp0)
write(6,'(a35,x,7(i8,x))') 'crystallite_subF: ', shape(crystallite_subF)
write(6,'(a35,x,7(i8,x))') 'crystallite_subTemperature0: ', shape(crystallite_subTemperature0)
write(6,'(a35,x,7(i8,x))') 'crystallite_symmetryID: ', shape(crystallite_symmetryID)
write(6,'(a35,x,7(i8,x))') 'crystallite_subF0: ', shape(crystallite_subF0)
write(6,'(a35,x,7(i8,x))') 'crystallite_subFe0: ', shape(crystallite_subFe0)
write(6,'(a35,x,7(i8,x))') 'crystallite_subFp0: ', shape(crystallite_subFp0)
write(6,'(a35,x,7(i8,x))') 'crystallite_subLp0: ', shape(crystallite_subLp0)
write(6,'(a35,x,7(i8,x))') 'crystallite_P: ', shape(crystallite_P)
write(6,'(a35,x,7(i8,x))') 'crystallite_Tstar_v: ', shape(crystallite_Tstar_v)
write(6,'(a35,x,7(i8,x))') 'crystallite_Tstar0_v: ', shape(crystallite_Tstar0_v)
write(6,'(a35,x,7(i8,x))') 'crystallite_partionedTstar0_v: ', shape(crystallite_partionedTstar0_v)
write(6,'(a35,x,7(i8,x))') 'crystallite_subTstar0_v: ', shape(crystallite_subTstar0_v)
write(6,'(a35,x,7(i8,x))') 'crystallite_dPdF: ', shape(crystallite_dPdF)
write(6,'(a35,x,7(i8,x))') 'crystallite_dPdF0: ', shape(crystallite_dPdF0)
write(6,'(a35,x,7(i8,x))') 'crystallite_partioneddPdF0: ', shape(crystallite_partioneddPdF0)
write(6,'(a35,x,7(i8,x))') 'crystallite_fallbackdPdF: ', shape(crystallite_fallbackdPdF)
write(6,'(a35,x,7(i8,x))') 'crystallite_orientation: ', shape(crystallite_orientation)
write(6,'(a35,x,7(i8,x))') 'crystallite_orientation0: ', shape(crystallite_orientation0)
write(6,'(a35,x,7(i8,x))') 'crystallite_rotation: ', shape(crystallite_rotation)
write(6,'(a35,x,7(i8,x))') 'crystallite_disorientation: ', shape(crystallite_disorientation)
write(6,'(a35,x,7(i8,x))') 'crystallite_dt: ', shape(crystallite_dt)
write(6,'(a35,x,7(i8,x))') 'crystallite_subdt: ', shape(crystallite_subdt)
write(6,'(a35,x,7(i8,x))') 'crystallite_subFrac: ', shape(crystallite_subFrac)
write(6,'(a35,x,7(i8,x))') 'crystallite_subStep: ', shape(crystallite_subStep)
write(6,'(a35,x,7(i8,x))') 'crystallite_stateDamper: ', shape(crystallite_stateDamper)
write(6,'(a35,x,7(i8,x))') 'crystallite_localConstitution: ', shape(crystallite_localConstitution)
write(6,'(a35,x,7(i8,x))') 'crystallite_requested: ', shape(crystallite_requested)
write(6,'(a35,x,7(i8,x))') 'crystallite_todo: ', shape(crystallite_todo)
write(6,'(a35,x,7(i8,x))') 'crystallite_converged: ', shape(crystallite_converged)
write(6,'(a35,x,7(i8,x))') 'crystallite_sizePostResults: ', shape(crystallite_sizePostResults)
write(6,'(a35,x,7(i8,x))') 'crystallite_sizePostResult: ', shape(crystallite_sizePostResult)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Temperature: ', shape(crystallite_Temperature)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_dotTemperature: ', shape(crystallite_dotTemperature)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fe: ', shape(crystallite_Fe)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fp: ', shape(crystallite_Fp)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Lp: ', shape(crystallite_Lp)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_F0: ', shape(crystallite_F0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fp0: ', shape(crystallite_Fp0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Lp0: ', shape(crystallite_Lp0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedF: ', shape(crystallite_partionedF)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedTemp0: ', shape(crystallite_partionedTemperature0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedF0: ', shape(crystallite_partionedF0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedFp0: ', shape(crystallite_partionedFp0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedLp0: ', shape(crystallite_partionedLp0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subF: ', shape(crystallite_subF)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subTemperature0: ', shape(crystallite_subTemperature0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_symmetryID: ', shape(crystallite_symmetryID)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subF0: ', shape(crystallite_subF0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subFe0: ', shape(crystallite_subFe0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subFp0: ', shape(crystallite_subFp0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subLp0: ', shape(crystallite_subLp0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_P: ', shape(crystallite_P)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Tstar_v: ', shape(crystallite_Tstar_v)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Tstar0_v: ', shape(crystallite_Tstar0_v)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedTstar0_v: ', shape(crystallite_partionedTstar0_v)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subTstar0_v: ', shape(crystallite_subTstar0_v)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_dPdF: ', shape(crystallite_dPdF)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_dPdF0: ', shape(crystallite_dPdF0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_partioneddPdF0: ', shape(crystallite_partioneddPdF0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_fallbackdPdF: ', shape(crystallite_fallbackdPdF)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_orientation: ', shape(crystallite_orientation)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_orientation0: ', shape(crystallite_orientation0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_rotation: ', shape(crystallite_rotation)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_disorientation: ', shape(crystallite_disorientation)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_dt: ', shape(crystallite_dt)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subdt: ', shape(crystallite_subdt)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subFrac: ', shape(crystallite_subFrac)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subStep: ', shape(crystallite_subStep)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_stateDamper: ', shape(crystallite_stateDamper)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_localConstitution: ', shape(crystallite_localConstitution)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_requested: ', shape(crystallite_requested)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_todo: ', shape(crystallite_todo)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_converged: ', shape(crystallite_converged)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_sizePostResults: ', shape(crystallite_sizePostResults)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_sizePostResult: ', shape(crystallite_sizePostResult)
write(6,*)
write(6,*) 'Number of nonlocal grains: ',count(.not. crystallite_localConstitution)
call flush(6)
@ -538,13 +538,13 @@ if (debug_verbosity > 4 .and. debug_e > 0 .and. debug_e <= mesh_NcpElems &
.and. debug_g > 0 .and. debug_g <= homogenization_maxNgrains) then
!$OMP CRITICAL (write2out)
write (6,*)
write (6,'(a,i8,x,i2,x,i3)') '<< CRYST >> crystallite start at el ip g ', debug_e, debug_i, debug_g
write (6,'(a,/,12(x),f14.9)') '<< CRYST >> Temp0', crystallite_partionedTemperature0(debug_g,debug_i,debug_e)
write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> F0 ', &
write (6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> crystallite start at el ip g ', debug_e, debug_i, debug_g
write (6,'(a,/,12x,f14.9)') '<< CRYST >> Temp0', crystallite_partionedTemperature0(debug_g,debug_i,debug_e)
write (6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> F0 ', &
math_transpose33(crystallite_partionedF0(1:3,1:3,debug_g,debug_i,debug_e))
write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> Fp0', &
write (6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Fp0', &
math_transpose33(crystallite_partionedFp0(1:3,1:3,debug_g,debug_i,debug_e))
write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> Lp0', &
write (6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Lp0', &
math_transpose33(crystallite_partionedLp0(1:3,1:3,debug_g,debug_i,debug_e))
!$OMP END CRITICAL (write2out)
endif
@ -706,11 +706,11 @@ enddo
#ifndef _OPENMP
if (debug_verbosity > 4 &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
write (6,'(a,i8,x,i2,x,i3)') '<< CRYST >> central solution of cryst_StressAndTangent at el ip g ',e,i,g
write (6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> central solution of cryst_StressAndTangent at el ip g ',e,i,g
write (6,*)
write (6,'(a,/,3(12(x),3(f12.4,x)/))') '<< CRYST >> P / MPa', math_transpose33(crystallite_P(1:3,1:3,g,i,e)) / 1e6
write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> Fp', math_transpose33(crystallite_Fp(1:3,1:3,g,i,e))
write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> Lp', math_transpose33(crystallite_Lp(1:3,1:3,g,i,e))
write (6,'(a,/,3(12x,3(f12.4,1x)/))') '<< CRYST >> P / MPa', math_transpose33(crystallite_P(1:3,1:3,g,i,e)) / 1e6
write (6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Fp', math_transpose33(crystallite_Fp(1:3,1:3,g,i,e))
write (6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Lp', math_transpose33(crystallite_Lp(1:3,1:3,g,i,e))
write (6,*)
endif
#endif
@ -763,7 +763,7 @@ if(updateJaco) then
#ifndef _OPENMP
if (debug_verbosity> 5) then
!$OMP CRITICAL (write2out)
write(6,'(a,2(x,i1),x,a)') '<< CRYST >> [[[[[[ Stiffness perturbation',k,l,']]]]]]'
write(6,'(a,2(x,i1),1x,a)') '<< CRYST >> [[[[[[ Stiffness perturbation',k,l,']]]]]]'
write(6,*)
!$OMP END CRITICAL (write2out)
endif
@ -1082,11 +1082,11 @@ do n = 1,4
if (debug_verbosity > 5 &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
mySizeDotState = constitutive_sizeDotState(g,i,e)
write(6,'(a,i8,x,i2,x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState)
write(6,'(a,/,(12x,12(e12.5,1x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState)
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState)
write(6,'(a,/,(12x,12(e12.5,1x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState)
write(6,*)
endif
#endif
@ -1310,7 +1310,7 @@ endif
! --- FIRST RUNGE KUTTA STEP ---
#ifndef _OPENMP
if (debug_verbosity > 5) then
write(6,'(a,x,i1)') '<< CRYST >> RUNGE KUTTA STEP',1
write(6,'(a,1x,i1)') '<< CRYST >> RUNGE KUTTA STEP',1
endif
#endif
!$OMP DO
@ -1446,7 +1446,7 @@ do n = 1,5
! --- dot state and RK dot state---
#ifndef _OPENMP
if (debug_verbosity > 5) then
write(6,'(a,x,i1)') '<< CRYST >> RUNGE KUTTA STEP',n+1
write(6,'(a,1x,i1)') '<< CRYST >> RUNGE KUTTA STEP',n+1
endif
#endif
!$OMP DO
@ -1571,17 +1571,17 @@ relTemperatureResiduum = 0.0_pReal
#ifndef _OPENMP
if (debug_verbosity > 5 &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
write(6,'(a,i8,x,i3,x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g
write(6,'(a,i8,1x,i3,1x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g
write(6,*)
write(6,'(a,/,(12(x),12(f12.1,x)))') '<< CRYST >> absolute residuum tolerance', &
write(6,'(a,/,(12x,12(f12.1,1x)))') '<< CRYST >> absolute residuum tolerance', &
stateResiduum(1:mySizeDotState,g,i,e) / constitutive_aTolState(g,i,e)%p(1:mySizeDotState)
write(6,*)
write(6,'(a,/,(12(x),12(f12.1,x)))') '<< CRYST >> relative residuum tolerance', &
write(6,'(a,/,(12x,12(f12.1,1x)))') '<< CRYST >> relative residuum tolerance', &
relStateResiduum(1:mySizeDotState,g,i,e) / rTol_crystalliteState
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState)
write(6,'(a,/,(12x,12(e12.5,1x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState)
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState)
write(6,'(a,/,(12x,12(e12.5,1x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState)
write(6,*)
endif
#endif
@ -1883,18 +1883,18 @@ relTemperatureResiduum = 0.0_pReal
#ifndef _OPENMP
if (debug_verbosity > 5 &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
write(6,'(a,i8,x,i2,x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g
write(6,*)
write(6,'(a,/,(12(x),12(f12.1,x)))') '<< CRYST >> absolute residuum tolerance', &
write(6,'(a,/,(12x,12(f12.1,1x)))') '<< CRYST >> absolute residuum tolerance', &
stateResiduum(1:mySizeDotState,g,i,e) / constitutive_aTolState(g,i,e)%p(1:mySizeDotState)
write(6,*)
write(6,'(a,/,(12(x),12(f12.1,x)))') '<< CRYST >> relative residuum tolerance', &
write(6,'(a,/,(12x,12(f12.1,1x)))') '<< CRYST >> relative residuum tolerance', &
relStateResiduum(1:mySizeDotState,g,i,e) / rTol_crystalliteState
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState) &
write(6,'(a,/,(12x,12(e12.5,1x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState) &
- 2.0_pReal * stateResiduum(1:mySizeDotState,g,i,e) / crystallite_subdt(g,i,e) ! calculate former dotstate from higher order solution and state residuum
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState)
write(6,'(a,/,(12x,12(e12.5,1x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState)
write(6,*)
endif
#endif
@ -2057,11 +2057,11 @@ if (numerics_integrationMode < 2) then
#ifndef _OPENMP
if (debug_verbosity > 5 &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
write(6,'(a,i8,x,i2,x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState)
write(6,'(a,/,(12x,12(e12.5,1x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState)
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState)
write(6,'(a,/,(12x,12(e12.5,1x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState)
write(6,*)
endif
#endif
@ -2463,7 +2463,7 @@ residuum = constitutive_state(g,i,e)%p(1:mySize) - constitutive_subState0(g,i,e)
if (any(residuum /= residuum)) then ! if NaN occured then return without changing the state
#ifndef _OPENMP
if (debug_verbosity > 4) then
write(6,'(a,i8,x,i2,x,i3)') '<< CRYST >> updateState encountered NaN at el ip g ',e,i,g
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> updateState encountered NaN at el ip g ',e,i,g
endif
#endif
return
@ -2480,16 +2480,16 @@ converged = all( abs(residuum) < constitutive_aTolState(g,i,e)%p(1:mySize) &
#ifndef _OPENMP
if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
if (converged) then
write(6,'(a,i8,x,i2,x,i3)') '<< CRYST >> updateState converged at el ip g ',e,i,g
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> updateState converged at el ip g ',e,i,g
else
write(6,'(a,i8,x,i2,x,i3)') '<< CRYST >> updateState did not converge at el ip g ',e,i,g
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> updateState did not converge at el ip g ',e,i,g
endif
write(6,*)
write(6,'(a,f6.1)') '<< CRYST >> crystallite_statedamper ',crystallite_statedamper(g,i,e)
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> dotState',dotState(1:mySize)
write(6,'(a,/,(12x,12(e12.5,1x)))') '<< CRYST >> dotState',dotState(1:mySize)
write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> new state',state(1:mySize)
write(6,'(a,/,(12x,12(e12.5,1x)))') '<< CRYST >> new state',state(1:mySize)
write(6,*)
endif
#endif
@ -2546,7 +2546,7 @@ residuum = crystallite_Temperature(g,i,e) - crystallite_subTemperature0(g,i,e) &
if (residuum /= residuum) then
#ifndef _OPENMP
if (debug_verbosity > 4) then
write(6,'(a,i8,x,i2,x,i3)') '<< CRYST >> updateTemperature encountered NaN at el ip g ',e,i,g
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> updateTemperature encountered NaN at el ip g ',e,i,g
endif
#endif
return
@ -2681,7 +2681,7 @@ integer(pLongInt) tick, &
crystallite_integrateStress = .false.
#ifndef _OPENMP
if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
write(6,'(a,i8,x,i2,x,i3)') '<< CRYST >> integrateStress at el ip g ',e,i,g
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> integrateStress at el ip g ',e,i,g
endif
#endif
@ -2710,10 +2710,10 @@ invFp_current = math_inv33(Fp_current)
if (all(invFp_current == 0.0_pReal)) then ! ... failed?
#ifndef _OPENMP
if (debug_verbosity > 4) then
write(6,'(a,i8,x,i2,x,i3)') '<< CRYST >> integrateStress failed on invFp_current inversion at el ip g ',e,i,g
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on invFp_current inversion at el ip g ',e,i,g
if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
write(6,*)
write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> invFp_new',math_transpose33(invFp_new(1:3,1:3))
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> invFp_new',math_transpose33(invFp_new(1:3,1:3))
endif
endif
#endif
@ -2745,7 +2745,7 @@ LpLoop: do
if (NiterationStress > nStress) then
#ifndef _OPENMP
if (debug_verbosity > 4) then
write(6,'(a,i8,x,i2,x,i3)') '<< CRYST >> integrateStress reached loop limit at el ip g ',e,i,g
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> integrateStress reached loop limit at el ip g ',e,i,g
write(6,*)
endif
#endif
@ -2786,8 +2786,8 @@ LpLoop: do
.and. numerics_integrationMode == 1_pInt) then
write(6,'(a,i3)') '<< CRYST >> iteration ', NiterationStress
write(6,*)
write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> Lp_constitutive', math_transpose33(Lp_constitutive)
write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> Lpguess', math_transpose33(Lpguess)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lp_constitutive', math_transpose33(Lp_constitutive)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lpguess', math_transpose33(Lpguess)
endif
#endif
@ -2810,7 +2810,7 @@ LpLoop: do
if (steplength >= steplength0 .and. any(residuum /= residuum)) then
#ifndef _OPENMP
if (debug_verbosity > 4) then
write(6,'(a,i8,x,i2,x,i3,a,i3,a)') '<< CRYST >> integrateStress encountered NaN at el ip g ',e,i,g,&
write(6,'(a,i8,1x,i2,1x,i3,a,i3,a)') '<< CRYST >> integrateStress encountered NaN at el ip g ',e,i,g,&
' ; iteration ', NiterationStress,&
' >> returning..!'
endif
@ -2858,7 +2858,7 @@ LpLoop: do
else
#ifndef _OPENMP
if (debug_verbosity > 5) then
write(6,'(a,i8,x,i2,x,i3,x,a,i3)') '<< CRYST >> integrateStress encountered high-speed crash at el ip g ',e,i,g,&
write(6,'(a,i8,1x,i2,1x,i3,1x,a,i3)') '<< CRYST >> integrateStress encountered high-speed crash at el ip g ',e,i,g,&
'; iteration ', NiterationStress
endif
#endif
@ -2895,17 +2895,17 @@ LpLoop: do
if (error) then
#ifndef _OPENMP
if (debug_verbosity > 4) then
write(6,'(a,i8,x,i2,x,i3,a,i3)') '<< CRYST >> integrateStress failed on dR/dLp inversion at el ip g ',e,i,g
write(6,'(a,i8,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on dR/dLp inversion at el ip g ',e,i,g
if (debug_verbosity > 5 &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
write(6,*)
write(6,'(a,/,9(12(x),9(e15.3,x)/))') '<< CRYST >> dR_dLp',transpose(dR_dLp)
write(6,'(a,/,9(12(x),9(e15.3,x)/))') '<< CRYST >> dT_dLp',transpose(dT_dLp)
write(6,'(a,/,9(12(x),9(e15.3,x)/))') '<< CRYST >> dLp_dT_constitutive',transpose(dLp_dT_constitutive)
write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> AB',math_transpose33(AB)
write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> BTA',math_transpose33(BTA)
write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> Lp_constitutive',math_transpose33(Lp_constitutive)
write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> Lpguess',math_transpose33(Lpguess)
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dR_dLp',transpose(dR_dLp)
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dT_dLp',transpose(dT_dLp)
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dLp_dT_constitutive',transpose(dLp_dT_constitutive)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> AB',math_transpose33(AB)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> BTA',math_transpose33(BTA)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lp_constitutive',math_transpose33(Lp_constitutive)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lpguess',math_transpose33(Lpguess)
endif
endif
#endif
@ -2940,11 +2940,11 @@ call math_invert33(invFp_new,Fp_new,det,error)
if (error) then
#ifndef _OPENMP
if (debug_verbosity > 4) then
write(6,'(a,i8,x,i2,x,i3,a,i3)') '<< CRYST >> integrateStress failed on invFp_new inversion at el ip g ',e,i,g, &
write(6,'(a,i8,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on invFp_new inversion at el ip g ',e,i,g, &
' ; iteration ', NiterationStress
if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
write(6,*)
write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> invFp_new',math_transpose33(invFp_new)
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> invFp_new',math_transpose33(invFp_new)
endif
endif
#endif
@ -2974,12 +2974,12 @@ crystallite_integrateStress = .true.
#ifndef _OPENMP
if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger) &
.and. numerics_integrationMode == 1_pInt) then
write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> P / MPa',math_transpose33(crystallite_P(1:3,1:3,g,i,e))/1e6
write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> Cauchy / MPa', &
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> P / MPa',math_transpose33(crystallite_P(1:3,1:3,g,i,e))/1e6
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Cauchy / MPa', &
math_mul33x33(crystallite_P(1:3,1:3,g,i,e), math_transpose33(Fg_new)) / 1e6 / math_det33(Fg_new)
write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> Fe Lp Fe^-1', &
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fe Lp Fe^-1', &
math_transpose33(math_mul33x33(Fe_new, math_mul33x33(crystallite_Lp(1:3,1:3,g,i,e), math_inv33(Fe_new)))) ! transpose to get correct print out order
write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> Fp',math_transpose33(crystallite_Fp(1:3,1:3,g,i,e))
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fp',math_transpose33(crystallite_Fp(1:3,1:3,g,i,e))
endif
#endif

56
code/debug.f90
View File

@ -163,16 +163,16 @@ subroutine debug_init()
if (debug_verbosity > 0) then
!$OMP CRITICAL (write2out)
write(6,'(a24,x,i1)') 'verbose: ',debug_verbosity
write(6,'(a24,x,l)') 'selective: ',debug_selectiveDebugger
write(6,'(a24,1x,i1)') 'verbose: ',debug_verbosity
write(6,'(a24,1x,l1)') 'selective: ',debug_selectiveDebugger
!$OMP END CRITICAL (write2out)
endif
if (debug_selectiveDebugger) then
if (debug_verbosity > 0) then
!$OMP CRITICAL (write2out)
write(6,'(a24,x,i8)') 'element: ',debug_e
write(6,'(a24,x,i8)') 'ip: ',debug_i
write(6,'(a24,x,i8)') 'grain: ',debug_g
write(6,'(a24,1x,i8)') 'element: ',debug_e
write(6,'(a24,1x,i8)') 'ip: ',debug_i
write(6,'(a24,1x,i8)') 'grain: ',debug_g
!$OMP END CRITICAL (write2out)
endif
else
@ -245,24 +245,24 @@ subroutine debug_info()
write(6,*)
write(6,*) 'DEBUG Info (from previous cycle)'
write(6,*)
write(6,'(a33,x,i12)') 'total calls to LpAndItsTangent :',debug_cumLpCalls
write(6,'(a33,1x,i12)') 'total calls to LpAndItsTangent :',debug_cumLpCalls
if (debug_cumLpCalls > 0_pInt) then
write(6,'(a33,x,f12.3)') 'total CPU time/s :',dble(debug_cumLpTicks)/tickrate
write(6,'(a33,x,f12.6)') 'avg CPU time/microsecs per call :',&
write(6,'(a33,1x,f12.3)') 'total CPU time/s :',dble(debug_cumLpTicks)/tickrate
write(6,'(a33,1x,f12.6)') 'avg CPU time/microsecs per call :',&
dble(debug_cumLpTicks)*1.0e6_pReal/tickrate/debug_cumLpCalls
endif
write(6,*)
write(6,'(a33,x,i12)') 'total calls to collectDotState :',debug_cumDotStateCalls
write(6,'(a33,1x,i12)') 'total calls to collectDotState :',debug_cumDotStateCalls
if (debug_cumdotStateCalls > 0_pInt) then
write(6,'(a33,x,f12.3)') 'total CPU time/s :',dble(debug_cumDotStateTicks)/tickrate
write(6,'(a33,x,f12.6)') 'avg CPU time/microsecs per call :',&
write(6,'(a33,1x,f12.3)') 'total CPU time/s :',dble(debug_cumDotStateTicks)/tickrate
write(6,'(a33,1x,f12.6)') 'avg CPU time/microsecs per call :',&
dble(debug_cumDotStateTicks)*1.0e6_pReal/tickrate/debug_cumDotStateCalls
endif
write(6,*)
write(6,'(a33,x,i12)') 'total calls to dotTemperature :',debug_cumDotTemperatureCalls
write(6,'(a33,1x,i12)') 'total calls to dotTemperature :',debug_cumDotTemperatureCalls
if (debug_cumdotTemperatureCalls > 0_pInt) then
write(6,'(a33,x,f12.3)') 'total CPU time/s :', dble(debug_cumDotTemperatureTicks)/tickrate
write(6,'(a33,x,f12.6)') 'avg CPU time/microsecs per call :',&
write(6,'(a33,1x,f12.3)') 'total CPU time/s :', dble(debug_cumDotTemperatureTicks)/tickrate
write(6,'(a33,1x,f12.6)') 'avg CPU time/microsecs per call :',&
dble(debug_cumDotTemperatureTicks)*1.0e6_pReal/tickrate/debug_cumDotTemperatureCalls
endif
@ -274,11 +274,11 @@ subroutine debug_info()
if (any(debug_StressLoopDistribution(i,:) /= 0_pInt ) .or. &
any(debug_LeapfrogBreakDistribution(i,:) /= 0_pInt ) ) then
integral = integral + i*debug_StressLoopDistribution(i,1) + i*debug_StressLoopDistribution(i,2)
write(6,'(i25,x,i10,x,i10,x,i10,x,i10)') i,debug_StressLoopDistribution(i,1),debug_LeapfrogBreakDistribution(i,1), &
write(6,'(i25,1x,i10,1x,i10,1x,i10,1x,i10)') i,debug_StressLoopDistribution(i,1),debug_LeapfrogBreakDistribution(i,1), &
debug_StressLoopDistribution(i,2),debug_LeapfrogBreakDistribution(i,2)
endif
enddo
write(6,'(a15,i10,x,i10,12x,i10)') ' total',integral,&
write(6,'(a15,i10,1x,i10,12x,i10)') ' total',integral,&
sum(debug_StressLoopDistribution(:,1)), &
sum(debug_StressLoopDistribution(:,2))
@ -288,10 +288,10 @@ subroutine debug_info()
do i=1,nState
if (any(debug_StateLoopDistribution(i,:) /= 0)) then
integral = integral + i*debug_StateLoopDistribution(i,1) + i*debug_StateLoopDistribution(i,2)
write(6,'(i25,x,i10,12x,i10)') i,debug_StateLoopDistribution(i,1),debug_StateLoopDistribution(i,2)
write(6,'(i25,1x,i10,12x,i10)') i,debug_StateLoopDistribution(i,1),debug_StateLoopDistribution(i,2)
endif
enddo
write(6,'(a15,i10,x,i10,12x,i10)') ' total',integral,&
write(6,'(a15,i10,1x,i10,12x,i10)') ' total',integral,&
sum(debug_StateLoopDistribution(:,1)), &
sum(debug_StateLoopDistribution(:,2))
@ -302,13 +302,13 @@ subroutine debug_info()
if (debug_CrystalliteLoopDistribution(i) /= 0) then
integral = integral + i*debug_CrystalliteLoopDistribution(i)
if (i <= nCryst) then
write(6,'(i25,x,i10)') i,debug_CrystalliteLoopDistribution(i)
write(6,'(i25,1x,i10)') i,debug_CrystalliteLoopDistribution(i)
else
write(6,'(i25,a1,i10)') i-1,'+',debug_CrystalliteLoopDistribution(i)
endif
endif
enddo
write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_CrystalliteLoopDistribution)
write(6,'(a15,i10,1x,i10)') ' total',integral,sum(debug_CrystalliteLoopDistribution)
!$OMP END CRITICAL (write2out)
endif
@ -322,10 +322,10 @@ subroutine debug_info()
do i=1,nMPstate
if (debug_MaterialpointStateLoopDistribution(i) /= 0) then
integral = integral + i*debug_MaterialpointStateLoopDistribution(i)
write(6,'(i25,x,i10)') i,debug_MaterialpointStateLoopDistribution(i)
write(6,'(i25,1x,i10)') i,debug_MaterialpointStateLoopDistribution(i)
endif
enddo
write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_MaterialpointStateLoopDistribution)
write(6,'(a15,i10,1x,i10)') ' total',integral,sum(debug_MaterialpointStateLoopDistribution)
integral = 0_pInt
write(6,*)
@ -334,23 +334,23 @@ subroutine debug_info()
if (debug_MaterialpointLoopDistribution(i) /= 0) then
integral = integral + i*debug_MaterialpointLoopDistribution(i)
if (i <= nHomog) then
write(6,'(i25,x,i10)') i,debug_MaterialpointLoopDistribution(i)
write(6,'(i25,1x,i10)') i,debug_MaterialpointLoopDistribution(i)
else
write(6,'(i25,a1,i10)') i-1,'+',debug_MaterialpointLoopDistribution(i)
endif
endif
enddo
write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_MaterialpointLoopDistribution)
write(6,'(a15,i10,1x,i10)') ' total',integral,sum(debug_MaterialpointLoopDistribution)
write(6,*)
write(6,*)
write(6,*) 'Extreme values of returned stress and jacobian'
write(6,*)
write(6,'(a39)') ' value el ip'
write(6,'(a14,x,e12.3,x,i6,x,i4)') 'stress min :', debug_stressMin, debug_stressMinLocation
write(6,'(a14,x,e12.3,x,i6,x,i4)') ' max :', debug_stressMax, debug_stressMaxLocation
write(6,'(a14,x,e12.3,x,i6,x,i4)') 'jacobian min :', debug_jacobianMin, debug_jacobianMinLocation
write(6,'(a14,x,e12.3,x,i6,x,i4)') ' max :', debug_jacobianMax, debug_jacobianMaxLocation
write(6,'(a14,1x,e12.3,1x,i6,1x,i4)') 'stress min :', debug_stressMin, debug_stressMinLocation
write(6,'(a14,1x,e12.3,1x,i6,1x,i4)') ' max :', debug_stressMax, debug_stressMaxLocation
write(6,'(a14,1x,e12.3,1x,i6,1x,i4)') 'jacobian min :', debug_jacobianMin, debug_jacobianMinLocation
write(6,'(a14,1x,e12.3,1x,i6,1x,i4)') ' max :', debug_jacobianMax, debug_jacobianMaxLocation
write(6,*)
!$OMP END CRITICAL (write2out)

54
code/homogenization.f90
View File

@ -209,30 +209,30 @@ allocate(materialpoint_results(materialpoint_sizeResults,mesh_maxNips,mesh_NcpEl
write(6,*) '$Id$'
#include "compilation_info.f90"
if (debug_verbosity > 0) then
write(6,'(a32,x,7(i8,x))') 'homogenization_state0: ', shape(homogenization_state0)
write(6,'(a32,x,7(i8,x))') 'homogenization_subState0: ', shape(homogenization_subState0)
write(6,'(a32,x,7(i8,x))') 'homogenization_state: ', shape(homogenization_state)
write(6,'(a32,x,7(i8,x))') 'homogenization_sizeState: ', shape(homogenization_sizeState)
write(6,'(a32,x,7(i8,x))') 'homogenization_sizePostResults: ', shape(homogenization_sizePostResults)
write(6,'(a32,1x,7(i8,1x))') 'homogenization_state0: ', shape(homogenization_state0)
write(6,'(a32,1x,7(i8,1x))') 'homogenization_subState0: ', shape(homogenization_subState0)
write(6,'(a32,1x,7(i8,1x))') 'homogenization_state: ', shape(homogenization_state)
write(6,'(a32,1x,7(i8,1x))') 'homogenization_sizeState: ', shape(homogenization_sizeState)
write(6,'(a32,1x,7(i8,1x))') 'homogenization_sizePostResults: ', shape(homogenization_sizePostResults)
write(6,*)
write(6,'(a32,x,7(i8,x))') 'materialpoint_dPdF: ', shape(materialpoint_dPdF)
write(6,'(a32,x,7(i8,x))') 'materialpoint_F0: ', shape(materialpoint_F0)
write(6,'(a32,x,7(i8,x))') 'materialpoint_F: ', shape(materialpoint_F)
write(6,'(a32,x,7(i8,x))') 'materialpoint_subF0: ', shape(materialpoint_subF0)
write(6,'(a32,x,7(i8,x))') 'materialpoint_subF: ', shape(materialpoint_subF)
write(6,'(a32,x,7(i8,x))') 'materialpoint_P: ', shape(materialpoint_P)
write(6,'(a32,x,7(i8,x))') 'materialpoint_Temperature: ', shape(materialpoint_Temperature)
write(6,'(a32,x,7(i8,x))') 'materialpoint_subFrac: ', shape(materialpoint_subFrac)
write(6,'(a32,x,7(i8,x))') 'materialpoint_subStep: ', shape(materialpoint_subStep)
write(6,'(a32,x,7(i8,x))') 'materialpoint_subdt: ', shape(materialpoint_subdt)
write(6,'(a32,x,7(i8,x))') 'materialpoint_requested: ', shape(materialpoint_requested)
write(6,'(a32,x,7(i8,x))') 'materialpoint_converged: ', shape(materialpoint_converged)
write(6,'(a32,x,7(i8,x))') 'materialpoint_doneAndHappy: ', shape(materialpoint_doneAndHappy)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_dPdF: ', shape(materialpoint_dPdF)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_F0: ', shape(materialpoint_F0)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_F: ', shape(materialpoint_F)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subF0: ', shape(materialpoint_subF0)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subF: ', shape(materialpoint_subF)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_P: ', shape(materialpoint_P)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_Temperature: ', shape(materialpoint_Temperature)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subFrac: ', shape(materialpoint_subFrac)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subStep: ', shape(materialpoint_subStep)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subdt: ', shape(materialpoint_subdt)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_requested: ', shape(materialpoint_requested)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_converged: ', shape(materialpoint_converged)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_doneAndHappy: ', shape(materialpoint_doneAndHappy)
write(6,*)
write(6,'(a32,x,7(i8,x))') 'materialpoint_results: ', shape(materialpoint_results)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_results: ', shape(materialpoint_results)
write(6,*)
write(6,'(a32,x,7(i8,x))') 'maxSizeState: ', homogenization_maxSizeState
write(6,'(a32,x,7(i8,x))') 'maxSizePostResults: ', homogenization_maxSizePostResults
write(6,'(a32,1x,7(i8,1x))') 'maxSizeState: ', homogenization_maxSizeState
write(6,'(a32,1x,7(i8,1x))') 'maxSizePostResults: ', homogenization_maxSizePostResults
endif
call flush(6)
!$OMP END CRITICAL (write2out)
@ -310,10 +310,10 @@ subroutine materialpoint_stressAndItsTangent(&
if (debug_verbosity > 2 .and. debug_e > 0 .and. debug_e <= mesh_NcpElems .and. debug_i > 0 .and. debug_i <= mesh_maxNips) then
!$OMP CRITICAL (write2out)
write (6,*)
write (6,'(a,i5,x,i2)') '<< HOMOG >> Material Point start at el ip ', debug_e, debug_i
write (6,'(a,/,12(x),f14.9)') '<< HOMOG >> Temp0', materialpoint_Temperature(debug_i,debug_e)
write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< HOMOG >> F0', math_transpose33(materialpoint_F0(1:3,1:3,debug_i,debug_e))
write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< HOMOG >> F', math_transpose33(materialpoint_F(1:3,1:3,debug_i,debug_e))
write (6,'(a,i5,1x,i2)') '<< HOMOG >> Material Point start at el ip ', debug_e, debug_i
write (6,'(a,/,12x,f14.9)') '<< HOMOG >> Temp0', materialpoint_Temperature(debug_i,debug_e)
write (6,'(a,/,3(12x,3(f14.9,1x)/))') '<< HOMOG >> F0', math_transpose33(materialpoint_F0(1:3,1:3,debug_i,debug_e))
write (6,'(a,/,3(12x,3(f14.9,1x)/))') '<< HOMOG >> F', math_transpose33(materialpoint_F(1:3,1:3,debug_i,debug_e))
!$OMP END CRITICAL (write2out)
endif
@ -360,7 +360,7 @@ subroutine materialpoint_stressAndItsTangent(&
if ( materialpoint_converged(i,e) ) then
#ifndef _OPENMP
if (debug_verbosity > 2 .and. ((e == debug_e .and. i == debug_i) .or. .not. debug_selectiveDebugger)) then
write(6,'(a,x,f10.8,x,a,x,f10.8,x,a,/)') '<< HOMOG >> winding forward from', &
write(6,'(a,1x,f10.8,1x,a,1x,f10.8,1x,a,/)') '<< HOMOG >> winding forward from', &
materialpoint_subFrac(i,e), 'to current materialpoint_subFrac', &
materialpoint_subFrac(i,e)+materialpoint_subStep(i,e),'in materialpoint_stressAndItsTangent'
endif
@ -411,7 +411,7 @@ subroutine materialpoint_stressAndItsTangent(&
#ifndef _OPENMP
if (debug_verbosity > 2 .and. ((e == debug_e .and. i == debug_i) .or. .not. debug_selectiveDebugger)) then
write(6,'(a,x,f10.8,/)') &
write(6,'(a,1x,f10.8,/)') &
'<< HOMOG >> cutback step in materialpoint_stressAndItsTangent with new materialpoint_subStep:',&
materialpoint_subStep(i,e)
endif

118
code/homogenization_RGC.f90
View File

@ -172,13 +172,13 @@ subroutine homogenization_RGC_init(&
100 if (debug_verbosity == 4) then
!$OMP CRITICAL (write2out)
do i = 1,maxNinstance
write(6,'(a15,x,i4)') 'instance: ', i
write(6,'(a15,1x,i4)') 'instance: ', i
write(6,*)
write(6,'(a25,3(x,i8))') 'cluster size: ',(homogenization_RGC_Ngrains(j,i),j=1,3)
write(6,'(a25,x,e10.3)') 'scaling parameter: ', homogenization_RGC_xiAlpha(i)
write(6,'(a25,x,e10.3)') 'over-proportionality: ', homogenization_RGC_ciAlpha(i)
write(6,'(a25,3(x,e10.3))') 'grain size: ',(homogenization_RGC_dAlpha(j,i),j=1,3)
write(6,'(a25,3(x,e10.3))') 'cluster orientation: ',(homogenization_RGC_angles(j,i),j=1,3)
write(6,'(a25,3(1x,i8))') 'cluster size: ',(homogenization_RGC_Ngrains(j,i),j=1,3)
write(6,'(a25,1x,e10.3)') 'scaling parameter: ', homogenization_RGC_xiAlpha(i)
write(6,'(a25,1x,e10.3)') 'over-proportionality: ', homogenization_RGC_ciAlpha(i)
write(6,'(a25,3(1x,e10.3))') 'grain size: ',(homogenization_RGC_dAlpha(j,i),j=1,3)
write(6,'(a25,3(1x,e10.3))') 'cluster orientation: ',(homogenization_RGC_angles(j,i),j=1,3)
enddo
!$OMP END CRITICAL (write2out)
endif
@ -278,10 +278,10 @@ subroutine homogenization_RGC_partitionDeformation(&
!* Debugging the overall deformation gradient
if (debug_verbosity == 4) then
!$OMP CRITICAL (write2out)
write(6,'(x,a,i3,a,i3,a)')'========== Increment: ',theInc,' Cycle: ',cycleCounter,' =========='
write(6,'(x,a32)')'Overall deformation gradient: '
write(6,'(1x,a,i3,a,i3,a)')'========== Increment: ',theInc,' Cycle: ',cycleCounter,' =========='
write(6,'(1x,a32)')'Overall deformation gradient: '
do i = 1,3
write(6,'(x,3(e14.8,x))')(avgF(i,j), j = 1,3)
write(6,'(1x,3(e14.8,1x))')(avgF(i,j), j = 1,3)
enddo
write(6,*)' '
call flush(6)
@ -305,9 +305,9 @@ subroutine homogenization_RGC_partitionDeformation(&
!* Debugging the grain deformation gradients
if (debug_verbosity == 4) then
!$OMP CRITICAL (write2out)
write(6,'(x,a32,x,i3)')'Deformation gradient of grain: ',iGrain
write(6,'(1x,a32,1x,i3)')'Deformation gradient of grain: ',iGrain
do i = 1,3
write(6,'(x,3(e14.8,x))')(F(i,j,iGrain), j = 1,3)
write(6,'(1x,3(e14.8,1x))')(F(i,j,iGrain), j = 1,3)
enddo
write(6,*)' '
call flush(6)
@ -392,9 +392,9 @@ function homogenization_RGC_updateState(&
!* Debugging the obtained state
if (debug_verbosity == 4) then
!$OMP CRITICAL (write2out)
write(6,'(x,a30)')'Obtained state: '
write(6,'(1x,a30)')'Obtained state: '
do i = 1,3*nIntFaceTot
write(6,'(x,2(e14.8,x))')state%p(i)
write(6,'(1x,2(e14.8,1x))')state%p(i)
enddo
write(6,*)' '
!$OMP END CRITICAL (write2out)
@ -410,11 +410,11 @@ function homogenization_RGC_updateState(&
if (debug_verbosity == 4) then
!$OMP CRITICAL (write2out)
do iGrain = 1,nGrain
write(6,'(x,a30,x,i3,x,a4,3(x,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,e14.8))')'Mismatch magnitude of grain(',iGrain,') :',NN(1,iGrain),NN(2,iGrain),NN(3,iGrain)
write(6,*)' '
write(6,'(x,a30,x,i3)')'Stress and penalties of grain: ',iGrain
write(6,'(1x,a30,1x,i3)')'Stress and penalties of grain: ',iGrain
do i = 1,3
write(6,'(x,3(e14.8,x),x,3(e14.8,x),x,3(e14.8,x))')(P(i,j,iGrain), j = 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), &
(R(i,j,iGrain), j = 1,3), &
(D(i,j,iGrain), j = 1,3)
enddo
@ -458,8 +458,8 @@ function homogenization_RGC_updateState(&
!* Debugging the residual stress
if (debug_verbosity == 4) then
!$OMP CRITICAL (write2out)
write(6,'(x,a30,x,i3)')'Traction at interface: ',iNum
write(6,'(x,3(e14.8,x))')(tract(iNum,j), j = 1,3)
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,*)' '
!$OMP END CRITICAL (write2out)
endif
@ -476,11 +476,11 @@ function homogenization_RGC_updateState(&
!* Debugging the convergent criteria
if (debug_verbosity == 4 .and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out)
write(6,'(x,a)')' '
write(6,'(x,a,x,i2,x,i4)')'RGC residual check ...',ip,el
write(6,'(x,a15,x,e14.8,x,a7,i3,x,a12,i2,i2)')'Max stress: ',stresMax, &
write(6,'(1x,a)')' '
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, &
'@ grain',stresLoc(3),'in component',stresLoc(1),stresLoc(2)
write(6,'(x,a15,x,e14.8,x,a7,i3,x,a12,i2)')'Max residual: ',residMax, &
write(6,'(1x,a15,1x,e14.8,1x,a7,i3,1x,a12,i2)')'Max residual: ',residMax, &
'@ iface',residLoc(1),'in direction',residLoc(2)
call flush(6)
!$OMP END CRITICAL (write2out)
@ -493,12 +493,12 @@ function homogenization_RGC_updateState(&
if (debug_verbosity == 4 .and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out)
write(6,'(x,a55)')'... done and happy'
write(6,'(1x,a55)')'... done and happy'
write(6,*)' '
call flush(6)
!$OMP END CRITICAL (write2out)
endif
! write(6,'(x,a,x,i3,x,a6,x,i3,x,a12)')'RGC_updateState: ip',ip,'| el',el,'converged :)'
! write(6,'(1x,a,1x,i3,1x,a6,1x,i3,1x,a12)')'RGC_updateState: ip',ip,'| el',el,'converged :)'
!* Then compute/update the state for postResult, i.e., ...
!* ... all energy densities computed by time-integration
@ -523,15 +523,15 @@ function homogenization_RGC_updateState(&
if (debug_verbosity == 4 .and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out)
write(6,'(x,a30,x,e14.8)')'Constitutive work: ',constitutiveWork
write(6,'(x,a30,3(x,e14.8))')'Magnitude mismatch: ',sum(NN(1,:))/dble(nGrain), &
write(6,'(1x,a30,1x,e14.8)')'Constitutive work: ',constitutiveWork
write(6,'(1x,a30,3(1x,e14.8))')'Magnitude mismatch: ',sum(NN(1,:))/dble(nGrain), &
sum(NN(2,:))/dble(nGrain), &
sum(NN(3,:))/dble(nGrain)
write(6,'(x,a30,x,e14.8)')'Penalty energy: ',penaltyEnergy
write(6,'(x,a30,x,e14.8)')'Volume discrepancy: ',volDiscrep
write(6,'(1x,a30,1x,e14.8)')'Penalty energy: ',penaltyEnergy
write(6,'(1x,a30,1x,e14.8)')'Volume discrepancy: ',volDiscrep
write(6,*)''
write(6,'(x,a30,x,e14.8)')'Maximum relaxation rate: ',maxval(abs(drelax))/dt
write(6,'(x,a30,x,e14.8)')'Average relaxation rate: ',sum(abs(drelax))/dt/dble(3*nIntFaceTot)
write(6,'(1x,a30,1x,e14.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,*)''
call flush(6)
!$OMP END CRITICAL (write2out)
@ -547,7 +547,7 @@ function homogenization_RGC_updateState(&
if (debug_verbosity == 4 .and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out)
write(6,'(x,a55)')'... broken'
write(6,'(1x,a55)')'... broken'
write(6,*)' '
call flush(6)
!$OMP END CRITICAL (write2out)
@ -561,7 +561,7 @@ function homogenization_RGC_updateState(&
if (debug_verbosity == 4 .and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out)
write(6,'(x,a55)')'... not yet done'
write(6,'(1x,a55)')'... not yet done'
write(6,*)' '
call flush(6)
!$OMP END CRITICAL (write2out)
@ -617,9 +617,9 @@ function homogenization_RGC_updateState(&
!* Debugging the global Jacobian matrix of stress tangent
if (debug_verbosity == 4) then
!$OMP CRITICAL (write2out)
write(6,'(x,a30)')'Jacobian matrix of stress'
write(6,'(1x,a30)')'Jacobian matrix of stress'
do i = 1,3*nIntFaceTot
write(6,'(x,100(e10.4,x))')(smatrix(i,j), j = 1,3*nIntFaceTot)
write(6,'(1x,100(e10.4,1x))')(smatrix(i,j), j = 1,3*nIntFaceTot)
enddo
write(6,*)' '
call flush(6)
@ -673,9 +673,9 @@ function homogenization_RGC_updateState(&
!* Debugging the global Jacobian matrix of penalty tangent
if (debug_verbosity == 4) then
!$OMP CRITICAL (write2out)
write(6,'(x,a30)')'Jacobian matrix of penalty'
write(6,'(1x,a30)')'Jacobian matrix of penalty'
do i = 1,3*nIntFaceTot
write(6,'(x,100(e10.4,x))')(pmatrix(i,j), j = 1,3*nIntFaceTot)
write(6,'(1x,100(e10.4,1x))')(pmatrix(i,j), j = 1,3*nIntFaceTot)
enddo
write(6,*)' '
call flush(6)
@ -693,9 +693,9 @@ function homogenization_RGC_updateState(&
!* Debugging the global Jacobian matrix of numerical viscosity tangent
if (debug_verbosity == 4) then
!$OMP CRITICAL (write2out)
write(6,'(x,a30)')'Jacobian matrix of penalty'
write(6,'(1x,a30)')'Jacobian matrix of penalty'
do i = 1,3*nIntFaceTot
write(6,'(x,100(e10.4,x))')(rmatrix(i,j), j = 1,3*nIntFaceTot)
write(6,'(1x,100(e10.4,1x))')(rmatrix(i,j), j = 1,3*nIntFaceTot)
enddo
write(6,*)' '
call flush(6)
@ -707,9 +707,9 @@ function homogenization_RGC_updateState(&
if (debug_verbosity == 4) then
!$OMP CRITICAL (write2out)
write(6,'(x,a30)')'Jacobian matrix (total)'
write(6,'(1x,a30)')'Jacobian matrix (total)'
do i = 1,3*nIntFaceTot
write(6,'(x,100(e10.4,x))')(jmatrix(i,j), j = 1,3*nIntFaceTot)
write(6,'(1x,100(e10.4,1x))')(jmatrix(i,j), j = 1,3*nIntFaceTot)
enddo
write(6,*)' '
call flush(6)
@ -726,9 +726,9 @@ function homogenization_RGC_updateState(&
!* Debugging the inverse Jacobian matrix
if (debug_verbosity == 4) then
!$OMP CRITICAL (write2out)
write(6,'(x,a30)')'Jacobian inverse'
write(6,'(1x,a30)')'Jacobian inverse'
do i = 1,3*nIntFaceTot
write(6,'(x,100(e10.4,x))')(jnverse(i,j), j = 1,3*nIntFaceTot)
write(6,'(1x,100(e10.4,1x))')(jnverse(i,j), j = 1,3*nIntFaceTot)
enddo
write(6,*)' '
call flush(6)
@ -747,8 +747,8 @@ function homogenization_RGC_updateState(&
if (any(abs(drelax(:)) > maxdRelax_RGC)) then ! Forcing cutback when the incremental change of relaxation vector becomes too large
homogenization_RGC_updateState = (/.true.,.false./)
!$OMP CRITICAL (write2out)
write(6,'(x,a,x,i3,x,a,x,i3,x,a)')'RGC_updateState: ip',ip,'| el',el,'enforces cutback'
write(6,'(x,a,x,e14.8)')'due to large relaxation change =',maxval(abs(drelax))
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))
call flush(6)
!$OMP END CRITICAL (write2out)
endif
@ -756,9 +756,9 @@ function homogenization_RGC_updateState(&
!* Debugging the return state
if (debug_verbosity == 4) then
!$OMP CRITICAL (write2out)
write(6,'(x,a30)')'Returned state: '
write(6,'(1x,a30)')'Returned state: '
do i = 1,3*nIntFaceTot
write(6,'(x,2(e14.8,x))')state%p(i)
write(6,'(1x,2(e14.8,1x))')state%p(i)
enddo
write(6,*)' '
call flush(6)
@ -808,9 +808,9 @@ subroutine homogenization_RGC_averageStressAndItsTangent(&
!$OMP CRITICAL (write2out)
do iGrain = 1,Ngrains
dPdF99 = math_Plain3333to99(dPdF(1:3,1:3,1:3,1:3,iGrain))
write(6,'(x,a30,x,i3)')'Stress tangent of grain: ',iGrain
write(6,'(1x,a30,1x,i3)')'Stress tangent of grain: ',iGrain
do i = 1,9
write(6,'(x,(e14.8,x))') (dPdF99(i,j), j = 1,9)
write(6,'(1x,(e14.8,1x))') (dPdF99(i,j), j = 1,9)
enddo
write(6,*)' '
enddo
@ -954,8 +954,8 @@ subroutine homogenization_RGC_stressPenalty(&
!* Debugging the surface correction factor
! if (ip == 1 .and. el == 1) then
! write(6,'(x,a20,2(x,i3))')'Correction factor: ',ip,el
! write(6,'(x,3(e10.4,x))')(surfCorr(i), i = 1,3)
! write(6,'(1x,a20,2(1x,i3))')'Correction factor: ',ip,el
! write(6,'(1x,3(e10.4,1x))')(surfCorr(i), i = 1,3)
! endif
!* -------------------------------------------------------------------------------------------------------------
@ -1003,11 +1003,11 @@ subroutine homogenization_RGC_stressPenalty(&
!* Debugging the mismatch tensor
! if (ip == 1 .and. el == 1) then
! write(6,'(x,a20,i2,x,a20,x,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
! write(6,'(x,3(e10.4,x))')(nDef(i,j), j = 1,3)
! write(6,'(1x,3(e10.4,1x))')(nDef(i,j), j = 1,3)
! enddo
! write(6,'(x,a20,e10.4))')'with magnitude: ',nDefNorm
! write(6,'(1x,a20,e10.4))')'with magnitude: ',nDefNorm
! endif
!* Compute the stress penalty of all interfaces
@ -1028,9 +1028,9 @@ subroutine homogenization_RGC_stressPenalty(&
!* Debugging the stress-like penalty
! if (ip == 1 .and. el == 1) then
! write(6,'(x,a20,i2)')'Penalty of grain: ',iGrain
! write(6,'(1x,a20,i2)')'Penalty of grain: ',iGrain
! do i = 1,3
! write(6,'(x,3(e10.4,x))')(rPen(i,j,iGrain), j = 1,3)
! write(6,'(1x,3(e10.4,1x))')(rPen(i,j,iGrain), j = 1,3)
! enddo
! endif
@ -1088,9 +1088,9 @@ subroutine homogenization_RGC_volumePenalty(&
!* Debugging the stress-like penalty of volume discrepancy
! if (ip == 1 .and. el == 1) then
! write(6,'(x,a30,i2)')'Volume penalty of grain: ',iGrain
! write(6,'(1x,a30,i2)')'Volume penalty of grain: ',iGrain
! do i = 1,3
! write(6,'(x,3(e10.4,x))')(vPen(i,j,iGrain), j = 1,3)
! write(6,'(1x,3(e10.4,1x))')(vPen(i,j,iGrain), j = 1,3)
! enddo
! endif
@ -1238,8 +1238,8 @@ function homogenization_RGC_interfaceNormal(&
! map the normal vector into sample coordinate system (basis)
! if (ip == 1 .and. el == 1) then
! write(6,'(x,a32,3(x,i3))')'Interface normal: ',intFace(1)
! write(6,'(x,3(e14.8,x))')(nVect(i), i = 1,3)
! 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,*)' '
! call flush(6)
! endif

4
code/lattice.f90
View File

@ -752,8 +752,8 @@ subroutine lattice_init()
if (debug_verbosity > 0) then
!$OMP CRITICAL (write2out)
write(6,'(a16,x,i5)') '# phases:',Nsections
write(6,'(a16,x,i5)') '# structures:',lattice_Nstructure
write(6,'(a16,1x,i5)') '# phases:',Nsections
write(6,'(a16,1x,i5)') '# structures:',lattice_Nstructure
write(6,*)
!$OMP END CRITICAL (write2out)
endif

23
code/makefile
View File

@ -31,6 +31,7 @@
# LAPACKROOT = pathinfo:LAPACK (will be adjusted by setup_code.py if present in pathinfo)
# PREFIX = arbitrary prefix
# SUFFIX = arbitrary suffix
# STANDARD_CHECK = checking for Fortran 2008, compiler dependend
########################################################################################
# Here are some useful debugging switches. 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/
@ -46,9 +47,6 @@ 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
#checks for standard
DEBUG5 =-stand std03/std95
#SUFFIX =$(DEBUG1) $(DEBUG2) $(DEBUG3)
########################################################################################
@ -115,20 +113,29 @@ endif
endif
endif
ifdef STANDARD_CHECK
STANDARD_CHECK_ifort = $(STANDARD_CHECK)
STANDARD_CHECK_gfortran = $(STANDARD_CHECK)
endif
STANDARD_CHECK_ifort ?= -stand f08
STANDARD_CHECK_gfortran ?=-std=f2008
OPTIMIZATION_OFF_ifort :=-O0
OPTIMIZATION_OFF_gfortran :=-O0
OPTIMIZATION_DEFENSIVE_ifort :=-O2
OPTIMIZATION_DEFENSIVE_gfortran :=-O2
OPTIMIZATION_AGGRESSIVE_ifort :=-O3 $(PORTABLE_SWITCH) -ip -static -fp-model fast=2 -no-prec-div
OPTIMIZATION_AGGRESSIVE_ifort :=-O3 $(PORTABLE_SWITCH) -ip -static -fp-model fast=2 -no-prec-div -xO
OPTIMIZATION_AGGRESSIVE_gfortran :=-O3 $(PORTABLE_SWITCH) -ffast-math -funroll-loops -ftree-vectorize
COMPILE_OPTIONS_ifort := -fpp -diag-disable 8291,8290
COMPILE_OPTIONS_gfortran := -xf95-cpp-input -ffixed-line-length-132 -fno-range-check
COMPILE_OPTIONS_ifort := -fpp -diag-disable 8291,8290,5268
#warning ID 9291,8290:
#warning ID 5268: Extension to standard: The text exceeds right hand column allowed on the line (we have only comments there)
COMPILE_OPTIONS_gfortran := -xf95-cpp-input -ffree-line-length-132 -fno-range-check
COMPILE = $(OPENMP_FLAG_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(OPTIMIZATION_$(OPTI)_$(F90)) -c
COMPILE_MAXOPTI = $(OPENMP_FLAG_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(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

14
code/material.f90
View File

@ -178,20 +178,20 @@ subroutine material_init()
write (6,*)
write (6,*) 'MATERIAL configuration'
write (6,*)
write (6,'(a32,x,a16,x,a6)') 'homogenization ','type ','grains'
write (6,'(a32,1x,a16,1x,a6)') 'homogenization ','type ','grains'
do i = 1,material_Nhomogenization
write (6,'(x,a32,x,a16,x,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
write (6,*)
write (6,'(a32,x,a11,x,a12,x,a13)') 'microstructure ','crystallite','constituents','homogeneous'
write (6,'(a32,1x,a11,1x,a12,1x,a13)') 'microstructure ','crystallite','constituents','homogeneous'
do i = 1,material_Nmicrostructure
write (6,'(a32,4x,i4,8x,i4,8x,l)') microstructure_name(i), &
write (6,'(a32,4x,i4,8x,i4,8x,l1)') microstructure_name(i), &
microstructure_crystallite(i), &
microstructure_Nconstituents(i), &
microstructure_elemhomo(i)
if (microstructure_Nconstituents(i) > 0_pInt) then
do j = 1,microstructure_Nconstituents(i)
write (6,'(a1,x,a32,x,a32,x,f6.4)') '>',phase_name(microstructure_phase(j,i)),&
write (6,'(a1,1x,a32,1x,a32,1x,f6.4)') '>',phase_name(microstructure_phase(j,i)),&
texture_name(microstructure_texture(j,i)),&
microstructure_fraction(j,i)
enddo
@ -667,7 +667,7 @@ subroutine material_populateGrains()
write (6,*)
write (6,*) 'MATERIAL grain population'
write (6,*)
write (6,'(a32,x,a32,x,a6)') 'homogenization_name','microstructure_name','grain#'
write (6,'(a32,1x,a32,1x,a6)') 'homogenization_name','microstructure_name','grain#'
!$OMP END CRITICAL (write2out)
endif
do homog = 1,material_Nhomogenization ! loop over homogenizations
@ -678,7 +678,7 @@ subroutine material_populateGrains()
if (debug_verbosity > 0) then
!$OMP CRITICAL (write2out)
write (6,*)
write (6,'(a32,x,a32,x,i6)') homogenization_name(homog),microstructure_name(micro),myNgrains
write (6,'(a32,1x,a32,1x,i6)') homogenization_name(homog),microstructure_name(micro),myNgrains
!$OMP END CRITICAL (write2out)
endif

28
code/math.f90
View File

@ -31,6 +31,7 @@
real(pReal), parameter :: pi = 3.14159265358979323846264338327950288419716939937510_pReal
real(pReal), parameter :: inDeg = 180.0_pReal/pi
real(pReal), parameter :: inRad = pi/180.0_pReal
complex(pReal), parameter :: two_pi_img = (0.0_pReal,2.0_pReal) * pi
! *** 3x3 Identity ***
real(pReal), dimension(3,3), parameter :: math_I3 = &
@ -153,7 +154,7 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = &
! comment the first random_seed call out, set randSize to 1, and use ifort
character(len=64) :: error_msg
!$OMP CRITICAL (write2out)
write(6,*)
write(6,*) ''
write(6,*) '<<<+- math init -+>>>'
write(6,*) '$Id$'
#include "compilation_info.f90"
@ -196,7 +197,7 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = &
q2 = math_AxisAngleToQuaternion(axisangle(1:3),axisangle(4))
if ( any(abs( q-q2) > tol_math_check) .and. &
any(abs(-q-q2) > tol_math_check) ) then
write (error_msg, '(a,e14.6)' ) 'maximum deviation',min(maxval(abs( q-q2)),maxval(abs(-q-q2)))
write (error_msg, '(a,e14.6)' ) 'maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2)))
call IO_error(670_pInt,ext_msg=error_msg)
endif
@ -205,7 +206,7 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = &
q2 = math_RToQuaternion(R)
if ( any(abs( q-q2) > tol_math_check) .and. &
any(abs(-q-q2) > tol_math_check) ) then
write (error_msg, '(a,e14.6)' ) 'maximum deviation',min(maxval(abs( q-q2)),maxval(abs(-q-q2)))
write (error_msg, '(a,e14.6)' ) 'maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2)))
call IO_error(671_pInt,ext_msg=error_msg)
endif
@ -214,7 +215,7 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = &
q2 = math_EulerToQuaternion(Eulers)
if ( any(abs( q-q2) > tol_math_check) .and. &
any(abs(-q-q2) > tol_math_check) ) then
write (error_msg, '(a,e14.6)' ) 'maximum deviation',min(maxval(abs( q-q2)),maxval(abs(-q-q2)))
write (error_msg, '(a,e14.6)' ) 'maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2)))
call IO_error(672_pInt,ext_msg=error_msg)
endif
@ -222,7 +223,7 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = &
Eulers = math_RToEuler(R);
R2 = math_EulerToR(Eulers)
if ( any(abs( R-R2) > tol_math_check) ) then
write (error_msg, '(a,e14.6)' ) 'maximum deviation',maxval(abs( R-R2))
write (error_msg, '(a,e14.6)' ) 'maximum deviation ',maxval(abs( R-R2))
call IO_error(673_pInt,ext_msg=error_msg)
endif
@ -3267,7 +3268,6 @@ subroutine deformed_fft(res,geomdim,defgrad_av,scaling,defgrad,coords)
! other variables
integer(pInt) :: i, j, k, res1_red
integer(pInt), dimension(3) :: k_s
complex(pReal), parameter :: integration_factor = cmplx(0.0_pReal,1.0_pReal)*pi*2.0_pReal
real(pReal), dimension(3) :: step, offset_coords
if (debug_verbosity > 0_pInt) then
@ -3317,11 +3317,11 @@ subroutine deformed_fft(res,geomdim,defgrad_av,scaling,defgrad,coords)
do i = 1_pInt, res1_red
k_s(1) = i-1_pInt
if(i/=1_pInt) coords_complex(i,j,k,1:3) = coords_complex(i,j,k,1:3)&
+ defgrad_complex(i,j,k,1:3,1)*geomdim(1)/(real(k_s(1),pReal)*integration_factor)
+ defgrad_complex(i,j,k,1:3,1)*geomdim(1)/(real(k_s(1),pReal)*two_pi_img)
if(j/=1_pInt) coords_complex(i,j,k,1:3) = coords_complex(i,j,k,1:3)&
+ defgrad_complex(i,j,k,1:3,2)*geomdim(2)/(real(k_s(2),pReal)*integration_factor)
+ defgrad_complex(i,j,k,1:3,2)*geomdim(2)/(real(k_s(2),pReal)*two_pi_img)
if(k/=1_pInt) coords_complex(i,j,k,1:3) = coords_complex(i,j,k,1:3)&
+ defgrad_complex(i,j,k,1:3,3)*geomdim(3)/(real(k_s(3),pReal)*integration_factor)
+ defgrad_complex(i,j,k,1:3,3)*geomdim(3)/(real(k_s(3),pReal)*two_pi_img)
enddo; enddo; enddo
call fftw_execute_dft_c2r(fftw_back,coords_complex,coords_real)
@ -3372,7 +3372,6 @@ subroutine curl_fft(res,geomdim,vec_tens,field,curl)
integer(pInt) i, j, k, l, res1_red
integer(pInt), dimension(3) :: k_s,cutting_freq
real(pReal) :: wgt
complex(pReal), parameter :: differentation_factor = cmplx(0.0_pReal,1.0_pReal)*2.0_pReal*pi ! cmplx(0.0_pReal, 2.0_pReal*pi) gets huge rounding error (casting to single prec?)
if (debug_verbosity > 0_pInt) then
print*, 'Calculating curl of vector/tensor field'
@ -3430,11 +3429,11 @@ subroutine curl_fft(res,geomdim,vec_tens,field,curl)
do k = 1, res(3); do j = 1, res(2);do i = 1, res1_red
do l = 1, vec_tens
curl_complex(i,j,k,l,1) = ( field_complex(i,j,k,l,3)*xi(i,j,k,2)&
-field_complex(i,j,k,l,2)*xi(i,j,k,3) )*differentation_factor
-field_complex(i,j,k,l,2)*xi(i,j,k,3) )*two_pi_img
curl_complex(i,j,k,l,2) = (-field_complex(i,j,k,l,3)*xi(i,j,k,1)&
+field_complex(i,j,k,l,1)*xi(i,j,k,3) )*differentation_factor
+field_complex(i,j,k,l,1)*xi(i,j,k,3) )*two_pi_img
curl_complex(i,j,k,l,3) = ( field_complex(i,j,k,l,2)*xi(i,j,k,1)&
-field_complex(i,j,k,l,1)*xi(i,j,k,2) )*differentation_factor
-field_complex(i,j,k,l,1)*xi(i,j,k,2) )*two_pi_img
enddo
enddo; enddo; enddo
@ -3477,7 +3476,6 @@ subroutine divergence_fft(res,geomdim,vec_tens,field,divergence)
integer(pInt) :: i, j, k, l, res1_red
real(pReal) :: wgt
integer(pInt), dimension(3) :: k_s,cutting_freq
complex(pReal), parameter :: differentation_factor = cmplx(0.0_pReal,1.0_pReal)*2.0_pReal*pi ! cmplx(0.0_pReal, 2.0_pReal*pi) gets huge rounding error (casting to single prec?)
if (debug_verbosity > 0_pInt) then
print '(a)', 'Calculating divergence of tensor/vector field using FFT'
@ -3533,7 +3531,7 @@ subroutine divergence_fft(res,geomdim,vec_tens,field,divergence)
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red
do l = 1_pInt, vec_tens
divergence_complex(i,j,k,l) = sum(field_complex(i,j,k,l,1:3)*xi(i,j,k,1:3))&
*differentation_factor
*two_pi_img
enddo
enddo; enddo; enddo
call fftw_execute_dft_c2r(fftw_back, divergence_complex, divergence_real)

68
code/mesh.f90
View File

@ -3360,32 +3360,32 @@ enddo
if (debug_verbosity > 0) then
!$OMP CRITICAL (write2out)
write (6,*)
write (6,*) "Input Parser: STATISTICS"
write (6,*) 'Input Parser: STATISTICS'
write (6,*)
write (6,*) mesh_Nelems, " : total number of elements in mesh"
write (6,*) mesh_NcpElems, " : total number of CP elements in mesh"
write (6,*) mesh_Nnodes, " : total number of nodes in mesh"
write (6,*) mesh_maxNnodes, " : max number of nodes in any CP element"
write (6,*) mesh_maxNips, " : max number of IPs in any CP element"
write (6,*) mesh_maxNipNeighbors, " : max number of IP neighbors in any CP element"
write (6,*) mesh_maxNsubNodes, " : max number of (additional) subnodes in any CP element"
write (6,*) mesh_maxNsharedElems, " : max number of CP elements sharing a node"
write (6,*) mesh_Nelems, ' : total number of elements in mesh'
write (6,*) mesh_NcpElems, ' : total number of CP elements in mesh'
write (6,*) mesh_Nnodes, ' : total number of nodes in mesh'
write (6,*) mesh_maxNnodes, ' : max number of nodes in any CP element'
write (6,*) mesh_maxNips, ' : max number of IPs in any CP element'
write (6,*) mesh_maxNipNeighbors, ' : max number of IP neighbors in any CP element'
write (6,*) mesh_maxNsubNodes, ' : max number of (additional) subnodes in any CP element'
write (6,*) mesh_maxNsharedElems, ' : max number of CP elements sharing a node'
write (6,*)
write (6,*) "Input Parser: HOMOGENIZATION/MICROSTRUCTURE"
write (6,*) 'Input Parser: HOMOGENIZATION/MICROSTRUCTURE'
write (6,*)
write (6,*) mesh_maxValStateVar(1), " : maximum homogenization index"
write (6,*) mesh_maxValStateVar(2), " : maximum microstructure index"
write (6,*) mesh_maxValStateVar(1), ' : maximum homogenization index'
write (6,*) mesh_maxValStateVar(2), ' : maximum microstructure index'
write (6,*)
write (fmt,"(a,i32.32,a)") "(9(x),a2,x,",mesh_maxValStateVar(2),"(i8))"
write (6,fmt) "+-",math_range(mesh_maxValStateVar(2))
write (fmt,"(a,i32.32,a)") "(i8,x,a2,x,",mesh_maxValStateVar(2),"(i8))"
write (fmt,'(a,i32.32,a)') '(9x,a2,1x,',mesh_maxValStateVar(2),'(i8))'
write (6,fmt) '+-',math_range(mesh_maxValStateVar(2))
write (fmt,'(a,i32.32,a)') '(i8,1x,a2,1x,',mesh_maxValStateVar(2),'(i8))'
do i=1,mesh_maxValStateVar(1) ! loop over all (possibly assigned) homogenizations
write (6,fmt) i,"| ",mesh_HomogMicro(i,:) ! loop over all (possibly assigned) microstructures
write (6,fmt) i,'| ',mesh_HomogMicro(i,:) ! loop over all (possibly assigned) microstructures
enddo
write(6,*)
write(6,*) "Input Parser: ADDITIONAL MPIE OPTIONS"
write(6,*) 'Input Parser: ADDITIONAL MPIE OPTIONS'
write(6,*)
write(6,*) "periodic surface : ", mesh_periodicSurface
write(6,*) 'periodic surface : ', mesh_periodicSurface
write(6,*)
call flush(6)
!$OMP END CRITICAL (write2out)
@ -3394,9 +3394,9 @@ endif
if (debug_verbosity > 1) then
!$OMP CRITICAL (write2out)
write (6,*)
write (6,*) "Input Parser: SUBNODE COORDINATES"
write (6,*) 'Input Parser: SUBNODE COORDINATES'
write (6,*)
write(6,'(a8,x,a5,x,a15,x,a15,x,a20,3(x,a12))') 'elem','IP','IP neighbor','IPFaceNodes','subNodeOnIPFace','x','y','z'
write(6,'(a8,1x,a5,1x,a15,1x,a15,1x,a20,3(1x,a12))') 'elem','IP','IP neighbor','IPFaceNodes','subNodeOnIPFace','x','y','z'
do e = 1,mesh_NcpElems ! loop over cpElems
if (debug_selectiveDebugger .and. debug_e /= e) cycle
t = mesh_element(2,e) ! get elemType
@ -3404,7 +3404,7 @@ if (debug_verbosity > 1) then
if (debug_selectiveDebugger .and. debug_i /= i) cycle
do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP
do n = 1,FE_NipFaceNodes ! loop over nodes on interface
write(6,'(i8,x,i5,x,i15,x,i15,x,i20,3(x,f12.8))') e,i,f,n,FE_subNodeOnIPFace(n,f,i,t),&
write(6,'(i8,1x,i5,1x,i15,1x,i15,1x,i20,3(1x,f12.8))') e,i,f,n,FE_subNodeOnIPFace(n,f,i,t),&
mesh_subNodeCoord(1,FE_subNodeOnIPFace(n,f,i,t),e),&
mesh_subNodeCoord(2,FE_subNodeOnIPFace(n,f,i,t),e),&
mesh_subNodeCoord(3,FE_subNodeOnIPFace(n,f,i,t),e)
@ -3414,52 +3414,52 @@ if (debug_verbosity > 1) then
enddo
write(6,*)
write(6,*) 'Input Parser: IP COORDINATES'
write(6,'(a8,x,a5,3(x,a12))') 'elem','IP','x','y','z'
write(6,'(a8,1x,a5,3(1x,a12))') 'elem','IP','x','y','z'
do e = 1,mesh_NcpElems
if (debug_selectiveDebugger .and. debug_e /= e) cycle
do i = 1,FE_Nips(mesh_element(2,e))
if (debug_selectiveDebugger .and. debug_i /= i) cycle
write (6,'(i8,x,i5,3(x,f12.8))') e, i, mesh_ipCenterOfGravity(:,i,e)
write (6,'(i8,1x,i5,3(1x,f12.8))') e, i, mesh_ipCenterOfGravity(:,i,e)
enddo
enddo
write (6,*)
write (6,*) "Input Parser: ELEMENT VOLUME"
write (6,*) 'Input Parser: ELEMENT VOLUME'
write (6,*)
write (6,"(a13,x,e15.8)") "total volume", sum(mesh_ipVolume)
write (6,'(a13,1x,e15.8)') 'total volume', sum(mesh_ipVolume)
write (6,*)
write (6,"(a8,x,a5,x,a15,x,a5,x,a15,x,a16)") "elem","IP","volume","face","area","-- normal --"
write (6,'(a8,1x,a5,1x,a15,1x,a5,1x,a15,1x,a16)') 'elem','IP','volume','face','area','-- normal --'
do e = 1,mesh_NcpElems
if (debug_selectiveDebugger .and. debug_e /= e) cycle
do i = 1,FE_Nips(mesh_element(2,e))
if (debug_selectiveDebugger .and. debug_i /= i) cycle
write (6,"(i8,x,i5,x,e15.8)") e,i,mesh_IPvolume(i,e)
write (6,'(i8,1x,i5,1x,e15.8)') e,i,mesh_IPvolume(i,e)
do f = 1,FE_NipNeighbors(mesh_element(2,e))
write (6,"(i33,x,e15.8,x,3(f6.3,x))") f,mesh_ipArea(f,i,e),mesh_ipAreaNormal(:,f,i,e)
write (6,'(i33,1x,e15.8,1x,3(f6.3,1x))') f,mesh_ipArea(f,i,e),mesh_ipAreaNormal(:,f,i,e)
enddo
enddo
enddo
write (6,*)
write (6,*) "Input Parser: NODE TWINS"
write (6,*) 'Input Parser: NODE TWINS'
write (6,*)
write(6,'(a6,3(3(x),a6))') ' node','twin_x','twin_y','twin_z'
write(6,'(a6,3(3x,a6))') ' node','twin_x','twin_y','twin_z'
do n = 1,mesh_Nnodes ! loop over cpNodes
if (debug_e <= mesh_NcpElems) then
if (any(mesh_element(5:,debug_e) == n)) then
write(6,'(i6,3(3(x),i6))') n, mesh_nodeTwins(1:3,n)
write(6,'(i6,3(3x,i6))') n, mesh_nodeTwins(1:3,n)
endif
endif
enddo
write(6,*)
write(6,*) "Input Parser: IP NEIGHBORHOOD"
write(6,*) 'Input Parser: IP NEIGHBORHOOD'
write(6,*)
write(6,"(a8,x,a10,x,a10,x,a3,x,a13,x,a13)") "elem","IP","neighbor","","elemNeighbor","ipNeighbor"
write(6,'(a8,1x,a10,1x,a10,1x,a3,1x,a13,1x,a13)') 'elem','IP','neighbor','','elemNeighbor','ipNeighbor'
do e = 1,mesh_NcpElems ! loop over cpElems
if (debug_selectiveDebugger .and. debug_e /= e) cycle
t = mesh_element(2,e) ! get elemType
do i = 1,FE_Nips(t) ! loop over IPs of elem
if (debug_selectiveDebugger .and. debug_i /= i) cycle
do n = 1,FE_NipNeighbors(t) ! loop over neighbors of IP
write (6,"(i8,x,i10,x,i10,x,a3,x,i13,x,i13)") e,i,n,'-->',mesh_ipNeighborhood(1,n,i,e),mesh_ipNeighborhood(2,n,i,e)
write (6,'(i8,1x,i10,1x,i10,1x,a3,1x,i13,1x,i13)') e,i,n,'-->',mesh_ipNeighborhood(1,n,i,e),mesh_ipNeighborhood(2,n,i,e)
enddo
enddo
enddo

98
code/numerics.f90
View File

@ -278,67 +278,67 @@ subroutine numerics_init()
! writing parameters to output file
!$OMP CRITICAL (write2out)
write(6,'(a24,x,e8.1)') ' relevantStrain: ',relevantStrain
write(6,'(a24,x,e8.1)') ' defgradTolerance: ',defgradTolerance
write(6,'(a24,x,i8)') ' iJacoStiffness: ',iJacoStiffness
write(6,'(a24,x,i8)') ' iJacoLpresiduum: ',iJacoLpresiduum
write(6,'(a24,x,e8.1)') ' pert_Fg: ',pert_Fg
write(6,'(a24,x,i8)') ' pert_method: ',pert_method
write(6,'(a24,x,i8)') ' nCryst: ',nCryst
write(6,'(a24,x,e8.1)') ' subStepMinCryst: ',subStepMinCryst
write(6,'(a24,x,e8.1)') ' subStepSizeCryst: ',subStepSizeCryst
write(6,'(a24,x,e8.1)') ' stepIncreaseCryst: ',stepIncreaseCryst
write(6,'(a24,x,i8)') ' nState: ',nState
write(6,'(a24,x,i8)') ' nStress: ',nStress
write(6,'(a24,x,e8.1)') ' rTol_crystalliteState: ',rTol_crystalliteState
write(6,'(a24,x,e8.1)') ' rTol_crystalliteTemp: ',rTol_crystalliteTemperature
write(6,'(a24,x,e8.1)') ' rTol_crystalliteStress: ',rTol_crystalliteStress
write(6,'(a24,x,e8.1)') ' aTol_crystalliteStress: ',aTol_crystalliteStress
write(6,'(a24,2(x,i8),/)')' integrator: ',numerics_integrator
write(6,'(a24,1x,e8.1)') ' relevantStrain: ',relevantStrain
write(6,'(a24,1x,e8.1)') ' defgradTolerance: ',defgradTolerance
write(6,'(a24,1x,i8)') ' iJacoStiffness: ',iJacoStiffness
write(6,'(a24,1x,i8)') ' iJacoLpresiduum: ',iJacoLpresiduum
write(6,'(a24,1x,e8.1)') ' pert_Fg: ',pert_Fg
write(6,'(a24,1x,i8)') ' pert_method: ',pert_method
write(6,'(a24,1x,i8)') ' nCryst: ',nCryst
write(6,'(a24,1x,e8.1)') ' subStepMinCryst: ',subStepMinCryst
write(6,'(a24,1x,e8.1)') ' subStepSizeCryst: ',subStepSizeCryst
write(6,'(a24,1x,e8.1)') ' stepIncreaseCryst: ',stepIncreaseCryst
write(6,'(a24,1x,i8)') ' nState: ',nState
write(6,'(a24,1x,i8)') ' nStress: ',nStress
write(6,'(a24,1x,e8.1)') ' rTol_crystalliteState: ',rTol_crystalliteState
write(6,'(a24,1x,e8.1)') ' rTol_crystalliteTemp: ',rTol_crystalliteTemperature
write(6,'(a24,1x,e8.1)') ' rTol_crystalliteStress: ',rTol_crystalliteStress
write(6,'(a24,1x,e8.1)') ' aTol_crystalliteStress: ',aTol_crystalliteStress
write(6,'(a24,2(1x,i8),/)')' integrator: ',numerics_integrator
write(6,'(a24,x,i8)') ' nHomog: ',nHomog
write(6,'(a24,x,e8.1)') ' subStepMinHomog: ',subStepMinHomog
write(6,'(a24,x,e8.1)') ' subStepSizeHomog: ',subStepSizeHomog
write(6,'(a24,x,e8.1)') ' stepIncreaseHomog: ',stepIncreaseHomog
write(6,'(a24,x,i8,/)') ' nMPstate: ',nMPstate
write(6,'(a24,1x,i8)') ' nHomog: ',nHomog
write(6,'(a24,1x,e8.1)') ' subStepMinHomog: ',subStepMinHomog
write(6,'(a24,1x,e8.1)') ' subStepSizeHomog: ',subStepSizeHomog
write(6,'(a24,1x,e8.1)') ' stepIncreaseHomog: ',stepIncreaseHomog
write(6,'(a24,1x,i8,/)') ' nMPstate: ',nMPstate
!* RGC parameters
write(6,'(a24,x,e8.1)') ' aTol_RGC: ',absTol_RGC
write(6,'(a24,x,e8.1)') ' rTol_RGC: ',relTol_RGC
write(6,'(a24,x,e8.1)') ' aMax_RGC: ',absMax_RGC
write(6,'(a24,x,e8.1)') ' rMax_RGC: ',relMax_RGC
write(6,'(a24,x,e8.1)') ' perturbPenalty_RGC: ',pPert_RGC
write(6,'(a24,x,e8.1)') ' relevantMismatch_RGC: ',xSmoo_RGC
write(6,'(a24,x,e8.1)') ' viscosityrate_RGC: ',viscPower_RGC
write(6,'(a24,x,e8.1)') ' viscositymodulus_RGC: ',viscModus_RGC
write(6,'(a24,x,e8.1)') ' maxrelaxation_RGC: ',maxdRelax_RGC
write(6,'(a24,x,e8.1)') ' maxVolDiscrepancy_RGC: ',maxVolDiscr_RGC
write(6,'(a24,x,e8.1)') ' volDiscrepancyMod_RGC: ',volDiscrMod_RGC
write(6,'(a24,x,e8.1,/)') ' discrepancyPower_RGC: ',volDiscrPow_RGC
write(6,'(a24,1x,e8.1)') ' aTol_RGC: ',absTol_RGC
write(6,'(a24,1x,e8.1)') ' rTol_RGC: ',relTol_RGC
write(6,'(a24,1x,e8.1)') ' aMax_RGC: ',absMax_RGC
write(6,'(a24,1x,e8.1)') ' rMax_RGC: ',relMax_RGC
write(6,'(a24,1x,e8.1)') ' perturbPenalty_RGC: ',pPert_RGC
write(6,'(a24,1x,e8.1)') ' relevantMismatch_RGC: ',xSmoo_RGC
write(6,'(a24,1x,e8.1)') ' viscosityrate_RGC: ',viscPower_RGC
write(6,'(a24,1x,e8.1)') ' viscositymodulus_RGC: ',viscModus_RGC
write(6,'(a24,1x,e8.1)') ' maxrelaxation_RGC: ',maxdRelax_RGC
write(6,'(a24,1x,e8.1)') ' maxVolDiscrepancy_RGC: ',maxVolDiscr_RGC
write(6,'(a24,1x,e8.1)') ' volDiscrepancyMod_RGC: ',volDiscrMod_RGC
write(6,'(a24,1x,e8.1,/)') ' discrepancyPower_RGC: ',volDiscrPow_RGC
!* spectral parameters
write(6,'(a24,x,e8.1)') ' err_div_tol: ',err_div_tol
write(6,'(a24,x,e8.1)') ' err_stress_tolrel: ',err_stress_tolrel
write(6,'(a24,x,i8)') ' itmax: ',itmax
write(6,'(a24,x,L8)') ' memory_efficient: ',memory_efficient
write(6,'(a24,1x,e8.1)') ' err_div_tol: ',err_div_tol
write(6,'(a24,1x,e8.1)') ' err_stress_tolrel: ',err_stress_tolrel
write(6,'(a24,1x,i8)') ' itmax: ',itmax
write(6,'(a24,1x,L8)') ' memory_efficient: ',memory_efficient
if(fftw_timelimit<0.0_pReal) then
write(6,'(a24,x,L8)') ' fftw_timelimit: ',.false.
write(6,'(a24,1x,L8)') ' fftw_timelimit: ',.false.
else
write(6,'(a24,x,e8.1)') ' fftw_timelimit: ',fftw_timelimit
write(6,'(a24,1x,e8.1)') ' fftw_timelimit: ',fftw_timelimit
endif
write(6,'(a24,x,a)') ' fftw_planner_string: ',trim(fftw_planner_string)
write(6,'(a24,x,i8)') ' fftw_planner_flag: ',fftw_planner_flag
write(6,'(a24,x,e8.1)') ' rotation_tol: ',rotation_tol
write(6,'(a24,x,L8,/)') ' divergence_correction: ',divergence_correction
write(6,'(a24,x,L8,/)') ' update_gamma: ',update_gamma
write(6,'(a24,x,L8,/)') ' simplified_algorithm: ',simplified_algorithm
write(6,'(a24,x,e8.1)') ' cut_off_value: ',cut_off_value
write(6,'(a24,1x,a)') ' fftw_planner_string: ',trim(fftw_planner_string)
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,L8,/)') ' divergence_correction: ',divergence_correction
write(6,'(a24,1x,L8,/)') ' update_gamma: ',update_gamma
write(6,'(a24,1x,L8,/)') ' simplified_algorithm: ',simplified_algorithm
write(6,'(a24,1x,e8.1)') ' cut_off_value: ',cut_off_value
!* Random seeding parameters
write(6,'(a24,x,i16,/)') ' fixed_seed: ',fixedSeed
write(6,'(a24,1x,i16,/)') ' fixed_seed: ',fixedSeed
!$OMP END CRITICAL (write2out)
!* openMP parameter
!$ write(6,'(a24,x,i8,/)') ' number of threads: ',DAMASK_NumThreadsInt
!$ write(6,'(a24,1x,i8,/)') ' number of threads: ',DAMASK_NumThreadsInt
! sanity check
if (relevantStrain <= 0.0_pReal) call IO_error(260)

12
code/prec.f90
View File

@ -34,7 +34,7 @@ real(pReal), parameter :: tol_gravityNodePos = 1.0e-100_pReal
! NaN is precision dependent
! from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html
! copy can be found in documentation/Code/Fortran
real(pReal), parameter :: DAMASK_NaN = Z'7FF0000000000001'
real(pReal), parameter :: DAMASK_NaN = real(Z'7FF0000000000001', pReal)
type :: p_vec
real(pReal), dimension(:), pointer :: p
end type p_vec
@ -49,11 +49,11 @@ implicit none
write(6,*) '<<<+- prec init -+>>>'
write(6,*) '$Id$'
#include "compilation_info.f90"
write(6,'(a,i3)'), ' Bytes for pReal: ',pReal
write(6,'(a,i3)'), ' Bytes for pInt: ',pInt
write(6,'(a,i3)'), ' Bytes for pLongInt: ',pLongInt
write(6,'(a,e3.3)'), ' NaN: ',DAMASK_NAN
write(6,'(a,l3)'), ' NaN /= NaN: ',DAMASK_NaN/=DAMASK_NaN
write(6,'(a,i3)') ' Bytes for pReal: ',pReal
write(6,'(a,i3)') ' Bytes for pInt: ',pInt
write(6,'(a,i3)') ' Bytes for pLongInt: ',pLongInt
write(6,'(a,e3.3)') ' NaN: ',DAMASK_NAN
write(6,'(a,l3)') ' NaN /= NaN: ',DAMASK_NaN/=DAMASK_NaN
write(6,*)
!$OMP END CRITICAL (write2out)

12
code/prec_single.f90
View File

@ -34,7 +34,7 @@ real(pReal), parameter :: tol_gravityNodePos = 1.0e-36_pReal
! NaN is precision dependent
! from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html
! copy can be found in documentation/Code/Fortran
real(pReal), parameter :: DAMASK_NaN = Z'7F800001'
real(pReal), parameter :: DAMASK_NaN = real(Z'7F800001', pReal)
type :: p_vec
real(pReal), dimension(:), pointer :: p
end type p_vec
@ -49,11 +49,11 @@ implicit none
write(6,*) '<<<+- prec_single init -+>>>'
write(6,*) '$Id$'
#include "compilation_info.f90"
write(6,'(a,i3)'), ' Bytes for pReal: ',pReal
write(6,'(a,i3)'), ' Bytes for pInt: ',pInt
write(6,'(a,i3)'), ' Bytes for pLongInt: ',pLongInt
write(6,'(a,e3.3)'), ' NaN: ',DAMASK_NAN
write(6,'(a,l3)'), ' NaN /= NaN: ',DAMASK_NaN/=DAMASK_NaN
write(6,'(a,i3)') ' Bytes for pReal: ',pReal
write(6,'(a,i3)') ' Bytes for pInt: ',pInt
write(6,'(a,i3)') ' Bytes for pLongInt: ',pLongInt
write(6,'(a,e3.3)') ' NaN: ',DAMASK_NAN
write(6,'(a,l3)') ' NaN /= NaN: ',DAMASK_NaN/=DAMASK_NaN
write(6,*)
!$OMP END CRITICAL (write2out)