From 5ea01396782c87508d697c07d4bf526c4023da3c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 28 Nov 2012 15:04:05 +0000 Subject: [PATCH] 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 --- code/DAMASK_spectral_driver.f90 | 15 ++++++------ code/DAMASK_spectral_solverAL.f90 | 4 ++-- code/DAMASK_spectral_solverBasic.f90 | 29 ++++++++++++----------- code/DAMASK_spectral_solverBasicPETSc.f90 | 2 +- code/DAMASK_spectral_utilities.f90 | 20 ++++++++++------ code/debug.f90 | 6 ++++- 6 files changed, 43 insertions(+), 33 deletions(-) diff --git a/code/DAMASK_spectral_driver.f90 b/code/DAMASK_spectral_driver.f90 index 11f7879d6..d01de374b 100644 --- a/code/DAMASK_spectral_driver.f90 +++ b/code/DAMASK_spectral_driver.f90 @@ -196,18 +196,17 @@ program DAMASK_spectral_Driver loadCases(currentLoadCase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory case('euler') ! rotation of currentLoadCase given in euler angles temp_valueVector = 0.0_pReal - l = 0_pInt ! assuming values given in degrees - k = 0_pInt ! assuming keyword indicating degree/radians + l = 1_pInt ! assuming values given in degrees + k = 1_pInt ! assuming keyword indicating degree/radians present select case (IO_lc(IO_stringValue(line,positions,i+1_pInt))) - case('rad','radian') - l = 1_pInt - k = 1_pInt case('deg','degree') - l = 1_pInt ! for conversion from degree to radian - case default + case('rad','radian') ! don't convert from degree to radian + l = 0_pInt + case default + k = 0_pInt end select 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 case('rotation','rot') ! assign values for the rotation of currentLoadCase matrix temp_valueVector = 0.0_pReal diff --git a/code/DAMASK_spectral_solverAL.f90 b/code/DAMASK_spectral_solverAL.f90 index a7f0eae2a..768de4d98 100644 --- a/code/DAMASK_spectral_solverAL.f90 +++ b/code/DAMASK_spectral_solverAL.f90 @@ -291,9 +291,9 @@ type(tSolutionState) function & mesh_ipCoordinates = 0.0_pReal !reshape(mesh_deformedCoordsFFT(geomdim,& !reshape(F,[3,3,res(1),res(2),res(3)])),[3,1,mesh_NcpElems]) 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), & - 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_lambda_lastInc = reshape(F_lambda,[3,3,res(1),res(2),res(3)]) diff --git a/code/DAMASK_spectral_solverBasic.f90 b/code/DAMASK_spectral_solverBasic.f90 index f0330fbda..c23826a4e 100644 --- a/code/DAMASK_spectral_solverBasic.f90 +++ b/code/DAMASK_spectral_solverBasic.f90 @@ -171,7 +171,8 @@ type(tSolutionState) function & Utilities_FFTbackward, & Utilities_updateGamma, & Utilities_constitutiveResponse, & - Utilities_calculateRate + Utilities_calculateRate, & + debugRotation use FEsolving, only: & restartWrite, & restartRead, & @@ -198,10 +199,9 @@ type(tSolutionState) function & real(pReal), dimension(3,3,3,3) :: & S !< current average compliance - real(pReal), dimension(3,3) :: & - F_aim_lab, & - F_aim_lab_lastIter, & !< aim of last iteration - P_av + real(pReal), dimension(3,3) :: & + F_aim_lastIter, & !< aim of last iteration + P_av real(pReal), dimension(3,3,res(1),res(2),res(3)) :: P !-------------------------------------------------------------------------------------------------- ! loop variables, convergence etc. @@ -258,12 +258,12 @@ type(tSolutionState) function & endif if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old F_aim_lastInc = F_aim - Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & - timeinc,timeinc_old,guess,F_lastInc,F) + Fdot = utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & + timeinc_old,guess,F_lastInc,F) F_lastInc = F endif 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) @@ -280,12 +280,14 @@ type(tSolutionState) function & ! report begin of new iteration write(6,'(/,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & ' @ 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 =', & - math_transpose33(F_aim) - + math_transpose33(F_aim) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - F_aim_lab_lastIter = math_rotate_backward33(F_aim,rotation_BC) + F_aim_lastIter = F_aim basic_solution%termIll = .false. call Utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,& P,C,P_av,ForwardData,rotation_BC) @@ -297,8 +299,7 @@ type(tSolutionState) function & ! stress BC handling 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 - 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 field_real = 0.0_pReal @@ -306,7 +307,7 @@ type(tSolutionState) function & order=[4,5,1,2,3]) ! field real has a different order call Utilities_FFTforward() 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() 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) diff --git a/code/DAMASK_spectral_solverBasicPETSc.f90 b/code/DAMASK_spectral_solverBasicPETSc.f90 index 4e28f6425..98fb3f620 100644 --- a/code/DAMASK_spectral_solverBasicPETSc.f90 +++ b/code/DAMASK_spectral_solverBasicPETSc.f90 @@ -270,7 +270,7 @@ type(tSolutionState) function & !-------------------------------------------------------------------------------------------------- ! update coordinates and rate and forward last inc 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)]) endif F_aim = F_aim + f_aimDot * timeinc diff --git a/code/DAMASK_spectral_utilities.f90 b/code/DAMASK_spectral_utilities.f90 index df2d6a537..b24663fd7 100644 --- a/code/DAMASK_spectral_utilities.f90 +++ b/code/DAMASK_spectral_utilities.f90 @@ -46,7 +46,8 @@ module DAMASK_spectral_utilities debugGeneral, & !< general debugging of spectral solver debugDivergence, & !< debugging of divergence calculation (comparison to function used for post processing) 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 @@ -105,7 +106,8 @@ subroutine utilities_init() debug_levelBasic, & debug_spectralDivergence, & debug_spectralRestart, & - debug_spectralFFTW + debug_spectralFFTW, & + debug_spectralRotation use mesh, only: & res, & res1_red, & @@ -132,6 +134,7 @@ subroutine utilities_init() debugDivergence = iand(debug_level(debug_spectral),debug_spectralDivergence) /= 0 debugRestart = iand(debug_level(debug_spectral),debug_spectralRestart) /= 0 debugFFTW = iand(debug_level(debug_spectral),debug_spectralFFTW) /= 0 + debugRotation = iand(debug_level(debug_spectral),debug_spectralRotation) /= 0 !-------------------------------------------------------------------------------------------------- ! allocation @@ -721,10 +724,14 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,& C = C * wgt 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 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 =',& math_transpose33(P_av)/1.e6_pReal end subroutine utilities_constitutiveResponse @@ -733,14 +740,13 @@ end subroutine utilities_constitutiveResponse !-------------------------------------------------------------------------------------------------- !> @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: & res 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) :: & - timeinc, & !< timeinc of current step timeinc_old !< timeinc of last step logical, intent(in) :: & guess !< guess along former trajectory @@ -752,7 +758,7 @@ pure function utilities_calculateRate(delta_aim,timeinc,timeinc_old,guess,field_ if(guess) then utilities_calculateRate = (field-field_lastInc) / timeinc_old 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 end function utilities_calculateRate diff --git a/code/debug.f90 b/code/debug.f90 index 72f90a243..39f19c3c4 100644 --- a/code/debug.f90 +++ b/code/debug.f90 @@ -43,7 +43,8 @@ module debug integer(pInt), parameter, public :: & debug_spectralRestart = debug_maxGeneral*2_pInt**1_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 :: & debug_debug = 1_pInt, & @@ -232,6 +233,8 @@ subroutine debug_init debug_level(what) = ior(debug_level(what), debug_spectralFFTW) case('divergence') debug_level(what) = ior(debug_level(what), debug_spectralDivergence) + case('rotation') + debug_level(what) = ior(debug_level(what), debug_spectralRotation) end select enddo 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_spectralFFTW) /= 0) write(6,'(a)') ' FFTW' 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) endif enddo