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

debug.config, debug.f90, DAMASK_abaqus_exp.f, DAMASK_abaqus_std.f: changed to new debug scheme (wasn't working) 

lattice.f90, FEsolving.f90: explicitly defined public functions and variables, all others are now private
numerics.f90: changed output format of real numbers, now instead of 0.1eX 1.0e(X-1) is printed to screen
Makefile: now using correct Optimization flags for OPTIMIZATION=AGGRESSIVE
DAMASK_spectral_AL.f90: improved, but still testing. Stress BCs now seem to be handled correctly
This commit is contained in:
Martin Diehl 2012-03-20 18:01:31 +00:00
parent 04d83ec9cd
commit dbdc7ddfa2
9 changed files with 771 additions and 713 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
code/DAMASK_abaqus_exp.f
View File

@ -177,7 +177,9 @@ subroutine vumat (jblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, &
use math, only: invnrmMandel use math, only: invnrmMandel
use debug, only: debug_info, & use debug, only: debug_info, &
debug_reset, & debug_reset, &
debug_verbosity debug_levelBasic, &
debug_what, &
debug_abaqus
use mesh, only: mesh_FEasCP use mesh, only: mesh_FEasCP
use CPFEM, only: CPFEM_general,CPFEM_init_done, CPFEM_initAll use CPFEM, only: CPFEM_general,CPFEM_init_done, CPFEM_initAll
use homogenization, only: materialpoint_sizeResults, materialpoint_results use homogenization, only: materialpoint_sizeResults, materialpoint_results
@ -216,7 +218,7 @@ subroutine vumat (jblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, &
call CPFEM_initAll(temp,nElement(n),nMatPoint(n)) call CPFEM_initAll(temp,nElement(n),nMatPoint(n))
outdatedByNewInc = .false. outdatedByNewInc = .false.
if ( debug_verbosity > 1 ) then if (iand(debug_what(debug_abaqus),debug_levelBasic) /= 0) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(i8,x,i2,x,a)') nElement(n),nMatPoint(n),'first call special case..!'; call flush(6) write(6,'(i8,x,i2,x,a)') nElement(n),nMatPoint(n),'first call special case..!'; call flush(6)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
@ -225,7 +227,7 @@ subroutine vumat (jblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, &
else if (theTime < totalTime) then ! reached convergence else if (theTime < totalTime) then ! reached convergence
outdatedByNewInc = .true. outdatedByNewInc = .true.
if ( debug_verbosity > 1 ) then if (iand(debug_what(debug_abaqus),debug_levelBasic) /= 0) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write (6,'(i8,x,i2,x,a)') nElement(n),nMatPoint(n),'lastIncConverged + outdated'; call flush(6) write (6,'(i8,x,i2,x,a)') nElement(n),nMatPoint(n),'lastIncConverged + outdated'; call flush(6)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
@ -247,7 +249,7 @@ subroutine vumat (jblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, &
theTime = totalTime ! record current starting time theTime = totalTime ! record current starting time
if ( debug_verbosity > 1 ) then if (iand(debug_what(debug_abaqus),debug_levelBasic) /= 0) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(a16,x,i2,x,a,i8,x,i5,a)') 'computationMode',computationMode,'(',nElement(n),nMatPoint(n),')'; call flush(6) write(6,'(a16,x,i2,x,a,i8,x,i5,a)') 'computationMode',computationMode,'(',nElement(n),nMatPoint(n),')'; call flush(6)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)

8
code/DAMASK_abaqus_std.f
View File

@ -141,7 +141,9 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,&
use math, only: invnrmMandel use math, only: invnrmMandel
use debug, only: debug_info, & use debug, only: debug_info, &
debug_reset, & debug_reset, &
debug_verbosity debug_levelBasic, &
debug_what, &
debug_abaqus
use mesh, only: mesh_FEasCP use mesh, only: mesh_FEasCP
use CPFEM, only: CPFEM_general,CPFEM_init_done, CPFEM_initAll use CPFEM, only: CPFEM_general,CPFEM_init_done, CPFEM_initAll
use homogenization, only: materialpoint_sizeResults, materialpoint_results use homogenization, only: materialpoint_sizeResults, materialpoint_results
@ -167,7 +169,7 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,&
real(pReal), dimension(6,6) :: ddsdde_h real(pReal), dimension(6,6) :: ddsdde_h
integer(pInt) computationMode, i, cp_en integer(pInt) computationMode, i, cp_en
if (debug_verbosity > 1 .and. noel == 1 .and. npt == 1) then if (iand(debug_what(debug_abaqus),debug_levelBasic) /= 0 .and. noel == 1 .and. npt == 1) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) 'el',noel,'ip',npt write(6,*) 'el',noel,'ip',npt
write(6,*) 'got kinc as',kinc write(6,*) 'got kinc as',kinc
@ -255,7 +257,7 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,&
theInc = kinc ! record current increment number theInc = kinc ! record current increment number
lastMode = calcMode(npt,cp_en) ! record calculationMode lastMode = calcMode(npt,cp_en) ! record calculationMode
if ( debug_verbosity > 1 ) then if (iand(debug_what(debug_abaqus),debug_levelBasic) /= 0) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(a16,x,i2,x,a,i8,a,i8,x,i5,a)') 'computationMode',computationMode,'(',cp_en,':',noel,npt,')'; call flush(6) write(6,'(a16,x,i2,x,a,i8,a,i8,x,i5,a)') 'computationMode',computationMode,'(',cp_en,':',noel,npt,')'; call flush(6)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)

49
code/DAMASK_spectral_AL.f90
View File

@ -36,7 +36,7 @@
!################################################################################################## !##################################################################################################
! used modules ! used modules
!################################################################################################## !##################################################################################################
program DAMASK_spectral program DAMASK_spectral_AL
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment)
use DAMASK_interface use DAMASK_interface
@ -83,7 +83,7 @@ program DAMASK_spectral
! variable storing information from load case file ! variable storing information from load case file
type bc_type type bc_type
real(pReal), dimension (3,3) :: deformation = 0.0_pReal, & ! applied velocity gradient or time derivative of deformation gradient real(pReal), dimension (3,3) :: deformation = 0.0_pReal, & ! applied velocity gradient or time derivative of deformation gradient
stress = 0.0_pReal, & ! stress BC (if applicable) P = 0.0_pReal, & ! stress BC (if applicable)
rotation = math_I3 ! rotation of BC (if applicable) rotation = math_I3 ! rotation of BC (if applicable)
real(pReal) :: time = 0.0_pReal, & ! length of increment real(pReal) :: time = 0.0_pReal, & ! length of increment
temperature = 300.0_pReal ! isothermal starting conditions temperature = 300.0_pReal ! isothermal starting conditions
@ -230,10 +230,10 @@ program DAMASK_spectral
IO_stringValue(line,positions,j+k) /= '*' IO_stringValue(line,positions,j+k) /= '*'
do k = 1_pInt,9_pInt do k = 1_pInt,9_pInt
if (bc(loadcase)%maskStressVector(k)) temp_valueVector(k) =& if (bc(loadcase)%maskStressVector(k)) temp_valueVector(k) =&
IO_floatValue(line,positions,j+k) ! assign values for the bc(loadcase)%stress matrix IO_floatValue(line,positions,j+k) ! assign values for the bc(loadcase)%P matrix
enddo enddo
bc(loadcase)%maskStress = transpose(reshape(bc(loadcase)%maskStressVector,[ 3,3])) bc(loadcase)%maskStress = transpose(reshape(bc(loadcase)%maskStressVector,[ 3,3]))
bc(loadcase)%stress = math_plain9to33(temp_valueVector) bc(loadcase)%P = math_plain9to33(temp_valueVector)
case('t','time','delta') ! increment time case('t','time','delta') ! increment time
bc(loadcase)%time = IO_floatValue(line,positions,j+1_pInt) bc(loadcase)%time = IO_floatValue(line,positions,j+1_pInt)
case('temp','temperature') ! starting temperature case('temp','temperature') ! starting temperature
@ -379,7 +379,7 @@ program DAMASK_spectral
write (*,'(3(3(f12.7,1x)/))',advance='no') merge(math_transpose33(bc(loadcase)%deformation),& write (*,'(3(3(f12.7,1x)/))',advance='no') merge(math_transpose33(bc(loadcase)%deformation),&
reshape(spread(DAMASK_NaN,1,9),[ 3,3]),transpose(bc(loadcase)%maskDeformation)) reshape(spread(DAMASK_NaN,1,9),[ 3,3]),transpose(bc(loadcase)%maskDeformation))
write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') 'stress / GPa:',& write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') 'stress / GPa:',&
1e-9_pReal*merge(math_transpose33(bc(loadcase)%stress),& 1e-9_pReal*merge(math_transpose33(bc(loadcase)%P),&
reshape(spread(DAMASK_NaN,1,9),[ 3,3]),transpose(bc(loadcase)%maskStress)) reshape(spread(DAMASK_NaN,1,9),[ 3,3]),transpose(bc(loadcase)%maskStress))
if (any(bc(loadcase)%rotation /= math_I3)) & if (any(bc(loadcase)%rotation /= math_I3)) &
write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') ' rotation of loadframe:',& write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') ' rotation of loadframe:',&
@ -672,7 +672,7 @@ program DAMASK_spectral
guessmode = 1.0_pReal ! keep guessing along former trajectory during same loadcase guessmode = 1.0_pReal ! keep guessing along former trajectory during same loadcase
CPFEM_mode = 1_pInt ! winding forward CPFEM_mode = 1_pInt ! winding forward
iter = 0_pInt iter = 0_pInt
err_crit = 2.0_pReal * err_div_tol ! go into loop err_crit = huge(err_div_tol) ! go into loop
!################################################################################################## !##################################################################################################
! convergence loop (looping over iterations) ! convergence loop (looping over iterations)
@ -691,21 +691,19 @@ program DAMASK_spectral
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress BC handling ! stress BC handling
if(size_reduced > 0_pInt) then ! calculate stress BC if applied if(size_reduced > 0_pInt) then ! calculate stress BC if applied
err_stress = maxval(abs(mask_stress * (P_av - bc(loadcase)%stress))) ! maximum deviaton (tensor norm not applicable) err_stress = maxval(abs(mask_stress * (P_av - bc(loadcase)%P))) ! maximum deviaton (tensor norm not applicable)
write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'stress deviation =',& write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'stress deviation =',&
math_transpose33(mask_stress * (P_av - bc(loadcase)%stress))/1.0e6_pReal math_transpose33(mask_stress * (P_av - bc(loadcase)%P))/1.0e6_pReal
F_aim = F_aim + math_mul3333xx33(S_lastInc,bc(loadcase)%stress- P_av) F_aim = F_aim + math_mul3333xx33(S_lastInc,bc(loadcase)%P- P_av)
err_stress_tol = maxval(abs(P_av)) * err_stress_tolrel ! don't use any tensor norm because the comparison should be coherent err_stress_tol = maxval(abs(P_av)) * err_stress_tolrel ! don't use any tensor norm because the comparison should be coherent
else else
err_stress_tol = + huge(1.0_pReal) err_stress_tol = + huge(1.0_pReal)
endif endif
F_aim_lab = math_rotate_backward33(F_aim,bc(loadcase)%rotation) F_aim_lab = math_rotate_backward33(F_aim,bc(loadcase)%rotation)
write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'F =',& write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'F aim =',&
math_transpose33(F_aim) math_transpose33(F_aim)
temp33_real = 0.0_pReal write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'F* =',&
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) math_transpose33(F_star_av)
temp33_real = temp33_real + F_star(i,j,k,1:3,1:3)
enddo; enddo; enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! doing Fourier transform ! doing Fourier transform
@ -750,6 +748,11 @@ program DAMASK_spectral
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! doing inverse Fourier transform ! doing inverse Fourier transform
call fftw_execute_dft_c2r(plan_correction,F_fourier,F_real) ! back transform of fluct deformation gradient call fftw_execute_dft_c2r(plan_correction,F_fourier,F_real) ! back transform of fluct deformation gradient
! do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
! write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'delta F real =',&
! math_transpose33(F_real(i,j,k,1:3,1:3)*wgt)
! enddo; enddo; enddo
F_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = F_real(1:res(1),1:res(2),1:res(3),1:3,1:3) * wgt + & F_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = F_real(1:res(1),1:res(2),1:res(3),1:3,1:3) * wgt + &
F_star(1:res(1),1:res(2),1:res(3),1:3,1:3) F_star(1:res(1),1:res(2),1:res(3),1:3,1:3)
@ -757,18 +760,21 @@ program DAMASK_spectral
! !
print '(a)', '... update stress field P(F*) .....................................' print '(a)', '... update stress field P(F*) .....................................'
ielem = 0_pInt ielem = 0_pInt
temp33_Real = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
ielem = ielem + 1_pInt ielem = ielem + 1_pInt
call CPFEM_general(3_pInt,& ! collect cycle call CPFEM_general(3_pInt,& ! collect cycle
coordinates(i,j,k,1:3), F_lastInc(i,j,k,1:3,1:3),& coordinates(i,j,k,1:3), F_lastInc(i,j,k,1:3,1:3),&
F_star(i,j,k,1:3,1:3),temperature(i,j,k),timeinc,ielem,1_pInt,& F_star(i,j,k,1:3,1:3),temperature(i,j,k),timeinc,ielem,1_pInt,&
sigma,dsde, P, dPdF) sigma,dsde, P, dPdF)
temp33_Real = temp33_Real + F_real(i,j,k,1:3,1:3)
enddo; enddo; enddo enddo; enddo; enddo
write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'F =',&
math_transpose33(temp33_Real*wgt)
ielem = 0_pInt ielem = 0_pInt
err_f = 0.0_pReal err_f = 0.0_pReal
F_star_av = 0.0_pReal F_star_av = 0.0_pReal
P_av =0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
ielem = ielem + 1_pInt ielem = ielem + 1_pInt
call CPFEM_general(CPFEM_mode,& call CPFEM_general(CPFEM_mode,&
@ -776,9 +782,12 @@ program DAMASK_spectral
F_star(i,j,k,1:3,1:3),temperature(i,j,k),timeinc,ielem,1_pInt,& F_star(i,j,k,1:3,1:3),temperature(i,j,k),timeinc,ielem,1_pInt,&
sigma,dsde, P,dPdF) sigma,dsde, P,dPdF)
CPFEM_mode = 2_pInt ! winding forward CPFEM_mode = 2_pInt ! winding forward
P_av = P_av + P
if (iter == 1_pInt) lambda(i,j,k,1:3,1:3) = P
temp33_Real = lambda(i,j,k,1:3,1:3) - P & temp33_Real = lambda(i,j,k,1:3,1:3) - P &
+ math_mul3333xx33(C_inc0,F_real(i,j,k,1:3,1:3)- F_star(i,j,k,1:3,1:3)) + math_mul3333xx33(C_inc0,F_real(i,j,k,1:3,1:3)- F_star(i,j,k,1:3,1:3))
! write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'F - F* =',&
! math_transpose33(F_real(i,j,k,1:3,1:3)- F_star(i,j,k,1:3,1:3))
F_star(i,j,k,1:3,1:3) = F_star(i,j,k,1:3,1:3) + & F_star(i,j,k,1:3,1:3) = F_star(i,j,k,1:3,1:3) + &
math_mul3333xx33(math_invSym3333(C_inc0 + dPdF), temp33_Real) math_mul3333xx33(math_invSym3333(C_inc0 + dPdF), temp33_Real)
lambda(i,j,k,1:3,1:3) = lambda(i,j,k,1:3,1:3) + math_mul3333xx33(C_inc0,F_real(i,j,k,1:3,1:3) & lambda(i,j,k,1:3,1:3) = lambda(i,j,k,1:3,1:3) + math_mul3333xx33(C_inc0,F_real(i,j,k,1:3,1:3) &
@ -787,9 +796,7 @@ program DAMASK_spectral
temp33_real = F_star(i,j,k,1:3,1:3) - F_real(i,j,k,1:3,1:3) temp33_real = F_star(i,j,k,1:3,1:3) - F_real(i,j,k,1:3,1:3)
err_f = max(err_f, sqrt(math_mul33xx33(temp33_real,temp33_real))) err_f = max(err_f, sqrt(math_mul33xx33(temp33_real,temp33_real)))
enddo; enddo; enddo enddo; enddo; enddo
P_av = math_rotate_forward33(P_av * wgt,bc(loadcase)%rotation)
write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'P(F*) =',&
math_transpose33(P_av)/1.e6_pReal
F_star_av = F_star_av *wgt F_star_av = F_star_av *wgt
write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'F* =',& write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'F* =',&
math_transpose33(F_star_av) math_transpose33(F_star_av)
@ -832,7 +839,7 @@ program DAMASK_spectral
enddo ! end looping when convergency is achieved enddo ! end looping when convergency is achieved
print '(a)', '' print '(a)', ''
print '(a)', '==================================================================' print '(a)', '=================================================================='
if(err_f > err_div_tol .or. err_stress > err_stress_tol) then if(err_crit > err_div_tol .or. err_stress > err_stress_tol) then
print '(A,I5.5,A)', 'increment ', totalIncsCounter, ' NOT converged' print '(A,I5.5,A)', 'increment ', totalIncsCounter, ' NOT converged'
notConvergedCounter = notConvergedCounter + 1_pInt notConvergedCounter = notConvergedCounter + 1_pInt
else else
@ -869,7 +876,7 @@ program DAMASK_spectral
close(538) close(538)
call fftw_destroy_plan(plan_lambda); call fftw_destroy_plan(plan_correction) call fftw_destroy_plan(plan_lambda); call fftw_destroy_plan(plan_correction)
call quit(0_pInt) call quit(0_pInt)
end program DAMASK_spectral end program DAMASK_spectral_AL
!******************************************************************** !********************************************************************
! quit subroutine to satisfy IO_error ! quit subroutine to satisfy IO_error

46
code/FEsolving.f90
View File

@ -21,42 +21,47 @@
!############################################################## !##############################################################
module FEsolving module FEsolving
!############################################################## !##############################################################
use prec, only: pInt,pReal use prec, only: &
pInt, &
pReal
implicit none implicit none
integer(pInt) :: & private
integer(pInt), public :: &
cycleCounter = 0_pInt, & cycleCounter = 0_pInt, &
theInc = -1_pInt, & theInc = -1_pInt, &
restartInc = 1_pInt restartInc = 1_pInt
real(pReal) :: & real(pReal), public :: &
theTime = 0.0_pReal, & theTime = 0.0_pReal, &
theDelta = 0.0_pReal theDelta = 0.0_pReal
logical :: & logical, public :: &
lastIncConverged = .false., &
outdatedByNewInc = .false., &
outdatedFFN1 = .false., & outdatedFFN1 = .false., &
terminallyIll = .false., &
symmetricSolver = .false., & symmetricSolver = .false., &
parallelExecution = .true., &
restartWrite = .false., & restartWrite = .false., &
restartRead = .false., & restartRead = .false., &
lastMode = .true., & terminallyIll = .false., &
cutBack = .false. parallelExecution = .true., &
lastMode = .true.
integer(pInt), dimension(:,:), allocatable :: & integer(pInt), dimension(:,:), allocatable, public :: &
FEsolving_execIP FEsolving_execIP
integer(pInt), dimension(2) :: & integer(pInt), dimension(2), public :: &
FEsolving_execElem FEsolving_execElem
character(len=1024) :: & character(len=1024), public :: &
FEmodelGeometry FEmodelGeometry
logical, dimension(:,:), allocatable :: & logical, dimension(:,:), allocatable, public :: &
calcMode calcMode
logical, private :: &
lastIncConverged = .false., &
outdatedByNewInc = .false., &
cutBack = .false.
public :: FE_init public :: FE_init
contains contains
@ -68,21 +73,26 @@ contains
subroutine FE_init subroutine FE_init
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use debug, only: debug_what, & use debug, only: &
debug_what, &
debug_FEsolving, & debug_FEsolving, &
debug_levelBasic debug_levelBasic
use IO, only: IO_open_inputFile, &
use IO, only: &
IO_open_inputFile, &
IO_stringPos, & IO_stringPos, &
IO_stringValue, & IO_stringValue, &
IO_intValue, & IO_intValue, &
IO_lc, & IO_lc, &
IO_open_logFile, & IO_open_logFile, &
IO_warning IO_warning
use DAMASK_interface use DAMASK_interface
implicit none implicit none
integer(pInt), parameter :: fileunit = 222_pInt integer(pInt), parameter :: &
integer(pInt), parameter :: maxNchunks = 6_pInt fileunit = 222_pInt, &
maxNchunks = 6_pInt
integer :: i, start = 0, length ! is save for FE_init (only called once) integer :: i, start = 0, length ! is save for FE_init (only called once)
integer(pInt) :: j integer(pInt) :: j

2
code/Makefile
View File

@ -343,7 +343,7 @@ DAMASK_spectral_interface.o: DAMASK_spectral_interface.f90 \
prec.o: prec.f90 prec.o: prec.f90
%.o : %.f90 %.o : %.f90
$(PREFIX) $(COMPILERNAME) $(COMPILE_MAXOPTI) -c $< $(SUFFIX) $(PREFIX) $(COMPILERNAME) $(COMPILE) -c $< $(SUFFIX)
tidy: tidy:
rm -rf *.o rm -rf *.o

1
code/config/debug.config
View File

@ -12,6 +12,7 @@ crystallite # crystallite.f90 possible keys: basi
homogenization # homogenization_*.f90 possible keys: basic, extensive, selective homogenization # homogenization_*.f90 possible keys: basic, extensive, selective
CPFEM # CPFEM.f90 possible keys: basic, selective CPFEM # CPFEM.f90 possible keys: basic, selective
spectral # DAMASK_spectral.f90 possible keys: basic, fft, restart, divergence spectral # DAMASK_spectral.f90 possible keys: basic, fft, restart, divergence
abaqus # ABAQUS FEM solver possible keys: basic
# #
# Parameters for selective # Parameters for selective
element 1 # selected element for debugging (synonymous: "el", "e") element 1 # selected element for debugging (synonymous: "el", "e")

8
code/debug.f90
View File

@ -51,9 +51,10 @@ module debug
debug_crystallite = 8_pInt, & debug_crystallite = 8_pInt, &
debug_homogenization = 9_pInt, & debug_homogenization = 9_pInt, &
debug_CPFEM = 10_pInt, & debug_CPFEM = 10_pInt, &
debug_spectral = 11_pInt debug_spectral = 11_pInt, &
debug_abaqus = 12_pInt
integer(pInt), dimension(11+2), public :: & ! 11 for specific, and 2 for "all" and "other" integer(pInt), dimension(12+2), public :: & ! 11 for specific, and 2 for "all" and "other"
debug_what = 0_pInt debug_what = 0_pInt
integer(pInt), public :: & integer(pInt), public :: &
@ -193,6 +194,8 @@ subroutine debug_init
what = debug_CPFEM what = debug_CPFEM
case ('spectral') case ('spectral')
what = debug_spectral what = debug_spectral
case ('abaqus')
what = debug_abaqus
case ('all') case ('all')
what = 12_pInt what = 12_pInt
case ('other') case ('other')
@ -257,6 +260,7 @@ subroutine debug_init
if(i == debug_homogenization) write(6,'(a)') 'Homogenization debugging:' if(i == debug_homogenization) write(6,'(a)') 'Homogenization debugging:'
if(i == debug_CPFEM) write(6,'(a)') 'CPFEM debugging:' if(i == debug_CPFEM) write(6,'(a)') 'CPFEM debugging:'
if(i == debug_spectral) write(6,'(a)') 'Spectral solver debugging:' if(i == debug_spectral) write(6,'(a)') 'Spectral solver debugging:'
if(i == debug_abaqus) write(6,'(a)') 'ABAQUS FEM solver debugging:'
if(iand(debug_what(i),debug_levelBasic) /= 0) write(6,'(a)') ' basic' if(iand(debug_what(i),debug_levelBasic) /= 0) write(6,'(a)') ' basic'
if(iand(debug_what(i),debug_levelExtensive) /= 0) write(6,'(a)') ' extensive' if(iand(debug_what(i),debug_levelExtensive) /= 0) write(6,'(a)') ' extensive'

236
code/lattice.f90
View File

@ -29,74 +29,83 @@
module lattice module lattice
use prec, only: pReal,pInt use prec, only: pReal, &
pInt
implicit none implicit none
private
!************************************ !************************************
!* Lattice structures * !* Lattice structures *
!************************************ !************************************
integer(pInt) :: & integer(pInt), parameter, public :: &
lattice_Nhexagonal, & !> # of hexagonal lattice structure (from tag CoverA_ratio)
lattice_Nstructure !> # of lattice structures (1: fcc,2: bcc,3+: hexagonal)
integer(pInt), parameter :: &
lattice_maxNslipFamily = 5_pInt, & !> max # of slip system families over lattice structures lattice_maxNslipFamily = 5_pInt, & !> max # of slip system families over lattice structures
lattice_maxNtwinFamily = 4_pInt, & !> max # of twin system families over lattice structures lattice_maxNtwinFamily = 4_pInt, & !> max # of twin system families over lattice structures
lattice_maxNslip = 54_pInt, & !> max # of slip systems over lattice structures lattice_maxNslip = 54_pInt, & !> max # of slip systems over lattice structures
lattice_maxNtwin = 24_pInt, & !> max # of twin systems over lattice structures lattice_maxNtwin = 24_pInt, & !> max # of twin systems over lattice structures
lattice_maxNinteraction = 30_pInt !> max # of interaction types (in hardening matrix part) lattice_maxNinteraction = 30_pInt !> max # of interaction types (in hardening matrix part)
integer(pInt), pointer, dimension(:,:) :: & integer(pInt), allocatable, dimension(:,:), public :: &
interactionSlipSlip, &
interactionSlipTwin, &
interactionTwinSlip, &
interactionTwinTwin
real(pReal), allocatable, dimension(:,:,:,:) :: &
lattice_Sslip ! Schmid matrices, normal, shear direction and d x n of slip systems
real(pReal), allocatable, dimension(:,:,:) :: &
lattice_Sslip_v, &
lattice_sn, &
lattice_sd, &
lattice_st
! rotation and Schmid matrices, normal, shear direction and d x n of twin systems
real(pReal), allocatable, dimension(:,:,:,:) :: &
lattice_Qtwin, &
lattice_Stwin
real(pReal), allocatable, dimension(:,:,:) :: &
lattice_Stwin_v, &
lattice_tn, &
lattice_td, &
lattice_tt
real(pReal), allocatable, dimension(:,:) :: &
lattice_shearTwin !> characteristic twin shear
integer(pInt), allocatable, dimension(:,:) :: &
lattice_NslipSystem, & !> number of slip systems in each family lattice_NslipSystem, & !> number of slip systems in each family
lattice_NtwinSystem !> number of twin systems in each family lattice_NtwinSystem !> number of twin systems in each family
integer(pInt), allocatable, dimension(:,:,:) :: & integer(pInt), allocatable, dimension(:,:,:), public :: &
lattice_interactionSlipSlip, & !> interaction type between slip/slip lattice_interactionSlipSlip, & !> interaction type between slip/slip
lattice_interactionSlipTwin, & !> interaction type between slip/twin lattice_interactionSlipTwin, & !> interaction type between slip/twin
lattice_interactionTwinSlip, & !> interaction type between twin/slip lattice_interactionTwinSlip, & !> interaction type between twin/slip
lattice_interactionTwinTwin !> interaction type between twin/twin lattice_interactionTwinTwin !> interaction type between twin/twin
real(pReal), allocatable, dimension(:,:,:,:), public :: &
lattice_Sslip !> Schmid matrices, normal, shear direction and d x n of slip systems
real(pReal), allocatable, dimension(:,:,:), public :: &
lattice_Sslip_v, &
lattice_sn, &
lattice_sd, &
lattice_st
! rotation and Schmid matrices, normal, shear direction and d x n of twin systems
real(pReal), allocatable, dimension(:,:,:,:), public :: &
lattice_Stwin, &
lattice_Qtwin
real(pReal), allocatable, dimension(:,:,:), public :: &
lattice_Stwin_v, &
lattice_tn, &
lattice_td, &
lattice_tt
real(pReal), allocatable, dimension(:,:), public :: &
lattice_shearTwin !> characteristic twin shear
integer(pInt), private :: &
lattice_Nhexagonal, & !> # of hexagonal lattice structure (from tag CoverA_ratio)
lattice_Nstructure !> # of lattice structures (1: fcc,2: bcc,3+: hexagonal)
integer(pInt), dimension(:,:), pointer, private :: &
interactionSlipSlip, &
interactionSlipTwin, &
interactionTwinSlip, &
interactionTwinTwin
!============================== fcc (1) ================================= !============================== fcc (1) =================================
integer(pInt), parameter, dimension(lattice_maxNslipFamily) :: lattice_fcc_NslipSystem = int([12, 0, 0, 0, 0],pInt) integer(pInt), dimension(lattice_maxNslipFamily), parameter, private :: &
integer(pInt), parameter, dimension(lattice_maxNtwinFamily) :: lattice_fcc_NtwinSystem = int([12, 0, 0, 0],pInt) lattice_fcc_NslipSystem = int([12, 0, 0, 0, 0],pInt)
integer(pInt), parameter :: lattice_fcc_Nslip = 12_pInt ! sum(lattice_fcc_NslipSystem)
integer(pInt), parameter :: lattice_fcc_Ntwin = 12_pInt ! sum(lattice_fcc_NtwinSystem)
integer(pInt) :: lattice_fcc_Nstructure = 0_pInt
real(pReal), dimension(3+3,lattice_fcc_Nslip), parameter :: lattice_fcc_systemSlip = & integer(pInt), dimension(lattice_maxNtwinFamily), parameter, private :: &
reshape(real([& lattice_fcc_NtwinSystem = int([12, 0, 0, 0],pInt)
integer(pInt), parameter, private :: &
lattice_fcc_Nslip = 12_pInt, & ! sum(lattice_fcc_NslipSystem)
lattice_fcc_Ntwin = 12_pInt ! sum(lattice_fcc_NtwinSystem)
integer(pInt), private :: &
lattice_fcc_Nstructure = 0_pInt
real(pReal), dimension(3+3,lattice_fcc_Nslip), parameter, private :: &
lattice_fcc_systemSlip = reshape(real([&
! Slip system <110>{111} Sorted according to Eisenlohr & Hantcherli ! Slip system <110>{111} Sorted according to Eisenlohr & Hantcherli
0, 1,-1, 1, 1, 1, & 0, 1,-1, 1, 1, 1, &
-1, 0, 1, 1, 1, 1, & -1, 0, 1, 1, 1, 1, &
@ -112,8 +121,8 @@ module lattice
-1,-1, 0, -1, 1,-1 & -1,-1, 0, -1, 1,-1 &
],pReal),[ 3_pInt + 3_pInt,lattice_fcc_Nslip]) ],pReal),[ 3_pInt + 3_pInt,lattice_fcc_Nslip])
real(pReal), dimension(3+3,lattice_fcc_Ntwin), parameter :: lattice_fcc_systemTwin = & real(pReal), dimension(3+3,lattice_fcc_Ntwin), parameter, private :: &
reshape(real( [& lattice_fcc_systemTwin = reshape(real( [&
! Twin system <112>{111} Sorted according to Eisenlohr & Hantcherli ! Twin system <112>{111} Sorted according to Eisenlohr & Hantcherli
-2, 1, 1, 1, 1, 1, & -2, 1, 1, 1, 1, 1, &
1,-2, 1, 1, 1, 1, & 1,-2, 1, 1, 1, 1, &
@ -129,8 +138,8 @@ module lattice
-1, 1, 2, -1, 1,-1 & -1, 1, 2, -1, 1,-1 &
],pReal),[ 3_pInt + 3_pInt ,lattice_fcc_Ntwin]) ],pReal),[ 3_pInt + 3_pInt ,lattice_fcc_Ntwin])
real(pReal), dimension(lattice_fcc_Ntwin), parameter :: lattice_fcc_shearTwin = & real(pReal), dimension(lattice_fcc_Ntwin), parameter, private :: &
reshape([& lattice_fcc_shearTwin = reshape([&
! Twin system <112>{111} Sorted according to Eisenlohr & Hantcherli ! Twin system <112>{111} Sorted according to Eisenlohr & Hantcherli
0.7071067812_pReal, & 0.7071067812_pReal, &
0.7071067812_pReal, & 0.7071067812_pReal, &
@ -146,8 +155,8 @@ module lattice
0.7071067812_pReal & 0.7071067812_pReal &
],[lattice_fcc_Ntwin]) ],[lattice_fcc_Ntwin])
integer(pInt), target, dimension(lattice_fcc_Nslip,lattice_fcc_Nslip) :: lattice_fcc_interactionSlipSlip = & integer(pInt), dimension(lattice_fcc_Nslip,lattice_fcc_Nslip), target, private :: &
reshape(int( [& lattice_fcc_interactionSlipSlip = reshape(int( [&
! Interaction types ! Interaction types
! 1 --- self interaction ! 1 --- self interaction
! 2 --- coplanar interaction ! 2 --- coplanar interaction
@ -169,8 +178,8 @@ module lattice
6,5,4,5,6,4,5,5,3,2,2,1 & 6,5,4,5,6,4,5,5,3,2,2,1 &
],pInt),[lattice_fcc_Nslip,lattice_fcc_Nslip]) ],pInt),[lattice_fcc_Nslip,lattice_fcc_Nslip])
integer(pInt), target, dimension(lattice_fcc_Ntwin,lattice_fcc_Nslip) :: lattice_fcc_interactionSlipTwin = & integer(pInt), dimension(lattice_fcc_Ntwin,lattice_fcc_Nslip), target, private :: &
reshape(int( [& lattice_fcc_interactionSlipTwin = reshape(int( [&
1,1,1,2,2,1,1,2,2,2,1,2, & 1,1,1,2,2,1,1,2,2,2,1,2, &
1,1,1,2,2,1,1,2,2,2,1,2, & 1,1,1,2,2,1,1,2,2,2,1,2, &
1,1,1,2,2,1,1,2,2,2,1,2, & 1,1,1,2,2,1,1,2,2,2,1,2, &
@ -185,10 +194,11 @@ module lattice
2,1,2,1,2,2,2,2,1,1,1,1 & 2,1,2,1,2,2,2,2,1,1,1,1 &
],pInt),[lattice_fcc_Ntwin,lattice_fcc_Nslip]) ],pInt),[lattice_fcc_Ntwin,lattice_fcc_Nslip])
integer(pInt), target, dimension(lattice_fcc_Nslip,lattice_fcc_Ntwin) :: lattice_fcc_interactionTwinSlip = 0_pInt integer(pInt), dimension(lattice_fcc_Nslip,lattice_fcc_Ntwin), target, private :: &
lattice_fcc_interactionTwinSlip = 0_pInt
integer(pInt), target, dimension(lattice_fcc_Ntwin,lattice_fcc_Ntwin) :: lattice_fcc_interactionTwinTwin = & integer(pInt), dimension(lattice_fcc_Ntwin,lattice_fcc_Ntwin), target, private :: &
reshape(int( [& lattice_fcc_interactionTwinTwin = reshape(int( [&
1,1,1,2,2,2,2,2,2,2,2,2, & 1,1,1,2,2,2,2,2,2,2,2,2, &
1,1,1,2,2,2,2,2,2,2,2,2, & 1,1,1,2,2,2,2,2,2,2,2,2, &
1,1,1,2,2,2,2,2,2,2,2,2, & 1,1,1,2,2,2,2,2,2,2,2,2, &
@ -206,14 +216,23 @@ module lattice
!============================== bcc (2) ================================= !============================== bcc (2) =================================
integer(pInt), parameter, dimension(lattice_maxNslipFamily) :: lattice_bcc_NslipSystem = int([ 12,12,24, 0, 0], pInt) integer(pInt), dimension(lattice_maxNslipFamily), parameter, private :: &
integer(pInt), parameter, dimension(lattice_maxNtwinFamily) :: lattice_bcc_NtwinSystem = int([ 12, 0, 0, 0], pInt) lattice_bcc_NslipSystem = int([ 12,12,24, 0, 0], pInt)
integer(pInt), parameter :: lattice_bcc_Nslip = 48_pInt ! sum(lattice_bcc_NslipSystem)
integer(pInt), parameter :: lattice_bcc_Ntwin = 12_pInt ! sum(lattice_bcc_NtwinSystem)
integer(pInt) :: lattice_bcc_Nstructure = 0_pInt
real(pReal), dimension(3+3,lattice_bcc_Nslip), parameter :: lattice_bcc_systemSlip = & integer(pInt), dimension(lattice_maxNtwinFamily), parameter, private :: &
reshape(real([& lattice_bcc_NtwinSystem = int([ 12, 0, 0, 0], pInt)
integer(pInt), parameter, private :: &
lattice_bcc_Nslip = 48_pInt ! sum(lattice_bcc_NslipSystem)
integer(pInt), parameter, private :: &
lattice_bcc_Ntwin = 12_pInt ! sum(lattice_bcc_NtwinSystem)
integer(pInt), private :: &
lattice_bcc_Nstructure = 0_pInt
real(pReal), dimension(3+3,lattice_bcc_Nslip), parameter, private :: &
lattice_bcc_systemSlip = reshape(real([&
! Slip system <111>{110} meaningful sorting? ! Slip system <111>{110} meaningful sorting?
1,-1, 1, 0, 1, 1, & 1,-1, 1, 0, 1, 1, &
-1,-1, 1, 0, 1, 1, & -1,-1, 1, 0, 1, 1, &
@ -269,8 +288,8 @@ module lattice
! twin system <111>{112} ! twin system <111>{112}
! MISSING: not implemented yet -- now dummy copy from fcc !! ! MISSING: not implemented yet -- now dummy copy from fcc !!
real(pReal), dimension(3+3,lattice_bcc_Ntwin), parameter :: lattice_bcc_systemTwin = & real(pReal), dimension(3+3,lattice_bcc_Ntwin), parameter, private :: &
reshape(real([& lattice_bcc_systemTwin = reshape(real([&
! Twin system <112>{111} Sorted according to Eisenlohr & Hantcherli ! Twin system <112>{111} Sorted according to Eisenlohr & Hantcherli
-2, 1, 1, 1, 1, 1, & -2, 1, 1, 1, 1, 1, &
1,-2, 1, 1, 1, 1, & 1,-2, 1, 1, 1, 1, &
@ -286,8 +305,8 @@ module lattice
-1, 1, 2, -1, 1,-1 & -1, 1, 2, -1, 1,-1 &
],pReal),[ 3_pInt + 3_pInt,lattice_bcc_Ntwin]) ],pReal),[ 3_pInt + 3_pInt,lattice_bcc_Ntwin])
real(pReal), dimension(lattice_bcc_Ntwin), parameter :: lattice_bcc_shearTwin = & real(pReal), dimension(lattice_bcc_Ntwin), parameter, private :: &
reshape([& lattice_bcc_shearTwin = reshape([&
! Twin system {111}<112> just a dummy ! Twin system {111}<112> just a dummy
0.123_pReal, & 0.123_pReal, &
0.123_pReal, & 0.123_pReal, &
@ -304,8 +323,8 @@ module lattice
],[lattice_bcc_Ntwin]) ],[lattice_bcc_Ntwin])
!*** slip--slip interactions for BCC structures (2) *** !*** slip--slip interactions for BCC structures (2) ***
integer(pInt), target, dimension(lattice_bcc_Nslip,lattice_bcc_Nslip) :: lattice_bcc_interactionSlipSlip = & integer(pInt), dimension(lattice_bcc_Nslip,lattice_bcc_Nslip), target, private :: &
reshape(int( [& lattice_bcc_interactionSlipSlip = reshape(int( [&
1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2, & 1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2, &
2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2, & 2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2, &
2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2, & 2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2, &
@ -358,8 +377,8 @@ module lattice
!*** slip--twin interactions for BCC structures (2) *** !*** slip--twin interactions for BCC structures (2) ***
! MISSING: not implemented yet ! MISSING: not implemented yet
integer(pInt), target, dimension(lattice_bcc_Ntwin,lattice_bcc_Nslip) :: lattice_bcc_interactionSlipTwin = & integer(pInt), dimension(lattice_bcc_Ntwin,lattice_bcc_Nslip), target, private :: &
reshape(int( [& lattice_bcc_interactionSlipTwin = reshape(int( [&
0,0,0,0,0,0,0,0,0,0,0,0, & 0,0,0,0,0,0,0,0,0,0,0,0, &
0,0,0,0,0,0,0,0,0,0,0,0, & 0,0,0,0,0,0,0,0,0,0,0,0, &
0,0,0,0,0,0,0,0,0,0,0,0, & 0,0,0,0,0,0,0,0,0,0,0,0, &
@ -413,8 +432,8 @@ module lattice
!*** twin--slip interactions for BCC structures (2) *** !*** twin--slip interactions for BCC structures (2) ***
! MISSING: not implemented yet ! MISSING: not implemented yet
integer(pInt), target, dimension(lattice_bcc_Nslip,lattice_bcc_Ntwin) :: lattice_bcc_interactionTwinSlip = & integer(pInt), dimension(lattice_bcc_Nslip,lattice_bcc_Ntwin), target, private :: &
reshape(int( [& lattice_bcc_interactionTwinSlip = reshape(int( [&
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, & 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, &
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, & 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, &
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, & 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, &
@ -431,8 +450,8 @@ module lattice
!*** twin-twin interactions for BCC structures (2) *** !*** twin-twin interactions for BCC structures (2) ***
! MISSING: not implemented yet ! MISSING: not implemented yet
integer(pInt), target, dimension(lattice_bcc_Ntwin,lattice_bcc_Ntwin) :: lattice_bcc_interactionTwinTwin = & integer(pInt), dimension(lattice_bcc_Ntwin,lattice_bcc_Ntwin), target, private :: &
reshape(int( [& lattice_bcc_interactionTwinTwin = reshape(int( [&
0,0,0,0,0,0,0,0,0,0,0,0, & 0,0,0,0,0,0,0,0,0,0,0,0, &
0,0,0,0,0,0,0,0,0,0,0,0, & 0,0,0,0,0,0,0,0,0,0,0,0, &
0,0,0,0,0,0,0,0,0,0,0,0, & 0,0,0,0,0,0,0,0,0,0,0,0, &
@ -450,15 +469,24 @@ module lattice
!============================== hex (3+) ================================= !============================== hex (3+) =================================
integer(pInt), parameter, dimension(lattice_maxNslipFamily) :: lattice_hex_NslipSystem = int([ 3, 3, 6,12, 6],pInt) integer(pInt), dimension(lattice_maxNslipFamily), parameter, private :: &
integer(pInt), parameter, dimension(lattice_maxNtwinFamily) :: lattice_hex_NtwinSystem = int([ 6, 6, 6, 6],pInt) lattice_hex_NslipSystem = int([ 3, 3, 6,12, 6],pInt)
integer(pInt), parameter :: lattice_hex_Nslip = 30_pInt ! sum(lattice_hex_NslipSystem)
integer(pInt), parameter :: lattice_hex_Ntwin = 24_pInt ! sum(lattice_hex_NtwinSystem) integer(pInt), dimension(lattice_maxNtwinFamily), parameter, private :: &
integer(pInt) :: lattice_hex_Nstructure = 0_pInt lattice_hex_NtwinSystem = int([ 6, 6, 6, 6],pInt)
integer(pInt), parameter , private :: &
lattice_hex_Nslip = 30_pInt ! sum(lattice_hex_NslipSystem)
integer(pInt), parameter, private :: &
lattice_hex_Ntwin = 24_pInt ! sum(lattice_hex_NtwinSystem)
integer(pInt), private :: &
lattice_hex_Nstructure = 0_pInt
!* sorted by A. Alankar & P. Eisenlohr !* sorted by A. Alankar & P. Eisenlohr
real(pReal), dimension(4+4,lattice_hex_Nslip), parameter :: lattice_hex_systemSlip = & real(pReal), dimension(4+4,lattice_hex_Nslip), parameter, private :: &
reshape(real([& lattice_hex_systemSlip = reshape(real([&
! Basal systems <1120>{0001} (independent of c/a-ratio, Bravais notation (4 coordinate base)) ! Basal systems <1120>{0001} (independent of c/a-ratio, Bravais notation (4 coordinate base))
2, -1, -1, 0, 0, 0, 0, 1, & 2, -1, -1, 0, 0, 0, 0, 1, &
-1, 2, -1, 0, 0, 0, 0, 1, & -1, 2, -1, 0, 0, 0, 0, 1, &
@ -496,8 +524,8 @@ module lattice
1, -2, 1, -3, 1, -2, 1, 2 & 1, -2, 1, -3, 1, -2, 1, 2 &
],pReal),[ 4_pInt + 4_pInt,lattice_hex_Nslip]) ],pReal),[ 4_pInt + 4_pInt,lattice_hex_Nslip])
real(pReal), dimension(4+4,lattice_hex_Ntwin), parameter :: lattice_hex_systemTwin = & real(pReal), dimension(4+4,lattice_hex_Ntwin), parameter, private :: &
reshape(real([& lattice_hex_systemTwin = reshape(real([&
0, 1, -1, 1, 0, -1, 1, 2, & ! <-10.1>{10.2} shear = (3-(c/a)^2)/(sqrt(3) c/a) 0, 1, -1, 1, 0, -1, 1, 2, & ! <-10.1>{10.2} shear = (3-(c/a)^2)/(sqrt(3) c/a)
-1, 1, 0, 1, 1, -1, 0, 2, & -1, 1, 0, 1, 1, -1, 0, 2, &
-1, 0, 1, 1, 1, 0, -1, 2, & !! -1, 0, 1, 1, 1, 0, -1, 2, & !!
@ -524,8 +552,8 @@ module lattice
-1, 1, 0, -2, -1, 1, 0, 1 & -1, 1, 0, -2, -1, 1, 0, 1 &
],pReal),[ 4_pInt + 4_pInt ,lattice_hex_Ntwin]) !* Sort? Numbering of twin system follows Prof. Tom Bieler's scheme (to be consistent with his work); but numbering in data was restarted from 1 & ],pReal),[ 4_pInt + 4_pInt ,lattice_hex_Ntwin]) !* Sort? Numbering of twin system follows Prof. Tom Bieler's scheme (to be consistent with his work); but numbering in data was restarted from 1 &
integer(pInt), dimension(lattice_hex_Ntwin), parameter :: lattice_hex_shearTwin = & ! indicator to formula further below integer(pInt), dimension(lattice_hex_Ntwin), parameter, private :: &
reshape(int( [& lattice_hex_shearTwin = reshape(int( [& ! indicator to formula further below
1, & ! {10.2}<-10.1> 1, & ! {10.2}<-10.1>
1, & 1, &
1, & 1, &
@ -558,8 +586,8 @@ module lattice
!* 3. twin-twin interaction - 20 types !* 3. twin-twin interaction - 20 types
!* 4. twin-slip interaction - 16 types !* 4. twin-slip interaction - 16 types
integer(pInt), target, dimension(lattice_hex_Nslip,lattice_hex_Nslip) :: lattice_hex_interactionSlipSlip = & integer(pInt), dimension(lattice_hex_Nslip,lattice_hex_Nslip), target, private :: &
reshape(int( [& lattice_hex_interactionSlipSlip = reshape(int( [&
1, 6, 6, 11,11,11, 15,15,15,15,15,15, 18,18,18,18,18,18,18,18,18,18,18,18, 20,20,20,20,20,20, & 1, 6, 6, 11,11,11, 15,15,15,15,15,15, 18,18,18,18,18,18,18,18,18,18,18,18, 20,20,20,20,20,20, &
6, 1, 6, 11,11,11, 15,15,15,15,15,15, 18,18,18,18,18,18,18,18,18,18,18,18, 20,20,20,20,20,20, & 6, 1, 6, 11,11,11, 15,15,15,15,15,15, 18,18,18,18,18,18,18,18,18,18,18,18, 20,20,20,20,20,20, &
6, 6, 1, 11,11,11, 15,15,15,15,15,15, 18,18,18,18,18,18,18,18,18,18,18,18, 20,20,20,20,20,20, & 6, 6, 1, 11,11,11, 15,15,15,15,15,15, 18,18,18,18,18,18,18,18,18,18,18,18, 20,20,20,20,20,20, &
@ -597,8 +625,8 @@ module lattice
],pInt),[lattice_hex_Nslip,lattice_hex_Nslip]) ],pInt),[lattice_hex_Nslip,lattice_hex_Nslip])
!* isotropic interaction at the moment !* isotropic interaction at the moment
integer(pInt), target, dimension(lattice_hex_Ntwin,lattice_hex_Nslip) :: lattice_hex_interactionSlipTwin = & integer(pInt), dimension(lattice_hex_Ntwin,lattice_hex_Nslip), target, private :: &
reshape(int( [& lattice_hex_interactionSlipTwin = reshape(int( [&
1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & ! --> twin 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & ! --> twin
1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & ! | 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & ! |
1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & ! | 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & ! |
@ -636,8 +664,8 @@ module lattice
],pInt),[lattice_hex_Ntwin,lattice_hex_Nslip]) ],pInt),[lattice_hex_Ntwin,lattice_hex_Nslip])
!* isotropic interaction at the moment !* isotropic interaction at the moment
integer(pInt), target, dimension(lattice_hex_Nslip,lattice_hex_Ntwin) :: lattice_hex_interactionTwinSlip = & integer(pInt), dimension(lattice_hex_Nslip,lattice_hex_Ntwin), target, private :: &
reshape(int( [& lattice_hex_interactionTwinSlip = reshape(int( [&
1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9, 13,13,13,13,13,13,13,13,13,13,13,13, 17,17,17,17,17,17, & ! --> slip 1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9, 13,13,13,13,13,13,13,13,13,13,13,13, 17,17,17,17,17,17, & ! --> slip
1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9, 13,13,13,13,13,13,13,13,13,13,13,13, 17,17,17,17,17,17, & ! | 1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9, 13,13,13,13,13,13,13,13,13,13,13,13, 17,17,17,17,17,17, & ! |
1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9, 13,13,13,13,13,13,13,13,13,13,13,13, 17,17,17,17,17,17, & ! | 1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9, 13,13,13,13,13,13,13,13,13,13,13,13, 17,17,17,17,17,17, & ! |
@ -668,8 +696,8 @@ module lattice
],pInt),[lattice_hex_Nslip,lattice_hex_Ntwin]) ],pInt),[lattice_hex_Nslip,lattice_hex_Ntwin])
integer(pInt), target, dimension(lattice_hex_Ntwin,lattice_hex_Ntwin) :: lattice_hex_interactionTwinTwin = & integer(pInt), dimension(lattice_hex_Ntwin,lattice_hex_Ntwin), target, private :: &
reshape(int( [& lattice_hex_interactionTwinTwin = reshape(int( [&
1, 5, 5, 5, 5, 5, 9, 9, 9, 9, 9, 9, 12,12,12,12,12,12, 14,14,14,14,14,14, & 1, 5, 5, 5, 5, 5, 9, 9, 9, 9, 9, 9, 12,12,12,12,12,12, 14,14,14,14,14,14, &
5, 1, 5, 5, 5, 5, 9, 9, 9, 9, 9, 9, 12,12,12,12,12,12, 14,14,14,14,14,14, & 5, 1, 5, 5, 5, 5, 9, 9, 9, 9, 9, 9, 12,12,12,12,12,12, 14,14,14,14,14,14, &
5, 5, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9, 12,12,12,12,12,12, 14,14,14,14,14,14, & 5, 5, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9, 12,12,12,12,12,12, 14,14,14,14,14,14, &
@ -699,12 +727,12 @@ module lattice
20,20,20,20,20,20, 19,19,19,19,19,19, 17,17,17,17,17,17, 8, 8, 8, 8, 8, 4 & 20,20,20,20,20,20, 19,19,19,19,19,19, 17,17,17,17,17,17, 8, 8, 8, 8, 8, 4 &
],pInt),[lattice_hex_Ntwin,lattice_hex_Ntwin]) ],pInt),[lattice_hex_Ntwin,lattice_hex_Ntwin])
public :: &
lattice_init, &
lattice_initializeStructure, &
lattice_symmetryType
CONTAINS contains
!****************************************
!* - lattice_init
!* - lattice_initializeStructure
!****************************************
integer(pInt) pure function lattice_symmetryType(structID) integer(pInt) pure function lattice_symmetryType(structID)
!************************************** !**************************************
@ -792,10 +820,14 @@ subroutine lattice_init
allocate(lattice_NslipSystem(lattice_maxNslipFamily,lattice_Nstructure)); lattice_NslipSystem = 0_pInt allocate(lattice_NslipSystem(lattice_maxNslipFamily,lattice_Nstructure)); lattice_NslipSystem = 0_pInt
allocate(lattice_NtwinSystem(lattice_maxNtwinFamily,lattice_Nstructure)); lattice_NtwinSystem = 0_pInt allocate(lattice_NtwinSystem(lattice_maxNtwinFamily,lattice_Nstructure)); lattice_NtwinSystem = 0_pInt
allocate(lattice_interactionSlipSlip(lattice_maxNslip,lattice_maxNslip,lattice_Nstructure)); lattice_interactionSlipSlip = 0_pInt ! other:me allocate(lattice_interactionSlipSlip(lattice_maxNslip,lattice_maxNslip,lattice_Nstructure))
allocate(lattice_interactionSlipTwin(lattice_maxNtwin,lattice_maxNslip,lattice_Nstructure)); lattice_interactionSlipTwin = 0_pInt ! other:me lattice_interactionSlipSlip = 0_pInt ! other:me
allocate(lattice_interactionTwinSlip(lattice_maxNslip,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionTwinSlip = 0_pInt ! other:me allocate(lattice_interactionSlipTwin(lattice_maxNtwin,lattice_maxNslip,lattice_Nstructure))
allocate(lattice_interactionTwinTwin(lattice_maxNtwin,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionTwinTwin = 0_pInt ! other:me lattice_interactionSlipTwin = 0_pInt ! other:me
allocate(lattice_interactionTwinSlip(lattice_maxNslip,lattice_maxNtwin,lattice_Nstructure))
lattice_interactionTwinSlip = 0_pInt ! other:me
allocate(lattice_interactionTwinTwin(lattice_maxNtwin,lattice_maxNtwin,lattice_Nstructure))
lattice_interactionTwinTwin = 0_pInt ! other:me
end subroutine lattice_init end subroutine lattice_init

58
code/numerics.f90
View File

@ -288,61 +288,61 @@ subroutine numerics_init
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(a24,1x,e8.1)') ' relevantStrain: ',relevantStrain write(6,'(a24,1x,es8.1)') ' relevantStrain: ',relevantStrain
write(6,'(a24,1x,e8.1)') ' defgradTolerance: ',defgradTolerance write(6,'(a24,1x,es8.1)') ' defgradTolerance: ',defgradTolerance
write(6,'(a24,1x,i8)') ' iJacoStiffness: ',iJacoStiffness write(6,'(a24,1x,i8)') ' iJacoStiffness: ',iJacoStiffness
write(6,'(a24,1x,i8)') ' iJacoLpresiduum: ',iJacoLpresiduum write(6,'(a24,1x,i8)') ' iJacoLpresiduum: ',iJacoLpresiduum
write(6,'(a24,1x,e8.1)') ' pert_Fg: ',pert_Fg write(6,'(a24,1x,es8.1)') ' pert_Fg: ',pert_Fg
write(6,'(a24,1x,i8)') ' pert_method: ',pert_method write(6,'(a24,1x,i8)') ' pert_method: ',pert_method
write(6,'(a24,1x,i8)') ' nCryst: ',nCryst write(6,'(a24,1x,i8)') ' nCryst: ',nCryst
write(6,'(a24,1x,e8.1)') ' subStepMinCryst: ',subStepMinCryst write(6,'(a24,1x,es8.1)') ' subStepMinCryst: ',subStepMinCryst
write(6,'(a24,1x,e8.1)') ' subStepSizeCryst: ',subStepSizeCryst write(6,'(a24,1x,es8.1)') ' subStepSizeCryst: ',subStepSizeCryst
write(6,'(a24,1x,e8.1)') ' stepIncreaseCryst: ',stepIncreaseCryst write(6,'(a24,1x,es8.1)') ' stepIncreaseCryst: ',stepIncreaseCryst
write(6,'(a24,1x,i8)') ' nState: ',nState write(6,'(a24,1x,i8)') ' nState: ',nState
write(6,'(a24,1x,i8)') ' nStress: ',nStress write(6,'(a24,1x,i8)') ' nStress: ',nStress
write(6,'(a24,1x,e8.1)') ' rTol_crystalliteState: ',rTol_crystalliteState write(6,'(a24,1x,es8.1)') ' rTol_crystalliteState: ',rTol_crystalliteState
write(6,'(a24,1x,e8.1)') ' rTol_crystalliteTemp: ',rTol_crystalliteTemperature write(6,'(a24,1x,es8.1)') ' rTol_crystalliteTemp: ',rTol_crystalliteTemperature
write(6,'(a24,1x,e8.1)') ' rTol_crystalliteStress: ',rTol_crystalliteStress write(6,'(a24,1x,es8.1)') ' rTol_crystalliteStress: ',rTol_crystalliteStress
write(6,'(a24,1x,e8.1)') ' aTol_crystalliteStress: ',aTol_crystalliteStress write(6,'(a24,1x,es8.1)') ' aTol_crystalliteStress: ',aTol_crystalliteStress
write(6,'(a24,2(1x,i8))') ' integrator: ',numerics_integrator write(6,'(a24,2(1x,i8))') ' integrator: ',numerics_integrator
write(6,'(a24,1x,L8,/)') ' analytic Jacobian: ',analyticJaco write(6,'(a24,1x,L8,/)') ' analytic Jacobian: ',analyticJaco
write(6,'(a24,1x,i8)') ' nHomog: ',nHomog write(6,'(a24,1x,i8)') ' nHomog: ',nHomog
write(6,'(a24,1x,e8.1)') ' subStepMinHomog: ',subStepMinHomog write(6,'(a24,1x,es8.1)') ' subStepMinHomog: ',subStepMinHomog
write(6,'(a24,1x,e8.1)') ' subStepSizeHomog: ',subStepSizeHomog write(6,'(a24,1x,es8.1)') ' subStepSizeHomog: ',subStepSizeHomog
write(6,'(a24,1x,e8.1)') ' stepIncreaseHomog: ',stepIncreaseHomog write(6,'(a24,1x,es8.1)') ' stepIncreaseHomog: ',stepIncreaseHomog
write(6,'(a24,1x,i8,/)') ' nMPstate: ',nMPstate write(6,'(a24,1x,i8,/)') ' nMPstate: ',nMPstate
!* RGC parameters !* RGC parameters
write(6,'(a24,1x,e8.1)') ' aTol_RGC: ',absTol_RGC write(6,'(a24,1x,es8.1)') ' aTol_RGC: ',absTol_RGC
write(6,'(a24,1x,e8.1)') ' rTol_RGC: ',relTol_RGC write(6,'(a24,1x,es8.1)') ' rTol_RGC: ',relTol_RGC
write(6,'(a24,1x,e8.1)') ' aMax_RGC: ',absMax_RGC write(6,'(a24,1x,es8.1)') ' aMax_RGC: ',absMax_RGC
write(6,'(a24,1x,e8.1)') ' rMax_RGC: ',relMax_RGC write(6,'(a24,1x,es8.1)') ' rMax_RGC: ',relMax_RGC
write(6,'(a24,1x,e8.1)') ' perturbPenalty_RGC: ',pPert_RGC write(6,'(a24,1x,es8.1)') ' perturbPenalty_RGC: ',pPert_RGC
write(6,'(a24,1x,e8.1)') ' relevantMismatch_RGC: ',xSmoo_RGC write(6,'(a24,1x,es8.1)') ' relevantMismatch_RGC: ',xSmoo_RGC
write(6,'(a24,1x,e8.1)') ' viscosityrate_RGC: ',viscPower_RGC write(6,'(a24,1x,es8.1)') ' viscosityrate_RGC: ',viscPower_RGC
write(6,'(a24,1x,e8.1)') ' viscositymodulus_RGC: ',viscModus_RGC write(6,'(a24,1x,es8.1)') ' viscositymodulus_RGC: ',viscModus_RGC
write(6,'(a24,1x,e8.1)') ' maxrelaxation_RGC: ',maxdRelax_RGC write(6,'(a24,1x,es8.1)') ' maxrelaxation_RGC: ',maxdRelax_RGC
write(6,'(a24,1x,e8.1)') ' maxVolDiscrepancy_RGC: ',maxVolDiscr_RGC write(6,'(a24,1x,es8.1)') ' maxVolDiscrepancy_RGC: ',maxVolDiscr_RGC
write(6,'(a24,1x,e8.1)') ' volDiscrepancyMod_RGC: ',volDiscrMod_RGC write(6,'(a24,1x,es8.1)') ' volDiscrepancyMod_RGC: ',volDiscrMod_RGC
write(6,'(a24,1x,e8.1,/)') ' discrepancyPower_RGC: ',volDiscrPow_RGC write(6,'(a24,1x,es8.1,/)') ' discrepancyPower_RGC: ',volDiscrPow_RGC
!* spectral parameters !* spectral parameters
write(6,'(a24,1x,e8.1)') ' err_div_tol: ',err_div_tol write(6,'(a24,1x,es8.1)') ' err_div_tol: ',err_div_tol
write(6,'(a24,1x,e8.1)') ' err_stress_tolrel: ',err_stress_tolrel write(6,'(a24,1x,es8.1)') ' err_stress_tolrel: ',err_stress_tolrel
write(6,'(a24,1x,i8)') ' itmax: ',itmax write(6,'(a24,1x,i8)') ' itmax: ',itmax
write(6,'(a24,1x,i8)') ' itmin: ',itmin write(6,'(a24,1x,i8)') ' itmin: ',itmin
write(6,'(a24,1x,L8)') ' memory_efficient: ',memory_efficient write(6,'(a24,1x,L8)') ' memory_efficient: ',memory_efficient
if(fftw_timelimit<0.0_pReal) then if(fftw_timelimit<0.0_pReal) then
write(6,'(a24,1x,L8)') ' fftw_timelimit: ',.false. write(6,'(a24,1x,L8)') ' fftw_timelimit: ',.false.
else else
write(6,'(a24,1x,e8.1)') ' fftw_timelimit: ',fftw_timelimit write(6,'(a24,1x,es8.1)') ' fftw_timelimit: ',fftw_timelimit
endif endif
write(6,'(a24,1x,a)') ' fftw_plan_mode: ',trim(fftw_plan_mode) write(6,'(a24,1x,a)') ' fftw_plan_mode: ',trim(fftw_plan_mode)
write(6,'(a24,1x,i8)') ' fftw_planner_flag: ',fftw_planner_flag write(6,'(a24,1x,i8)') ' fftw_planner_flag: ',fftw_planner_flag
write(6,'(a24,1x,e8.1)') ' rotation_tol: ',rotation_tol write(6,'(a24,1x,es8.1)') ' rotation_tol: ',rotation_tol
write(6,'(a24,1x,L8,/)') ' divergence_correction: ',divergence_correction write(6,'(a24,1x,L8,/)') ' divergence_correction: ',divergence_correction
write(6,'(a24,1x,L8,/)') ' update_gamma: ',update_gamma write(6,'(a24,1x,L8,/)') ' update_gamma: ',update_gamma