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

first try of implementing a rotation of the loadcase coordinate system. Seems to work so far for one loadcase and homogeneous materials.

This commit is contained in:
Martin Diehl 2011-10-24 18:26:34 +00:00
parent 1a68669ae9
commit c13aa2a829
4 changed files with 185 additions and 122 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

221
code/DAMASK_spectral.f90
View File

@ -47,14 +47,13 @@ program DAMASK_spectral
use DAMASK_interface use DAMASK_interface
use prec, only: pInt, pReal use prec, only: pInt, pReal
use IO use IO
use IFPORT, only: fseek, ftellI8
use debug, only: debug_Verbosity use debug, only: debug_Verbosity
use math use math
use mesh, only: mesh_ipCenterOfGravity use mesh, only: mesh_ipCenterOfGravity
use CPFEM, only: CPFEM_general, CPFEM_initAll use CPFEM, only: CPFEM_general, CPFEM_initAll
use numerics, only: err_div_tol, err_stress_tol, err_stress_tolrel,& use numerics, only: err_div_tol, err_stress_tol, err_stress_tolrel,&
relevantStrain, itmax, memory_efficient, DAMASK_NumThreadsInt,& relevantStrain, itmax, memory_efficient, DAMASK_NumThreadsInt,&
fftw_planner_flag, fftw_timelimit fftw_planner_flag, fftw_timelimit, tol_rotation
use homogenization, only: materialpoint_sizeResults, materialpoint_results use homogenization, only: materialpoint_sizeResults, materialpoint_results
!$ use OMP_LIB ! the openMP function library !$ use OMP_LIB ! the openMP function library
@ -62,10 +61,9 @@ program DAMASK_spectral
! variables to read from loadcase and geom file ! variables to read from loadcase and geom file
real(pReal), dimension(9) :: valueVector ! stores information temporarily from loadcase file real(pReal), dimension(9) :: valueVector ! stores information temporarily from loadcase file
integer(pInt), parameter :: maxNchunksLoadcase = & integer(pInt), parameter :: maxNchunksLoadcase = &
(1_pInt + 9_pInt)*2_pInt + & ! deformation and stress (1_pInt + 9_pInt)*3_pInt + & ! deformation, rotation, and stress
(1_pInt + 1_pInt)*4_pInt + & ! time, (log)incs, temp, and frequency (1_pInt + 1_pInt)*4_pInt + & ! time, (log)incs, temp, and frequency
1_pInt + 1_pInt +2_pInt*3_pInt + & ! rotation, keyword (deg,rad), axis + value 1_pInt ! dropguessing
1 ! dropguessing
integer(pInt), dimension (1 + maxNchunksLoadcase*2) :: posLoadcase integer(pInt), dimension (1 + maxNchunksLoadcase*2) :: posLoadcase
integer(pInt), parameter :: maxNchunksGeom = 7 ! 4 identifiers, 3 values integer(pInt), parameter :: maxNchunksGeom = 7 ! 4 identifiers, 3 values
integer(pInt), dimension (1 + maxNchunksGeom*2) :: posGeom integer(pInt), dimension (1 + maxNchunksGeom*2) :: posGeom
@ -74,6 +72,7 @@ program DAMASK_spectral
logical :: gotResolution, gotDimension, gotHomogenization logical :: gotResolution, gotDimension, gotHomogenization
! variables storing information from loadcase file ! variables storing information from loadcase file
!ToDo: create Data Structure loadcase
real(pReal), dimension (:,:,:), allocatable :: bc_deformation, & ! applied velocity gradient or time derivative of deformation gradient real(pReal), dimension (:,:,:), allocatable :: bc_deformation, & ! applied velocity gradient or time derivative of deformation gradient
bc_stress, & ! stress BC (if applicable) bc_stress, & ! stress BC (if applicable)
bc_rotation ! rotation of BC (if applicable) bc_rotation ! rotation of BC (if applicable)
@ -97,7 +96,8 @@ program DAMASK_spectral
! stress etc. ! stress etc.
real(pReal), dimension(3,3) :: pstress, pstress_av, defgrad_av, & real(pReal), dimension(3,3) :: pstress, pstress_av, defgrad_av, &
defgradAim, defgradAimOld, defgradAimCorr,& defgradAim, defgradAimOld, defgradAimCorr,&
mask_stress, mask_defgrad, fDot mask_stress, mask_defgrad, fDot, &
pstress_av_load, defgradAim_lab ! quantities rotated to other coordinate system
real(pReal), dimension(3,3,3,3) :: dPdF, c0_reference, c_current, s_prev, c_prev ! stiffness and compliance real(pReal), dimension(3,3,3,3) :: dPdF, c0_reference, c_current, s_prev, c_prev ! stiffness and compliance
real(pReal), dimension(6) :: cstress ! cauchy stress real(pReal), dimension(6) :: cstress ! cauchy stress
real(pReal), dimension(6,6) :: dsde ! small strain stiffness real(pReal), dimension(6,6) :: dsde ! small strain stiffness
@ -129,7 +129,6 @@ program DAMASK_spectral
integer(pInt) :: i, j, k, l, m, n, p integer(pInt) :: i, j, k, l, m, n, p
integer(pInt) :: N_Loadcases, loadcase, step, iter, ielem, CPFEM_mode, ierr, notConvergedCounter, writtenOutCounter integer(pInt) :: N_Loadcases, loadcase, step, iter, ielem, CPFEM_mode, ierr, notConvergedCounter, writtenOutCounter
logical errmatinv logical errmatinv
integer*8 :: incPosition
!Initializing !Initializing
!$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS !$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS
@ -149,6 +148,7 @@ program DAMASK_spectral
time = 0.0_pReal time = 0.0_pReal
notConvergedCounter = 0_pInt notConvergedCounter = 0_pInt
writtenOutCounter = 0_pInt
resolution = 1_pInt resolution = 1_pInt
geomdimension = 0.0_pReal geomdimension = 0.0_pReal
@ -160,6 +160,7 @@ program DAMASK_spectral
print '(a)', '******************************************************' print '(a)', '******************************************************'
print '(a,a)', 'Working Directory: ',trim(getSolverWorkingDirectoryName()) print '(a,a)', 'Working Directory: ',trim(getSolverWorkingDirectoryName())
print '(a,a)', 'Solver Job Name: ',trim(getSolverJobName()) print '(a,a)', 'Solver Job Name: ',trim(getSolverJobName())
print '(a)', '******************************************************'
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
if (.not. IO_open_file(myUnit,path)) call IO_error(30,ext_msg = trim(path)) if (.not. IO_open_file(myUnit,path)) call IO_error(30,ext_msg = trim(path))
@ -240,27 +241,22 @@ program DAMASK_spectral
bc_frequency(loadcase) = IO_intValue(line,posLoadcase,j+1) bc_frequency(loadcase) = IO_intValue(line,posLoadcase,j+1)
case('guessreset','dropguessing') case('guessreset','dropguessing')
bc_followFormerTrajectory(loadcase) = .false. ! do not continue to predict deformation along former trajectory bc_followFormerTrajectory(loadcase) = .false. ! do not continue to predict deformation along former trajectory
case('rotation','rot','euler') case('euler') ! rotation of loadcase given in euler angles
p = 0_pInt ! assuming values given in radians p = 0_pInt ! assuming values given in radians
k = 1_pInt ! assuming keyword indicating degree/radians l = 1_pInt ! assuming keyword indicating degree/radians
select case (IO_lc(IO_stringValue(line,posLoadcase,j+1))) select case (IO_lc(IO_stringValue(line,posLoadcase,j+1)))
case('deg','degree') case('deg','degree')
p = 1_pInt ! for conversion from degree to radian p = 1_pInt ! for conversion from degree to radian
case('rad','radian') case('rad','radian')
case default case default
k = 0_pInt ! imediately reading in angles, assuming radians l = 0_pInt ! imediately reading in angles, assuming radians
end select end select
do l = 0,4,2 ! looping to find keywords forall(k = 1:3) temp33_Real(k,1) = IO_floatValue(line,posLoadcase,j +l +k) * real(p,pReal) * inRad
select case (IO_lc(IO_stringValue(line,posLoadcase,j+1 +k +l)))
case('x')
temp33_Real(1,1) = IO_floatValue(line,posLoadcase,j+1 +k +l +1) * real(p,pReal) * inRad
case('y')
temp33_Real(2,1) = IO_floatValue(line,posLoadcase,j+1 +k +l +1) * real(p,pReal) * inRad
case('z')
temp33_Real(3,1) = IO_floatValue(line,posLoadcase,j+1 +k +l +1) * real(p,pReal) * inRad
end select
enddo
bc_rotation(:,:,loadcase) = math_EulerToR(temp33_Real(:,1)) bc_rotation(:,:,loadcase) = math_EulerToR(temp33_Real(:,1))
case('rotation','rot') ! assign values for the rotation of loadcase matrix
valueVector = 0.0_pReal
forall (k = 1:9) valueVector(k) = IO_floatValue(line,posLoadcase,j+k)
bc_rotation(:,:,loadcase) = math_plain9to33(valueVector)
end select end select
enddo; enddo enddo; enddo
@ -273,11 +269,6 @@ program DAMASK_spectral
spectralPictureMode = .false. spectralPictureMode = .false.
path = getModelName() path = getModelName()
!$OMP CRITICAL (write2out)
print '(a)', '******************************************************'
print '(a,a)', 'Geom File Name: ',trim(path)//'.geom'
print '(a)', '------------------------------------------------------'
!$OMP END CRITICAL (write2out)
if (.not. IO_open_file(myUnit,trim(path)//InputFileExtension))& if (.not. IO_open_file(myUnit,trim(path)//InputFileExtension))&
call IO_error(101,ext_msg = trim(path)//InputFileExtension) call IO_error(101,ext_msg = trim(path)//InputFileExtension)
@ -334,57 +325,6 @@ program DAMASK_spectral
mod(resolution(2),2_pInt)/=0_pInt .or.& mod(resolution(2),2_pInt)/=0_pInt .or.&
(mod(resolution(3),2_pInt)/=0_pInt .and. resolution(3)/= 1_pInt)) call IO_error(103) (mod(resolution(3),2_pInt)/=0_pInt .and. resolution(3)/= 1_pInt)) call IO_error(103)
!$OMP CRITICAL (write2out)
print '(a,/,i12,i12,i12)','resolution a b c:', resolution
print '(a,/,f12.5,f12.5,f12.5)','dimension x y z:', geomdimension
print '(a,i5)','homogenization: ',homog
print '(a,L)','spectralPictureMode: ',spectralPictureMode
print '(a)', '******************************************************'
print '(a,a)','Loadcase File Name: ',trim(getLoadcaseName())
if (bc_followFormerTrajectory(1)) then
call IO_warning(33) ! cannot guess along trajectory for first step of first loadcase
bc_followFormerTrajectory(1) = .false.
endif
! consistency checks and output of loadcase
do loadcase = 1, N_Loadcases
print '(a)', '------------------------------------------------------'
print '(a,i5)', 'Loadcase: ', loadcase
if (.not. bc_followFormerTrajectory(loadcase)) &
print '(a)', 'Drop Guessing Along Trajectory'
if (any(bc_mask(:,:,1,loadcase) .eqv. bc_mask(:,:,2,loadcase)))& ! exclusive or masking only
call IO_error(31,loadcase)
if (any(bc_mask(1:3,1:3,2,loadcase).and.transpose(bc_mask(1:3,1:3,2,loadcase)).and.& !checking if no rotation is allowed by stress BC
reshape((/.false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false./),(/3,3/))))&
call IO_error(38,loadcase)
if (bc_velGradApplied(loadcase)) then
do j = 1, 3
if (any(bc_mask(j,:,1,loadcase) .eqv. .true.) .and.&
any(bc_mask(j,:,1,loadcase) .eqv. .false.)) call IO_error(32,loadcase) ! each line should be either fully or not at all defined
enddo
print '(a,/,3(3(f12.6,x)/))','Velocity Gradient:', merge(math_transpose3x3(bc_deformation(:,:,loadcase)),&
reshape(spread(DAMASK_NaN,1,9),(/3,3/)),&
transpose(bc_mask(:,:,1,loadcase)))
else
print '(a,/,3(3(f12.6,x)/))','Change of Deformation Gradient:', merge(math_transpose3x3(bc_deformation(:,:,loadcase)),&
reshape(spread(DAMASK_NaN,1,9),(/3,3/)),&
transpose(bc_mask(:,:,1,loadcase)))
endif
print '(a,/,3(3(f12.6,x)/))','Stress Boundary Condition/MPa:',merge(math_transpose3x3(bc_stress(:,:,loadcase)),&
reshape(spread(DAMASK_NaN,1,9),(/3,3/)),&
transpose(bc_mask(:,:,2,loadcase)))*1e-6
if (any(bc_rotation(:,:,loadcase)/=math_I3)) &
print '(a,/,3(3(f12.6,x)/))','Rotation of BCs:',math_transpose3x3(bc_rotation(:,:,loadcase))
if (bc_timeIncrement(loadcase) < 0.0_pReal) call IO_error(34,loadcase) ! negative time increment
print '(a,f12.6)','Temperature: ',bc_temperature(loadcase)
print '(a,f12.6)','Time: ',bc_timeIncrement(loadcase)
if (bc_steps(loadcase) < 1_pInt) call IO_error(35,loadcase) ! non-positive increment count
print '(a,i5)','Increments: ',bc_steps(loadcase)
if (bc_frequency(loadcase) < 1_pInt) call IO_error(36,loadcase) ! non-positive result frequency
print '(a,i5)','Freq. of Output: ',bc_frequency(loadcase)
enddo
print '(a)', '******************************************************'
!$OMP END CRITICAL (write2out)
allocate (defgrad ( resolution(1),resolution(2),resolution(3),3,3)); defgrad = 0.0_pReal allocate (defgrad ( resolution(1),resolution(2),resolution(3),3,3)); defgrad = 0.0_pReal
allocate (defgradold ( resolution(1),resolution(2),resolution(3),3,3)); defgradold = 0.0_pReal allocate (defgradold ( resolution(1),resolution(2),resolution(3),3,3)); defgradold = 0.0_pReal
allocate (coordinates(3,resolution(1),resolution(2),resolution(3))); coordinates = 0.0_pReal allocate (coordinates(3,resolution(1),resolution(2),resolution(3))); coordinates = 0.0_pReal
@ -399,6 +339,67 @@ program DAMASK_spectral
! Initialization of CPFEM_general (= constitutive law) and of deformation gradient field ! Initialization of CPFEM_general (= constitutive law) and of deformation gradient field
call CPFEM_initAll(bc_temperature(1),1_pInt,1_pInt) call CPFEM_initAll(bc_temperature(1),1_pInt,1_pInt)
!Output of geom file
!$OMP CRITICAL (write2out)
print '(a)', '******************************************************'
print '(a,a)', 'Geom File Name: ',trim(path)//'.geom'
print '(a)', '------------------------------------------------------'
print '(a,/,i12,i12,i12)','resolution a b c:', resolution
print '(a,/,f12.5,f12.5,f12.5)','dimension x y z:', geomdimension
print '(a,i5)','homogenization: ',homog
print '(a,L)','spectralPictureMode: ',spectralPictureMode
print '(a)', '******************************************************'
print '(a,a)','Loadcase File Name: ',trim(getLoadcaseName())
if (bc_followFormerTrajectory(1)) then
call IO_warning(33) ! cannot guess along trajectory for first step of first loadcase
bc_followFormerTrajectory(1) = .false.
endif
! consistency checks and output of loadcase
do loadcase = 1, N_Loadcases
print '(a)', '------------------------------------------------------'
print '(a,i5)', 'Loadcase: ', loadcase
if (.not. bc_followFormerTrajectory(loadcase)) &
print '(a)', 'Drop Guessing Along Trajectory'
if (any(bc_mask(:,:,1,loadcase) .eqv. bc_mask(1:3,1:3,2,loadcase)))& ! exclusive or masking only
call IO_error(31,loadcase)
if (any(bc_mask(1:3,1:3,2,loadcase).and.transpose(bc_mask(1:3,1:3,2,loadcase)).and.& !checking if no rotation is allowed by stress BC
reshape((/.false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false./),(/3,3/))))&
call IO_error(38,loadcase)
if (bc_velGradApplied(loadcase)) then
do j = 1, 3
if (any(bc_mask(j,1:3,1,loadcase) .eqv. .true.) .and.&
any(bc_mask(j,1:3,1,loadcase) .eqv. .false.)) call IO_error(32,loadcase) ! each line should be either fully or not at all defined
enddo
print '(a,/,3(3(f12.6,x)/))','Velocity Gradient:', merge(math_transpose3x3(bc_deformation(1:3,1:3,loadcase)),&
reshape(spread(DAMASK_NaN,1,9),(/3,3/)),&
transpose(bc_mask(1:3,1:3,1,loadcase)))
else
print '(a,/,3(3(f12.6,x)/))','Change of Deformation Gradient:', merge(math_transpose3x3(bc_deformation(1:3,1:3,loadcase)),&
reshape(spread(DAMASK_NaN,1,9),(/3,3/)),&
transpose(bc_mask(1:3,1:3,1,loadcase)))
endif
print '(a,/,3(3(f12.6,x)/))','Stress Boundary Condition/MPa:',merge(math_transpose3x3(bc_stress(1:3,1:3,loadcase)),&
reshape(spread(DAMASK_NaN,1,9),(/3,3/)),&
transpose(bc_mask(:,:,2,loadcase)))*1e-6
if (any(abs(math_mul33x33(bc_rotation(1:3,1:3,loadcase),math_transpose3x3(bc_rotation(1:3,1:3,loadcase)))-math_I3)&
>reshape(spread(tol_rotation,1,9),(/3,3/)))&
.or. abs(math_det3x3(bc_rotation(1:3,1:3,loadcase)))>1.0_pReal + tol_rotation) call IO_error(46,loadcase)
if (any(bc_rotation(1:3,1:3,loadcase)/=math_I3)) then
print '(a,/,3(3(f12.6,x)/))','Rotation of BCs:',math_transpose3x3(bc_rotation(1:3,1:3,loadcase))
print '(a,f12.6)', 'Determinant of Rotation matrix', math_det3x3(bc_rotation(1:3,1:3,loadcase))
endif
if (bc_timeIncrement(loadcase) < 0.0_pReal) call IO_error(34,loadcase) ! negative time increment
print '(a,f12.6)','Temperature: ',bc_temperature(loadcase)
print '(a,f12.6)','Time: ',bc_timeIncrement(loadcase)
if (bc_steps(loadcase) < 1_pInt) call IO_error(35,loadcase) ! non-positive increment count
print '(a,i5)','Increments: ',bc_steps(loadcase)
if (bc_frequency(loadcase) < 1_pInt) call IO_error(36,loadcase) ! non-positive result frequency
print '(a,i5)','Freq. of Output: ',bc_frequency(loadcase)
enddo
print '(a)', '******************************************************'
!$OMP END CRITICAL (write2out)
ielem = 0_pInt ielem = 0_pInt
c_current = 0.0_pReal c_current = 0.0_pReal
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)
@ -410,8 +411,8 @@ program DAMASK_spectral
c_current = c_current + dPdF c_current = c_current + dPdF
enddo; enddo; enddo enddo; enddo; enddo
c0_reference = c_current * wgt ! linear reference material stiffness c0_reference = c_current * wgt ! linear reference material stiffness
c_prev = c0_reference c_prev = math_rotate_forward3x3x3x3(c0_reference,bc_rotation(1:3,1:3,loadcase))
! rotate_forward: lab -> load system
if (debug_verbosity > 1) then if (debug_verbosity > 1) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write (6,*) 'First Call to CPFEM_general finished' write (6,*) 'First Call to CPFEM_general finished'
@ -462,14 +463,18 @@ program DAMASK_spectral
allocate (workfft(resolution(1)+2,resolution(2),resolution(3),3,3)); workfft = 0.0_pReal allocate (workfft(resolution(1)+2,resolution(2),resolution(3),3,3)); workfft = 0.0_pReal
! Initialization of fftw (see manual on fftw.org for more details) ! Initialization of fftw (see manual on fftw.org for more details)
#ifdef _OPENMP
if(DAMASK_NumThreadsInt>0_pInt) then if(DAMASK_NumThreadsInt>0_pInt) then
call dfftw_init_threads(ierr) call dfftw_init_threads(ierr)
if(ierr == 0_pInt) call IO_error(104,ierr) if(ierr == 0_pInt) call IO_error(104,ierr)
call dfftw_plan_with_nthreads(DAMASK_NumThreadsInt) call dfftw_plan_with_nthreads(DAMASK_NumThreadsInt)
endif endif
#endif
!is not working, have to find out how it is working in FORTRAN
!call dfftw_timelimit(fftw_timelimit) !call dfftw_timelimit(fftw_timelimit)
! basically a translation from fftw.h
select case(IO_lc(fftw_planner_flag)) select case(IO_lc(fftw_planner_flag))
case('estimate','fftw_estimate') case('estimate','fftw_estimate')
fftw_flag = 64 fftw_flag = 64
@ -479,6 +484,11 @@ program DAMASK_spectral
fftw_flag= 32 fftw_flag= 32
case('exhaustive','fftw_exhaustive') case('exhaustive','fftw_exhaustive')
fftw_flag = 8 fftw_flag = 8
case default
!$OMP CRITICAL (write2out)
write (6,*) 'No valid parameter for FFTW given, using FFTW_PATIENT'
!$OMP END CRITICAL (write2out)
fftw_flag = 32
end select end select
call dfftw_plan_many_dft_r2c(fftw_plan(1),3,(/resolution(1),resolution(2),resolution(3)/),9,& call dfftw_plan_many_dft_r2c(fftw_plan(1),3,(/resolution(1),resolution(2),resolution(3)/),9,&
@ -505,15 +515,12 @@ program DAMASK_spectral
write(538), 'logscale', bc_logscale ! one entry per loadcase (0: linear, 1: log) write(538), 'logscale', bc_logscale ! one entry per loadcase (0: linear, 1: log)
write(538), 'frequencies', bc_frequency ! one entry per loadcase write(538), 'frequencies', bc_frequency ! one entry per loadcase
write(538), 'times', bc_timeIncrement ! one entry per loadcase write(538), 'times', bc_timeIncrement ! one entry per loadcase
write(538), 'increments' ! one entry per loadcase bc_timeIncrement(1)= bc_timeIncrement(1) + 1_pInt
incPosition = ftellI8(538) write(538), 'increments', bc_timeIncrement ! one entry per loadcase
write(538), 0_pInt ! end of header bc_timeIncrement(1)= bc_timeIncrement(1) - 1_pInt
write(538), 'startingIncrement', writtenOutCounter
write(538), 'eoh' ! end of header write(538), 'eoh' ! end of header
write(538), materialpoint_results(:,1,:) ! initial (non-deformed) results write(538), materialpoint_results(:,1,:) ! initial (non-deformed) results
ierr = fseek(538,incPosition,0)
if (ierr/=0_pInt ) call IO_error(107,ext_msg = 'inc 1')
writtenOutCounter = 1_pInt
write(538), writtenOutCounter ! 0 inc is written out
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
! Initialization done ! Initialization done
@ -558,7 +565,7 @@ program DAMASK_spectral
if (bc_velGradApplied(loadcase)) & ! calculate fDot from given L and current F if (bc_velGradApplied(loadcase)) & ! calculate fDot from given L and current F
fDot = math_mul33x33(bc_deformation(1:3,1:3,loadcase), defgradAim) fDot = math_mul33x33(bc_deformation(1:3,1:3,loadcase), defgradAim)
!winding forward of deformation aim !winding forward of deformation aim in loadcase system
temp33_Real = defgradAim temp33_Real = defgradAim
defgradAim = defgradAim & defgradAim = defgradAim &
+ guessmode * mask_stress * (defgradAim - defgradAimOld) & + guessmode * mask_stress * (defgradAim - defgradAimOld) &
@ -567,12 +574,14 @@ program DAMASK_spectral
! update local deformation gradient ! update local deformation gradient
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)
temp33_Real = defgrad(i,j,k,:,:) temp33_Real = defgrad(i,j,k,1:3,1:3)
if (bc_velGradApplied(loadcase)) & ! use velocity gradient to calculate new deformation gradient (if not guessing) if (bc_velGradApplied(loadcase)) & ! use velocity gradient to calculate new deformation gradient (if not guessing)
fDot = math_mul33x33(bc_deformation(1:3,1:3,loadcase),defgradold(i,j,k,1:3,1:3)) fDot = math_mul33x33(bc_deformation(1:3,1:3,loadcase),&
math_rotate_forward3x3(defgradold(i,j,k,1:3,1:3),bc_rotation(1:3,1:3,loadcase)))
defgrad(i,j,k,1:3,1:3) = defgrad(i,j,k,1:3,1:3) & ! decide if guessing along former trajectory or apply homogeneous addon defgrad(i,j,k,1:3,1:3) = defgrad(i,j,k,1:3,1:3) & ! decide if guessing along former trajectory or apply homogeneous addon
+ guessmode * (defgrad(i,j,k,1:3,1:3) - defgradold(i,j,k,1:3,1:3))& ! guessing... + guessmode * (defgrad(i,j,k,1:3,1:3) - defgradold(i,j,k,1:3,1:3))& ! guessing...
+ (1.0_pReal-guessmode) * mask_defgrad * fDot *timeinc ! apply the prescribed value where deformation is given if not guessing + math_rotate_backward3x3((1.0_pReal-guessmode) * mask_defgrad * fDot,&
bc_rotation(1:3,1:3,loadcase)) *timeinc ! apply the prescribed value where deformation is given if not guessing
defgradold(i,j,k,1:3,1:3) = temp33_Real defgradold(i,j,k,1:3,1:3) = temp33_Real
enddo; enddo; enddo enddo; enddo; enddo
guessmode = 1.0_pReal ! keep guessing along former trajectory during same loadcase guessmode = 1.0_pReal ! keep guessing along former trajectory during same loadcase
@ -654,14 +663,16 @@ program DAMASK_spectral
print '(a,/,3(3(f12.7,x)/))', 'Piola-Kirchhoff Stress / MPa: ',math_transpose3x3(pstress_av)/1.e6 print '(a,/,3(3(f12.7,x)/))', 'Piola-Kirchhoff Stress / MPa: ',math_transpose3x3(pstress_av)/1.e6
err_stress_tol = 0.0_pReal err_stress_tol = 0.0_pReal
pstress_av_load = math_rotate_forward3x3(pstress_av,bc_rotation(1:3,1:3,loadcase))
if(size_reduced > 0_pInt) then ! calculate stress BC if applied if(size_reduced > 0_pInt) then ! calculate stress BC if applied
err_stress = maxval(abs(mask_stress * (pstress_av - bc_stress(1:3,1:3,loadcase)))) ! maximum deviaton (tensor norm not applicable) err_stress = maxval(abs(mask_stress * (pstress_av_load - bc_stress(1:3,1:3,loadcase)))) ! maximum deviaton (tensor norm not applicable)
err_stress_tol = maxval(abs(mask_defgrad * pstress_av)) * err_stress_tolrel ! don't use any tensor norm because the comparison should be coherent err_stress_tol = maxval(abs(mask_defgrad * pstress_av_load)) * err_stress_tolrel ! don't use any tensor norm because the comparison should be coherent
print '(A,/)', '== Correcting Deformation Gradient to Fullfill BCs =========' print '(A,/)', '== Correcting Deformation Gradient to Fullfill BCs ========='
print '(2(a,E10.5)/)', 'Error Stress = ',err_stress, ', Tol. = ', err_stress_tol print '(2(a,E10.5)/)', 'Error Stress = ',err_stress, ', Tol. = ', err_stress_tol
defgradAimCorr = - math_mul3333xx33(s_prev, ((pstress_av - bc_stress(1:3,1:3,loadcase)))) ! residual on given stress components defgradAimCorr = - math_mul3333xx33(s_prev, ((pstress_av_load - bc_stress(1:3,1:3,loadcase)))) ! residual on given stress components
defgradAim = defgradAim + defgradAimCorr defgradAim = defgradAim + defgradAimCorr
print '(a,/,3(3(f12.7,x)/))', 'Deformation Aim: ',math_transpose3x3(defgradAim) print '(a,/,3(3(f12.7,x)/))', 'Deformation Aim: ',math_transpose3x3(math_rotate_backward3x3(&
defgradAim,bc_rotation(1:3,1:3,loadcase)))
print '(a,x,f12.7,/)' , 'Determinant of Deformation Aim: ', math_det3x3(defgradAim) print '(a,x,f12.7,/)' , 'Determinant of Deformation Aim: ', math_det3x3(defgradAim)
endif endif
print '(A,/)', '== Calculating Equilibrium Using Spectral Method ===========' print '(A,/)', '== Calculating Equilibrium Using Spectral Method ==========='
@ -718,24 +729,22 @@ program DAMASK_spectral
defgrad = defgrad + workfft(1:resolution(1),:,:,:,:)*wgt defgrad = defgrad + workfft(1:resolution(1),:,:,:,:)*wgt
do m = 1,3; do n = 1,3 do m = 1,3; do n = 1,3
defgrad_av(m,n) = sum(defgrad(:,:,:,m,n))*wgt defgrad_av(m,n) = sum(defgrad(:,:,:,m,n))*wgt
defgrad(:,:,:,m,n) = defgrad(:,:,:,m,n) + (defgradAim(m,n) - defgrad_av(m,n)) ! anticipated target minus current state enddo; enddo
defgradAim_lab = math_rotate_backward3x3(defgradAim,bc_rotation(1:3,1:3,loadcase))
do m = 1,3; do n = 1,3
defgrad(:,:,:,m,n) = defgrad(:,:,:,m,n) + (defgradAim_lab(m,n) - defgrad_av(m,n)) ! anticipated target minus current state
enddo; enddo enddo; enddo
print '(2(a,E10.5)/)', 'Error Divergence = ',err_div, ', Tol. = ', err_div_tol print '(2(a,E10.5)/)', 'Error Divergence = ',err_div, ', Tol. = ', err_div_tol
enddo ! end looping when convergency is achieved enddo ! end looping when convergency is achieved
c_prev = c_current*wgt ! calculate stiffness for next step c_prev = math_rotate_forward3x3x3x3(c_current*wgt,bc_rotation(1:3,1:3,loadcase)) ! calculate stiffness for next step
!ToDo: Incfluence for next loadcase
! if (mod(step,bc_frequency(loadcase)) == 0_pInt) then ! at output frequency if (mod(step,bc_frequency(loadcase)) == 0_pInt) then ! at output frequency
! ierr = fseek(538,0,2) write(538), materialpoint_results(:,1,:) ! write result to file
! if (ierr/=0_pInt ) call IO_error(107,ext_msg = 'searching backwards') writtenOutCounter = writtenOutCounter + 1_pInt
! write(538), materialpoint_results(:,1,:) ! write result to file endif
! writtenOutCounter = writtenOutCounter + 1_pInt
! ierr = fseek(538,incPosition,0)
! if (ierr/=0_pInt ) call IO_error(107,ext_msg = 'searching Position of Inc')
! write(538), writtenOutCounter ! write result to file
! endif
if(err_div<=err_div_tol .and. err_stress<=err_stress_tol) then if(err_div<=err_div_tol .and. err_stress<=err_stress_tol) then
print '(2(A,I5.5),A,/)', '== Step = ',step, ' of Loadcase = ',loadcase, ' Converged ==============' print '(2(A,I5.5),A,/)', '== Step = ',step, ' of Loadcase = ',loadcase, ' Converged =============='
else else

2
code/IO.f90
View File

@ -1162,6 +1162,8 @@ endfunction
msg = 'dimension in spectral mesh' msg = 'dimension in spectral mesh'
case (45) case (45)
msg = 'incomplete information in spectral mesh header' msg = 'incomplete information in spectral mesh header'
case (46)
msg = 'not a rotation defined for loadcase rotation'
case (50) case (50)
msg = 'writing constitutive output description' msg = 'writing constitutive output description'
case (100) case (100)

48
code/math.f90
View File

@ -2966,5 +2966,53 @@ end subroutine
endfunction math_volTetrahedron endfunction math_volTetrahedron
!**************************************************************************
! rotate 3x3 tensor forward
!**************************************************************************
pure function math_rotate_forward3x3(tensor,rot_tensor)
use prec, only: pReal
implicit none
real(pReal), dimension(3,3) :: math_rotate_forward3x3
real(pReal), dimension(3,3), intent(in) :: tensor, rot_tensor
math_rotate_forward3x3 = math_mul33x33(rot_tensor,&
math_mul33x33(tensor,math_transpose3x3(rot_tensor)))
endfunction math_rotate_forward3x3
!**************************************************************************
! rotate 3x3 tensor backward
!**************************************************************************
pure function math_rotate_backward3x3(tensor,rot_tensor)
use prec, only: pReal
implicit none
real(pReal), dimension(3,3) :: math_rotate_backward3x3
real(pReal), dimension(3,3), intent(in) :: tensor, rot_tensor
math_rotate_backward3x3 = math_mul33x33(math_transpose3x3(rot_tensor),&
math_mul33x33(tensor,rot_tensor))
endfunction math_rotate_backward3x3
!**************************************************************************
! rotate 3x3x3x3 tensor forward
! DUMMY FUNCTION, does nothing!
!**************************************************************************
pure function math_rotate_forward3x3x3x3(tensor,rot_tensor)
use prec, only: pReal
implicit none
real(pReal), dimension(3,3,3,3) :: math_rotate_forward3x3x3x3
real(pReal), dimension(3,3), intent(in) :: rot_tensor
real(pReal), dimension(3,3,3,3), intent(in) :: tensor
math_rotate_forward3x3x3x3 = tensor
endfunction math_rotate_forward3x3x3x3
END MODULE math END MODULE math

14
code/numerics.f90
View File

@ -26,7 +26,7 @@ use prec, only: pInt, pReal
implicit none implicit none
character(len=64), parameter :: numerics_configFile = 'numerics.config' ! name of configuration file character(len=64), parameter :: numerics_configFile = 'numerics.config' ! name of configuration file
integer(pInt) iJacoStiffness, & ! frequency of stiffness update integer(pInt) :: iJacoStiffness, & ! frequency of stiffness update
iJacoLpresiduum, & ! frequency of Jacobian update of residuum in Lp iJacoLpresiduum, & ! frequency of Jacobian update of residuum in Lp
nHomog, & ! homogenization loop limit (only for debugging info, loop limit is determined by "subStepMinHomog") nHomog, & ! homogenization loop limit (only for debugging info, loop limit is determined by "subStepMinHomog")
nMPstate, & ! materialpoint state loop limit nMPstate, & ! materialpoint state loop limit
@ -36,7 +36,7 @@ integer(pInt) iJacoStiffness, & ! freque
pert_method, & ! method used in perturbation technique for tangent pert_method, & ! method used in perturbation technique for tangent
numerics_integrationMode ! integration mode 1 = central solution ; integration mode 2 = perturbation numerics_integrationMode ! integration mode 1 = central solution ; integration mode 2 = perturbation
integer(pInt), dimension(2) :: numerics_integrator ! method used for state integration (central & perturbed state) integer(pInt), dimension(2) :: numerics_integrator ! method used for state integration (central & perturbed state)
real(pReal) relevantStrain, & ! strain increment considered significant (used by crystallite to determine whether strain inc is considered significant) real(pReal) :: relevantStrain, & ! strain increment considered significant (used by crystallite to determine whether strain inc is considered significant)
defgradTolerance, & ! deviation of deformation gradient that is still allowed (used by CPFEM to determine outdated ffn1) defgradTolerance, & ! deviation of deformation gradient that is still allowed (used by CPFEM to determine outdated ffn1)
pert_Fg, & ! strain perturbation for FEM Jacobi pert_Fg, & ! strain perturbation for FEM Jacobi
subStepMinCryst, & ! minimum (relative) size of sub-step allowed during cutback in crystallite subStepMinCryst, & ! minimum (relative) size of sub-step allowed during cutback in crystallite
@ -68,7 +68,8 @@ real(pReal) relevantStrain, & ! strain
err_div_tol, & ! error of divergence in fourier space err_div_tol, & ! error of divergence in fourier space
err_stress_tol, & ! absolut stress error, will be computed from err_stress_tolrel (dont prescribe a value) err_stress_tol, & ! absolut stress error, will be computed from err_stress_tolrel (dont prescribe a value)
err_stress_tolrel, & ! factor to multiply with highest stress to get err_stress_tol err_stress_tolrel, & ! factor to multiply with highest stress to get err_stress_tol
fftw_timelimit ! sets the timelimit of plan creation for FFTW, see manual on www.fftw.org fftw_timelimit, & ! sets the timelimit of plan creation for FFTW, see manual on www.fftw.org
tol_rotation ! tolerance of rotation specified in loadcase
character(len=64) fftw_planner_flag ! sets the planig-rigor flag, see manual on www.fftw.org character(len=64) fftw_planner_flag ! sets the planig-rigor flag, see manual on www.fftw.org
logical memory_efficient ! for fast execution (pre calculation of gamma_hat) logical memory_efficient ! for fast execution (pre calculation of gamma_hat)
integer(pInt) itmax , & ! maximum number of iterations integer(pInt) itmax , & ! maximum number of iterations
@ -169,8 +170,8 @@ subroutine numerics_init()
itmax = 20_pInt ! Maximum iteration number itmax = 20_pInt ! Maximum iteration number
memory_efficient = .true. ! Precalculate Gamma-operator (81 double per point) memory_efficient = .true. ! Precalculate Gamma-operator (81 double per point)
fftw_timelimit = -1.0_pReal ! no timelimit of plan creation for FFTW fftw_timelimit = -1.0_pReal ! no timelimit of plan creation for FFTW
fftw_planner_flag ='patient' fftw_planner_flag ='FFTW_PATIENT'
tol_rotation = 1.0e-4
!* Random seeding parameters: added <<<updated 27.08.2009>>> !* Random seeding parameters: added <<<updated 27.08.2009>>>
fixedSeed = 0_pInt fixedSeed = 0_pInt
@ -287,6 +288,8 @@ subroutine numerics_init()
fftw_timelimit = IO_floatValue(line,positions,2) fftw_timelimit = IO_floatValue(line,positions,2)
case ('fftw_planner_flag') case ('fftw_planner_flag')
fftw_planner_flag = IO_stringValue(line,positions,2) fftw_planner_flag = IO_stringValue(line,positions,2)
case ('tol_rotation')
tol_rotation = IO_floatValue(line,positions,2)
!* Random seeding parameters !* Random seeding parameters
case ('fixed_seed') case ('fixed_seed')
@ -359,6 +362,7 @@ subroutine numerics_init()
write(6,'(a24,x,e8.1)') 'fftw_timelimit: ',fftw_timelimit write(6,'(a24,x,e8.1)') 'fftw_timelimit: ',fftw_timelimit
endif endif
write(6,'(a24,x,a)') 'fftw_planner_flag: ',trim(fftw_planner_flag) write(6,'(a24,x,a)') 'fftw_planner_flag: ',trim(fftw_planner_flag)
write(6,'(a24,x,e8.1)') 'tol_rotation: ',tol_rotation
write(6,*) write(6,*)
!* Random seeding parameters !* Random seeding parameters