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

simplified the rate calculation interface in DAMASK_spectral_utilities.f90, corrected bug in rotation parsing and added debug option "rotation" to spectral solver options. This will show the current average stress and deformation additionally in lab coordinate system

This commit is contained in:
Martin Diehl 2012-11-28 15:04:05 +00:00
parent 85b0861893
commit 5ea0139678
6 changed files with 43 additions and 33 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
code/DAMASK_spectral_driver.f90
View File

@ -196,18 +196,17 @@ program DAMASK_spectral_Driver
loadCases(currentLoadCase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory loadCases(currentLoadCase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory
case('euler') ! rotation of currentLoadCase given in euler angles case('euler') ! rotation of currentLoadCase given in euler angles
temp_valueVector = 0.0_pReal temp_valueVector = 0.0_pReal
l = 0_pInt ! assuming values given in degrees l = 1_pInt ! assuming values given in degrees
k = 0_pInt ! assuming keyword indicating degree/radians k = 1_pInt ! assuming keyword indicating degree/radians present
select case (IO_lc(IO_stringValue(line,positions,i+1_pInt))) select case (IO_lc(IO_stringValue(line,positions,i+1_pInt)))
case('rad','radian')
l = 1_pInt
k = 1_pInt
case('deg','degree') case('deg','degree')
l = 1_pInt ! for conversion from degree to radian case('rad','radian') ! don't convert from degree to radian
l = 0_pInt
case default case default
k = 0_pInt
end select end select
forall(j = 1_pInt:3_pInt) temp_valueVector(j) = IO_floatValue(line,positions,i+k+j) forall(j = 1_pInt:3_pInt) temp_valueVector(j) = IO_floatValue(line,positions,i+k+j)
if (k == 1_pInt) temp_valueVector(1:3) = temp_valueVector(1:3) * inRad ! convert to rad if (l == 1_pInt) temp_valueVector(1:3) = temp_valueVector(1:3) * inRad ! convert to rad
loadCases(currentLoadCase)%rotation = math_EulerToR(temp_valueVector(1:3)) ! convert rad Eulers to rotation matrix loadCases(currentLoadCase)%rotation = math_EulerToR(temp_valueVector(1:3)) ! convert rad Eulers to rotation matrix
case('rotation','rot') ! assign values for the rotation of currentLoadCase matrix case('rotation','rot') ! assign values for the rotation of currentLoadCase matrix
temp_valueVector = 0.0_pReal temp_valueVector = 0.0_pReal

4
code/DAMASK_spectral_solverAL.f90
View File

@ -291,9 +291,9 @@ type(tSolutionState) function &
mesh_ipCoordinates = 0.0_pReal !reshape(mesh_deformedCoordsFFT(geomdim,& mesh_ipCoordinates = 0.0_pReal !reshape(mesh_deformedCoordsFFT(geomdim,&
!reshape(F,[3,3,res(1),res(2),res(3)])),[3,1,mesh_NcpElems]) !reshape(F,[3,3,res(1),res(2),res(3)])),[3,1,mesh_NcpElems])
Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
timeinc,timeinc_old,guess,F_lastInc,reshape(F,[3,3,res(1),res(2),res(3)])) timeinc_old,guess,F_lastInc,reshape(F,[3,3,res(1),res(2),res(3)]))
F_lambdaDot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & F_lambdaDot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
timeinc,timeinc_old,guess,F_lambda_lastInc,reshape(F_lambda,[3,3,res(1),res(2),res(3)])) timeinc_old,guess,F_lambda_lastInc,reshape(F_lambda,[3,3,res(1),res(2),res(3)]))
F_lastInc = reshape(F,[3,3,res(1),res(2),res(3)]) F_lastInc = reshape(F,[3,3,res(1),res(2),res(3)])
F_lambda_lastInc = reshape(F_lambda,[3,3,res(1),res(2),res(3)]) F_lambda_lastInc = reshape(F_lambda,[3,3,res(1),res(2),res(3)])

21
code/DAMASK_spectral_solverBasic.f90
View File

@ -171,7 +171,8 @@ type(tSolutionState) function &
Utilities_FFTbackward, & Utilities_FFTbackward, &
Utilities_updateGamma, & Utilities_updateGamma, &
Utilities_constitutiveResponse, & Utilities_constitutiveResponse, &
Utilities_calculateRate Utilities_calculateRate, &
debugRotation
use FEsolving, only: & use FEsolving, only: &
restartWrite, & restartWrite, &
restartRead, & restartRead, &
@ -199,8 +200,7 @@ type(tSolutionState) function &
real(pReal), dimension(3,3,3,3) :: & real(pReal), dimension(3,3,3,3) :: &
S !< current average compliance S !< current average compliance
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
F_aim_lab, & F_aim_lastIter, & !< aim of last iteration
F_aim_lab_lastIter, & !< aim of last iteration
P_av P_av
real(pReal), dimension(3,3,res(1),res(2),res(3)) :: P real(pReal), dimension(3,3,res(1),res(2),res(3)) :: P
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -258,12 +258,12 @@ type(tSolutionState) function &
endif endif
if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old
F_aim_lastInc = F_aim F_aim_lastInc = F_aim
Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & Fdot = utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
timeinc,timeinc_old,guess,F_lastInc,F) timeinc_old,guess,F_lastInc,F)
F_lastInc = F F_lastInc = F
endif endif
F_aim = F_aim + f_aimDot * timeinc F_aim = F_aim + f_aimDot * timeinc
F = Utilities_forwardField(timeinc,F_aim,F_lastInc,Fdot) !I think F aim should be rotated here F = Utilities_forwardField(timeinc,math_rotate_backward33(F_aim,rotation_BC),F_lastInc,Fdot) !I think F aim should be rotated here
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator) ! update stiffness (and gamma operator)
@ -280,12 +280,14 @@ type(tSolutionState) function &
! report begin of new iteration ! report begin of new iteration
write(6,'(/,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & write(6,'(/,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), &
' @ Iter. ', itmin, '≤',iter, '≤', itmax ' @ Iter. ', itmin, '≤',iter, '≤', itmax
if (debugRotation) &
write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim (lab)=', &
math_transpose33(math_rotate_backward33(F_aim,rotation_BC))
write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim =', & write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim =', &
math_transpose33(F_aim) math_transpose33(F_aim)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate constitutive response ! evaluate constitutive response
F_aim_lab_lastIter = math_rotate_backward33(F_aim,rotation_BC) F_aim_lastIter = F_aim
basic_solution%termIll = .false. basic_solution%termIll = .false.
call Utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,& call Utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,&
P,C,P_av,ForwardData,rotation_BC) P,C,P_av,ForwardData,rotation_BC)
@ -297,7 +299,6 @@ type(tSolutionState) function &
! stress BC handling ! stress BC handling
F_aim = F_aim - math_mul3333xx33(S, ((P_av - P_BC%values))) ! S = 0.0 for no bc F_aim = F_aim - math_mul3333xx33(S, ((P_av - P_BC%values))) ! S = 0.0 for no bc
err_stress = maxval(abs(P_BC%maskFloat * (P_av - P_BC%values))) ! mask = 0.0 for no bc err_stress = maxval(abs(P_BC%maskFloat * (P_av - P_BC%values))) ! mask = 0.0 for no bc
F_aim_lab = math_rotate_backward33(F_aim,rotation_BC) ! boundary conditions from load frame into lab (Fourier) frame
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! updated deformation gradient using fix point algorithm of basic scheme ! updated deformation gradient using fix point algorithm of basic scheme
@ -306,7 +307,7 @@ type(tSolutionState) function &
order=[4,5,1,2,3]) ! field real has a different order order=[4,5,1,2,3]) ! field real has a different order
call Utilities_FFTforward() call Utilities_FFTforward()
err_div = Utilities_divergenceRMS() err_div = Utilities_divergenceRMS()
call Utilities_fourierConvolution(F_aim_lab_lastIter - F_aim_lab) call Utilities_fourierConvolution(math_rotate_backward33(F_aim_lastIter-F_aim,rotation_BC))
call Utilities_FFTbackward() call Utilities_FFTbackward()
F = F - reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),shape(F),order=[3,4,5,1,2]) ! F(x)^(n+1) = F(x)^(n) + correction; *wgt: correcting for missing normalization F = F - reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),shape(F),order=[3,4,5,1,2]) ! F(x)^(n+1) = F(x)^(n) + correction; *wgt: correcting for missing normalization
basic_solution%converged = basic_Converged(err_div,P_av,err_stress,P_av) basic_solution%converged = basic_Converged(err_div,P_av,err_stress,P_av)

2
code/DAMASK_spectral_solverBasicPETSc.f90
View File

@ -270,7 +270,7 @@ type(tSolutionState) function &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update coordinates and rate and forward last inc ! update coordinates and rate and forward last inc
Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
timeinc,timeinc_old,guess,F_lastInc,reshape(F,[3,3,res(1),res(2),res(3)])) timeinc_old,guess,F_lastInc,reshape(F,[3,3,res(1),res(2),res(3)]))
F_lastInc = reshape(F,[3,3,res(1),res(2),res(3)]) F_lastInc = reshape(F,[3,3,res(1),res(2),res(3)])
endif endif
F_aim = F_aim + f_aimDot * timeinc F_aim = F_aim + f_aimDot * timeinc

20
code/DAMASK_spectral_utilities.f90
View File

@ -46,7 +46,8 @@ module DAMASK_spectral_utilities
debugGeneral, & !< general debugging of spectral solver debugGeneral, & !< general debugging of spectral solver
debugDivergence, & !< debugging of divergence calculation (comparison to function used for post processing) debugDivergence, & !< debugging of divergence calculation (comparison to function used for post processing)
debugRestart, & !< debbuging of restart features debugRestart, & !< debbuging of restart features
debugFFTW !< doing additional FFT on scalar field and compare to results of strided 3D FFT debugFFTW, & !< doing additional FFT on scalar field and compare to results of strided 3D FFT
debugRotation !< also printing out results in lab frame
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! derived types ! derived types
@ -105,7 +106,8 @@ subroutine utilities_init()
debug_levelBasic, & debug_levelBasic, &
debug_spectralDivergence, & debug_spectralDivergence, &
debug_spectralRestart, & debug_spectralRestart, &
debug_spectralFFTW debug_spectralFFTW, &
debug_spectralRotation
use mesh, only: & use mesh, only: &
res, & res, &
res1_red, & res1_red, &
@ -132,6 +134,7 @@ subroutine utilities_init()
debugDivergence = iand(debug_level(debug_spectral),debug_spectralDivergence) /= 0 debugDivergence = iand(debug_level(debug_spectral),debug_spectralDivergence) /= 0
debugRestart = iand(debug_level(debug_spectral),debug_spectralRestart) /= 0 debugRestart = iand(debug_level(debug_spectral),debug_spectralRestart) /= 0
debugFFTW = iand(debug_level(debug_spectral),debug_spectralFFTW) /= 0 debugFFTW = iand(debug_level(debug_spectral),debug_spectralFFTW) /= 0
debugRotation = iand(debug_level(debug_spectral),debug_spectralRotation) /= 0
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocation ! allocation
@ -721,10 +724,14 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,&
C = C * wgt C = C * wgt
call debug_info() call debug_info()
P_av = math_rotate_forward33(sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt,rotation_BC) ! average of P rotated
restartWrite = .false. ! reset restartWrite status restartWrite = .false. ! reset restartWrite status
cutBack = .false. ! reset cutBack status cutBack = .false. ! reset cutBack status
P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P
if (debugRotation) &
write(6,'(a,/,3(3(2x,f12.7,1x)/))',advance='no') 'Piola-Kirchhoff stress (lab) / MPa =',&
math_transpose33(P_av)/1.e6_pReal
P_av = math_rotate_forward33(P_av,rotation_BC)
write(6,'(a,/,3(3(2x,f12.7,1x)/))',advance='no') 'Piola-Kirchhoff stress / MPa =',& write(6,'(a,/,3(3(2x,f12.7,1x)/))',advance='no') 'Piola-Kirchhoff stress / MPa =',&
math_transpose33(P_av)/1.e6_pReal math_transpose33(P_av)/1.e6_pReal
end subroutine utilities_constitutiveResponse end subroutine utilities_constitutiveResponse
@ -733,14 +740,13 @@ end subroutine utilities_constitutiveResponse
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates forward rate, either guessing or just add delta/timeinc !> @brief calculates forward rate, either guessing or just add delta/timeinc
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function utilities_calculateRate(delta_aim,timeinc,timeinc_old,guess,field_lastInc,field) pure function utilities_calculateRate(avRate,timeinc_old,guess,field_lastInc,field)
use mesh, only: & use mesh, only: &
res res
implicit none implicit none
real(pReal), intent(in), dimension(3,3) :: delta_aim !< homogeneous addon real(pReal), intent(in), dimension(3,3) :: avRate !< homogeneous addon
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
timeinc, & !< timeinc of current step
timeinc_old !< timeinc of last step timeinc_old !< timeinc of last step
logical, intent(in) :: & logical, intent(in) :: &
guess !< guess along former trajectory guess !< guess along former trajectory
@ -752,7 +758,7 @@ pure function utilities_calculateRate(delta_aim,timeinc,timeinc_old,guess,field_
if(guess) then if(guess) then
utilities_calculateRate = (field-field_lastInc) / timeinc_old utilities_calculateRate = (field-field_lastInc) / timeinc_old
else else
utilities_calculateRate = spread(spread(spread(delta_aim,3,res(1)),4,res(2)),5,res(3))/timeinc utilities_calculateRate = spread(spread(spread(avRate,3,res(1)),4,res(2)),5,res(3))
endif endif
end function utilities_calculateRate end function utilities_calculateRate

6
code/debug.f90
View File

@ -43,7 +43,8 @@ module debug
integer(pInt), parameter, public :: & integer(pInt), parameter, public :: &
debug_spectralRestart = debug_maxGeneral*2_pInt**1_pInt, & debug_spectralRestart = debug_maxGeneral*2_pInt**1_pInt, &
debug_spectralFFTW = debug_maxGeneral*2_pInt**2_pInt, & debug_spectralFFTW = debug_maxGeneral*2_pInt**2_pInt, &
debug_spectralDivergence = debug_maxGeneral*2_pInt**3_pInt debug_spectralDivergence = debug_maxGeneral*2_pInt**3_pInt, &
debug_spectralRotation = debug_maxGeneral*2_pInt**4_pInt
integer(pInt), parameter, public :: & integer(pInt), parameter, public :: &
debug_debug = 1_pInt, & debug_debug = 1_pInt, &
@ -232,6 +233,8 @@ subroutine debug_init
debug_level(what) = ior(debug_level(what), debug_spectralFFTW) debug_level(what) = ior(debug_level(what), debug_spectralFFTW)
case('divergence') case('divergence')
debug_level(what) = ior(debug_level(what), debug_spectralDivergence) debug_level(what) = ior(debug_level(what), debug_spectralDivergence)
case('rotation')
debug_level(what) = ior(debug_level(what), debug_spectralRotation)
end select end select
enddo enddo
endif endif
@ -306,6 +309,7 @@ subroutine debug_init
if(iand(debug_level(i),debug_spectralRestart) /= 0) write(6,'(a)') ' restart' if(iand(debug_level(i),debug_spectralRestart) /= 0) write(6,'(a)') ' restart'
if(iand(debug_level(i),debug_spectralFFTW) /= 0) write(6,'(a)') ' FFTW' if(iand(debug_level(i),debug_spectralFFTW) /= 0) write(6,'(a)') ' FFTW'
if(iand(debug_level(i),debug_spectralDivergence)/= 0) write(6,'(a)') ' divergence' if(iand(debug_level(i),debug_spectralDivergence)/= 0) write(6,'(a)') ' divergence'
if(iand(debug_level(i),debug_spectralRotation) /= 0) write(6,'(a)') ' rotation'
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
enddo enddo