diff --git a/code/DAMASK_spectral.f90 b/code/DAMASK_spectral.f90 index 95c38e5d8..4a83013ac 100644 --- a/code/DAMASK_spectral.f90 +++ b/code/DAMASK_spectral.f90 @@ -195,7 +195,6 @@ program DAMASK_spectral real(pReal), dimension(:,:,:,:,:), allocatable :: F, F_lastInc real(pReal), dimension(:,:,:,:), allocatable :: coordinates real(pReal), dimension(:,:,:), allocatable :: temperature - real(pReal), dimension(:,:,:), allocatable :: phase_cont ! phase contrast field: C(x)/C_ref !-------------------------------------------------------------------------------------------------- ! variables storing information for spectral method and FFTW @@ -449,7 +448,6 @@ program DAMASK_spectral allocate (xi (3,res1_red,res(2),res(3)), source = 0.0_pReal) allocate (coordinates( res(1), res(2),res(3),3), source = 0.0_pReal) allocate (temperature( res(1), res(2),res(3)), source = bc(1)%temperature) ! start out isothermally - allocate (phase_cont ( res(1), res(2),res(3)), source = 0.0_pReal) tensorField = fftw_alloc_complex(int(res1_red*res(2)*res(3)*9_pInt,C_SIZE_T)) ! allocate continous data using a C function, C_SIZE_T is of type integer(8) call c_f_pointer(tensorField, P_real, [ res(1)+2_pInt,res(2),res(3),3,3]) ! place a pointer for a real representation on tensorField call c_f_pointer(tensorField, deltaF_real, [ res(1)+2_pInt,res(2),res(3),3,3]) ! place a pointer for a real representation on tensorField @@ -625,30 +623,30 @@ C_ref = C * wgt ! loop oper incs defined in input file for current loadcase !################################################################################################## do inc = 1_pInt, bc(loadcase)%incs - totalIncsCounter = totalIncsCounter + 1_pInt - if(totalIncsCounter >= restartInc) then ! do calculations (otherwise just forwarding) - + totalIncsCounter = totalIncsCounter + 1_pInt + !-------------------------------------------------------------------------------------------------- ! forwarding time - timeinc_old = timeinc - if (bc(loadcase)%logscale == 0_pInt) then ! linear scale - timeinc = bc(loadcase)%time/bc(loadcase)%incs ! only valid for given linear time scale. will be overwritten later in case loglinear scale is used - else - if (loadcase == 1_pInt) then ! 1st loadcase of logarithmic scale - if (inc == 1_pInt) then ! 1st inc of 1st loadcase of logarithmic scale - timeinc = bc(1)%time*(2.0_pReal**real( 1_pInt-bc(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd - else ! not-1st inc of 1st loadcase of logarithmic scale - timeinc = bc(1)%time*(2.0_pReal**real(inc-1_pInt-bc(1)%incs ,pReal)) - endif - else ! not-1st loadcase of logarithmic scale - timeinc = time0 *( (1.0_pReal + bc(loadcase)%time/time0 )**(real( inc,pReal)/& - real(bc(loadcase)%incs ,pReal))& - -(1.0_pReal + bc(loadcase)%time/time0 )**(real( (inc-1_pInt),pReal)/& - real(bc(loadcase)%incs ,pReal)) ) + timeinc_old = timeinc + if (bc(loadcase)%logscale == 0_pInt) then ! linear scale + timeinc = bc(loadcase)%time/bc(loadcase)%incs ! only valid for given linear time scale. will be overwritten later in case loglinear scale is used + else + if (loadcase == 1_pInt) then ! 1st loadcase of logarithmic scale + if (inc == 1_pInt) then ! 1st inc of 1st loadcase of logarithmic scale + timeinc = bc(1)%time*(2.0_pReal**real( 1_pInt-bc(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd + else ! not-1st inc of 1st loadcase of logarithmic scale + timeinc = bc(1)%time*(2.0_pReal**real(inc-1_pInt-bc(1)%incs ,pReal)) endif + else ! not-1st loadcase of logarithmic scale + timeinc = time0 *( (1.0_pReal + bc(loadcase)%time/time0 )**(real( inc,pReal)/& + real(bc(loadcase)%incs ,pReal))& + -(1.0_pReal + bc(loadcase)%time/time0 )**(real( (inc-1_pInt),pReal)/& + real(bc(loadcase)%incs ,pReal)) ) endif - time = time + timeinc + endif + time = time + timeinc + if(totalIncsCounter >= restartInc) then ! do calculations (otherwise just forwarding) if (bc(loadcase)%velGradApplied) then ! calculate deltaF_aim from given L and current F deltaF_aim = timeinc * mask_defgrad * math_mul33x33(bc(loadcase)%deformation, F_aim) else ! deltaF_aim = fDot *timeinc where applicable @@ -729,8 +727,8 @@ C_ref = C * wgt ! report begin of new iteration write(6,'(a)') '' write(6,'(a)') '==================================================================' - write(6,'(5(a,i6.6))') 'Loadcase ',loadcase,' Increment ',inc,'/',bc(loadcase)%incs,& - ' @ Iteration ',iter,'/',itmax + write(6,'(6(a,i6.6))') 'Loadcase ',loadcase,' Inc. ',inc,'/',bc(loadcase)%incs,& + ' @ Iter. ',itmin,' < ',iter,' < ',itmax write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim =',& math_transpose33(F_aim) write(6,'(a)') '' @@ -747,8 +745,6 @@ C_ref = C * wgt P_real(i,j,k,1:3,1:3),dPdF) enddo; enddo; enddo - - P_real = 0.0_pReal ! needed because of the padding for FFTW C = 0.0_pReal ielem = 0_pInt @@ -760,11 +756,10 @@ C_ref = C * wgt temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde, & P_real(i,j,k,1:3,1:3),dPdF) CPFEM_mode = 2_pInt - phase_cont(i,j,k) = sqrt(sum(dPdF*dPdF)/sum(C_ref*C_ref)) C = C + dPdF enddo; enddo; enddo call debug_info() - write(6,'(a,es11.4)') 'Max phase contrast = ',maxval(phase_cont) + !-------------------------------------------------------------------------------------------------- ! copy one component of the stress field to to a single FT and check for mismatch @@ -814,8 +809,6 @@ C_ref = C * wgt P_fourier (1:res1_red,1:res(2), res(3)/2_pInt+1_pInt,1:3,1:3)& = cmplx(0.0_pReal,0.0_pReal,pReal) - - !-------------------------------------------------------------------------------------------------- ! stress BC handling if(size_reduced > 0_pInt) then ! calculate stress BC if applied @@ -864,8 +857,9 @@ C_ref = C * wgt err_div_RMS = sqrt(err_div_RMS)*wgt ! RMS in real space calculated with Parsevals theorem from Fourier space - if (err_div_RMS/pstress_av_L2 > err_div & - .and. err_stress < err_stress_tol) then + if (err_div_RMS/pstress_av_L2 > err_div & + .and. err_stress < err_stress_tol & + .and. iter >= itmin ) then write(6,'(a)') 'Increasing divergence, stopping iterations' iter = itmax endif @@ -1000,8 +994,8 @@ C_ref = C * wgt ! updated deformation gradient do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - F(i,j,k,1:3,1:3) = F(i,j,k,1:3,1:3) - deltaF_real(i,j,k,1:3,1:3)*wgt/phase_cont(i,j,k) ! F(x)^(n+1) = F(x)^(n) + correction; *wgt: correcting for missing normalization - enddo; enddo; enddo ! preconditioning: F(x)^(n+1) = F(x)^(n) + correction/phase_contrast(x) + F(i,j,k,1:3,1:3) = F(i,j,k,1:3,1:3) - deltaF_real(i,j,k,1:3,1:3)*wgt ! F(x)^(n+1) = F(x)^(n) + correction; *wgt: correcting for missing normalization + enddo; enddo; enddo !-------------------------------------------------------------------------------------------------- @@ -1059,8 +1053,11 @@ C_ref = C * wgt enddo ! end looping over loadcases write(6,'(a)') '' write(6,'(a)') '##################################################################' - write(6,'(i6.6,a,i6.6,a)') notConvergedCounter, ' out of ', & - notConvergedCounter + convergedCounter, ' increments did not converge!' + write(6,'(i6.6,a,i6.6,a,f5.1,a)') convergedCounter, ' out of ', & + notConvergedCounter + convergedCounter, ' (', & + real(convergedCounter, pReal)/& + real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, & + ' %) increments converged!' close(538) call fftw_destroy_plan(plan_stress); call fftw_destroy_plan(plan_correction) if (debugDivergence) call fftw_destroy_plan(plan_divergence) @@ -1068,5 +1065,5 @@ C_ref = C * wgt call fftw_destroy_plan(plan_scalarField_forth) call fftw_destroy_plan(plan_scalarField_back) endif - call quit(1_pInt) + call quit(0_pInt) end program DAMASK_spectral diff --git a/code/DAMASK_spectral_AL.f90 b/code/DAMASK_spectral_AL.f90 index 6a6397006..c8b3610ff 100644 --- a/code/DAMASK_spectral_AL.f90 +++ b/code/DAMASK_spectral_AL.f90 @@ -89,7 +89,6 @@ program DAMASK_spectral_AL restartInc use numerics, only: & - err_div_tol, & err_stress_tolrel, & err_stress_tolabs, & rotation_tol, & @@ -158,11 +157,11 @@ program DAMASK_spectral_AL !-------------------------------------------------------------------------------------------------- ! stress, stiffness and compliance average etc. - real(pReal), dimension(3,3) :: P_star_av = 0.0_pReal, & + real(pReal), dimension(3,3) :: P_av = 0.0_pReal, & F_aim = math_I3, F_aim_lastInc = math_I3, lambda_av, & mask_stress, mask_defgrad, deltaF, F_star_av, & F_aim_lab ! quantities rotated to other coordinate system - real(pReal), dimension(3,3,3,3) :: C_inc0, C=0.0_pReal, S_lastInc, C_lastInc ! stiffness and compliance + real(pReal), dimension(3,3,3,3) :: C_inc0, C=0.0_pReal, S_inc0, S_lastInc, C_lastInc ! stiffness and compliance real(pReal), dimension(6) :: sigma ! cauchy stress real(pReal), dimension(6,6) :: dsde ! small strain stiffness real(pReal), dimension(9,9) :: s_prev99, c_prev99 ! compliance and stiffness in matrix notation @@ -172,16 +171,16 @@ program DAMASK_spectral_AL !-------------------------------------------------------------------------------------------------- ! pointwise data type(C_PTR) :: tensorField ! fields in real an fourier space - real(pReal), dimension(:,:,:,:,:), pointer :: lambda_real, F_real ! fields in real space (pointer) - complex(pReal), dimension(:,:,:,:,:), pointer :: lambda_fourier, F_fourier ! fields in fourier space (pointer) - real(pReal), dimension(:,:,:,:,:), allocatable :: F_lastInc, F_star, lambda, P, F_star_lastIter + real(pReal), dimension(:,:,:,:,:), pointer :: tau_real, F_real ! fields in real space (pointer) + complex(pReal), dimension(:,:,:,:,:), pointer :: tau_fourier, F_fourier ! fields in fourier space (pointer) + real(pReal), dimension(:,:,:,:,:), allocatable :: F_lastInc, lambda, P, F_star real(pReal), dimension(:,:,:,:,:,:,:), allocatable :: dPdF real(pReal), dimension(:,:,:,:), allocatable :: coordinates real(pReal), dimension(:,:,:), allocatable :: temperature !-------------------------------------------------------------------------------------------------- ! variables storing information for spectral method and FFTW - type(C_PTR) :: plan_correction, plan_lambda ! plans for fftw + type(C_PTR) :: plan_correction, plan_tau ! plans for fftw real(pReal), dimension(3,3) :: xiDyad ! product of wave vectors real(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat ! gamma operator (field) for spectral method real(pReal), dimension(:,:,:,:), allocatable :: xi ! wave vector field for divergence and for gamma operator @@ -191,7 +190,7 @@ program DAMASK_spectral_AL ! loop variables, convergence etc. real(pReal) :: time = 0.0_pReal, time0 = 0.0_pReal, timeinc = 1.0_pReal, timeinc_old = 0.0_pReal ! elapsed time, begin of interval, time interval real(pReal) :: guessmode, err_stress, err_stress_tol, err_f, err_p, err_crit, & - err_f_point, err_p_point, pstress_av_L2, err_div_rms, err_div + err_f_point, err_p_point, err_nl, err_nl_point real(pReal), dimension(3,3), parameter :: ones = 1.0_pReal, zeroes = 0.0_pReal complex(pReal), dimension(3,3) :: temp33_Complex real(pReal), dimension(3,3) :: temp33_Real @@ -325,6 +324,8 @@ program DAMASK_spectral_AL res = mesh_spectral_getResolution() geomdim = mesh_spectral_getDimension() homog = mesh_spectral_getHomogenization() + + res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c) Npoints = res(1)*res(2)*res(3) wgt = 1.0_pReal/real(Npoints, pReal) @@ -407,16 +408,15 @@ program DAMASK_spectral_AL allocate (P ( res(1), res(2),res(3),3,3), source = 0.0_pReal) allocate (dPdF ( res(1), res(2),res(3),3,3,3,3), source = 0.0_pReal) allocate (F_star ( res(1), res(2),res(3),3,3), source = 0.0_pReal) - allocate (F_star_lastIter ( res(1), res(2),res(3),3,3), source = 0.0_pReal) allocate (F_lastInc ( res(1), res(2),res(3),3,3), source = 0.0_pReal) allocate (lambda ( res(1), res(2),res(3),3,3), source = 0.0_pReal) allocate (xi (3,res1_red,res(2),res(3)), source = 0.0_pReal) allocate (coordinates( res(1), res(2),res(3),3), source = 0.0_pReal) allocate (temperature( res(1), res(2),res(3)), source = bc(1)%temperature) ! start out isothermally tensorField = fftw_alloc_complex(int(res1_red*res(2)*res(3)*9_pInt,C_SIZE_T)) ! allocate continous data using a C function, C_SIZE_T is of type integer(8) - call c_f_pointer(tensorField, lambda_real, [ res(1)+2_pInt,res(2),res(3),3,3]) ! place a pointer for the real representation + call c_f_pointer(tensorField, tau_real, [ res(1)+2_pInt,res(2),res(3),3,3]) ! place a pointer for the real representation call c_f_pointer(tensorField, F_real, [ res(1)+2_pInt,res(2),res(3),3,3]) ! place a pointer for the real representation - call c_f_pointer(tensorField, lambda_fourier, [ res1_red, res(2),res(3),3,3]) ! place a pointer for the complex representation + call c_f_pointer(tensorField, tau_fourier, [ res1_red, res(2),res(3),3,3]) ! place a pointer for the complex representation call c_f_pointer(tensorField, F_fourier, [ res1_red, res(2),res(3),3,3]) ! place a pointer for the complex representation !-------------------------------------------------------------------------------------------------- @@ -431,10 +431,10 @@ program DAMASK_spectral_AL !-------------------------------------------------------------------------------------------------- ! creating plans - plan_lambda = fftw_plan_many_dft_r2c(3,[ res(3),res(2) ,res(1)],9,& ! dimensions , length in each dimension in reversed order - lambda_real,[ res(3),res(2) ,res(1)+2_pInt],& ! input data , physical length in each dimension in reversed order + plan_tau = fftw_plan_many_dft_r2c(3,[ res(3),res(2) ,res(1)],9,& ! dimensions , length in each dimension in reversed order + tau_real,[ res(3),res(2) ,res(1)+2_pInt],& ! input data , physical length in each dimension in reversed order 1, res(3)*res(2)*(res(1)+2_pInt),& ! striding , product of physical lenght in the 3 dimensions - lambda_fourier,[ res(3),res(2) ,res1_red],& + tau_fourier,[ res(3),res(2) ,res1_red],& 1, res(3)*res(2)* res1_red,fftw_planner_flag) plan_correction = fftw_plan_many_dft_c2r(3,[ res(3),res(2) ,res(1)],9,& @@ -449,7 +449,8 @@ program DAMASK_spectral_AL ielem = 0_pInt do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) ielem = ielem + 1_pInt - F_real(i,j,k,1:3,1:3) = math_I3; F_lastInc(i,j,k,1:3,1:3) = math_I3 + F_star(i,j,k,1:3,1:3) = math_I3 + F_lastInc(i,j,k,1:3,1:3) = math_I3 coordinates(i,j,k,1:3) = geomdim/real(res * [i,j,k], pReal) - geomdim/real(2_pInt*res,pReal) call CPFEM_general(3_pInt,coordinates(i,j,k,1:3),math_I3,math_I3,temperature(i,j,k),& 0.0_pReal,ielem,1_pInt,sigma,dsde,P(i,j,k,1:3,1:3),dPdF(i,j,k,1:3,1:3,1:3,1:3)) @@ -462,8 +463,8 @@ program DAMASK_spectral_AL 0.0_pReal,ielem,1_pInt,sigma,dsde,P(i,j,k,1:3,1:3),dPdF(i,j,k,1:3,1:3,1:3,1:3)) C = C + dPdF(i,j,k,1:3,1:3,1:3,1:3) enddo; enddo; enddo - C_inc0 = C * wgt ! linear reference material stiffness - + C_inc0 = C * wgt *3.0_pReal ! linear reference material stiffness + !-------------------------------------------------------------------------------------------------- ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) and remove the given highest frequencies do k = 1_pInt, res(3) @@ -497,24 +498,6 @@ program DAMASK_spectral_AL gamma_hat(1,1,1, 1:3,1:3,1:3,1:3) = 0.0_pReal ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 endif -!-------------------------------------------------------------------------------------------------- -! possible restore deformation gradient from saved state - if (restartInc > 1_pInt) then ! using old values from file - if (debugRestart) write(6,'(a,i6,a)') 'Reading values of increment ',& - restartInc - 1_pInt,' from file' - call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',& - trim(getSolverJobName()),size(F_star)) - read (777,rec=1) F_star - close (777) - F_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = F_star - F_lastInc = F_star - F_aim = 0.0_pReal - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - F_aim = F_aim + F_real(i,j,k,1:3,1:3) ! calculating old average deformation - enddo; enddo; enddo - F_aim = F_aim * wgt - F_aim_lastInc = F_aim - endif !-------------------------------------------------------------------------------------------------- ! write header of output file @@ -601,15 +584,14 @@ program DAMASK_spectral_AL + guessmode * mask_stress * (F_aim - F_aim_lastInc)*timeinc/timeinc_old & + deltaF F_aim_lastInc = temp33_Real - F_star_av = F_aim !-------------------------------------------------------------------------------------------------- ! update local deformation gradient deltaF = math_rotate_backward33(deltaF,bc(loadcase)%rotation) do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - temp33_Real = F_real(i,j,k,1:3,1:3) - F_real(i,j,k,1:3,1:3) = F_real(i,j,k,1:3,1:3) & ! decide if guessing along former trajectory or apply homogeneous addon - + guessmode * (F_real(i,j,k,1:3,1:3) - F_lastInc(i,j,k,1:3,1:3))& ! guessing... + temp33_Real = 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) & ! decide if guessing along former trajectory or apply homogeneous addon + + guessmode * (F_star(i,j,k,1:3,1:3) - F_lastInc(i,j,k,1:3,1:3))& ! guessing... *timeinc/timeinc_old & + (1.0_pReal-guessmode) * deltaF ! if not guessing, use prescribed average deformation where applicable F_lastInc(i,j,k,1:3,1:3) = temp33_Real @@ -617,19 +599,23 @@ program DAMASK_spectral_AL !-------------------------------------------------------------------------------------------------- ! Initialize / Update lambda to useful value - P_star_av = P_star_av + math_mul3333xx33(C*wgt, F_aim-F_aim_lastInc) - lambda_av = 0.0_pReal + P_av = P_av + math_mul3333xx33(C*wgt, F_aim-F_aim_lastInc) + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - lambda(i,j,k,1:3,1:3) = P(i,j,k,1:3,1:3) + math_mul3333xx33(dPdF(i,j,k,1:3,1:3,1:3,1:3), & - F_real(i,j,k,1:3,1:3)-F_lastInc(i,j,k,1:3,1:3)) - lambda_av = lambda_av + lambda(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*wgt, F_aim-F_aim_lastInc) enddo; enddo; enddo - lambda_av=lambda_av*wgt - !-------------------------------------------------------------------------------------------------- !Initialize pointwise data for AL scheme: ToDo: good choice? - F_star(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) + F_star_av = 0.0_pReal + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + F_star_av = F_star_av + F_star(i,j,k,1:3,1:3) + enddo; enddo; enddo + + F_star_av = F_star_av *wgt + write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') 'F* =',& + math_transpose33(F_star_av) +! !-------------------------------------------------------------------------------------------------- ! calculate reduced compliance if(size_reduced > 0_pInt) then ! calculate compliance in case stress BC is applied @@ -660,7 +646,7 @@ program DAMASK_spectral_AL endif; enddo; endif; enddo S_lastInc = (math_Plain99to3333(s_prev99)) endif - + S_inc0 = math_invSym3333(C*wgt) !-------------------------------------------------------------------------------------------------- ! report begin of new increment write(6,'(a)') '##################################################################' @@ -690,9 +676,9 @@ program DAMASK_spectral_AL !-------------------------------------------------------------------------------------------------- ! stress BC handling if(size_reduced > 0_pInt) then ! calculate stress BC if applied - err_stress = maxval(abs(mask_stress * (P_star_av - bc(loadcase)%P))) ! maximum deviaton (tensor norm not applicable) - F_aim = F_aim + math_mul3333xx33(S_lastInc,bc(loadcase)%P- P_star_av) - err_stress_tol = maxval(abs(P_star_av)) * err_stress_tolrel ! don't use any tensor norm because the comparison should be coherent + err_stress = maxval(abs(mask_stress * (P_av - bc(loadcase)%P))) ! maximum deviaton (tensor norm not applicable) + 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 else err_stress_tol = + huge(1.0_pReal) endif @@ -703,42 +689,19 @@ program DAMASK_spectral_AL !-------------------------------------------------------------------------------------------------- ! doing Fourier transform write(6,'(a)') '... spectral method ...............................................' - lambda_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = lambda(1:res(1),1:res(2),1:res(3),1:3,1:3) - call fftw_execute_dft_r2c(plan_lambda,lambda_real,lambda_fourier) - lambda_fourier( res1_red,1:res(2) , 1:res(3) ,1:3,1:3)& + tau_real = 0.0_pReal + do k = 1_pInt, res(3); do j = 1_pInt, res(2) ;do i = 1_pInt, res(1) + tau_real(i,j,k,1:3,1:3) = lambda(i,j,k,1:3,1:3) - math_mul3333xx33(C_inc0,F_star(i,j,k,1:3,1:3)) + enddo; enddo; enddo + call fftw_execute_dft_r2c(plan_tau,tau_real,tau_fourier) + tau_fourier( res1_red,1:res(2) , 1:res(3) ,1:3,1:3)& = cmplx(0.0_pReal,0.0_pReal,pReal) - lambda_fourier(1:res1_red, res(2)/2_pInt+1_pInt,1:res(3) ,1:3,1:3)& + tau_fourier(1:res1_red, res(2)/2_pInt+1_pInt,1:res(3) ,1:3,1:3)& = cmplx(0.0_pReal,0.0_pReal,pReal) if(res(3)>1_pInt) & - lambda_fourier(1:res1_red,1:res(2), res(3)/2_pInt+1_pInt,1:3,1:3)& + tau_fourier(1:res1_red,1:res(2), res(3)/2_pInt+1_pInt,1:3,1:3)& = cmplx(0.0_pReal,0.0_pReal,pReal) !-------------------------------------------------------------------------------------------------- -! calculating RMS divergence criterion in Fourier space - pstress_av_L2 = sqrt(maxval(math_eigenvalues33(math_mul33x33(lambda_av,& ! L_2 norm of average stress (http://mathworld.wolfram.com/SpectralNorm.html) - math_transpose33(lambda_av))))) - err_div_RMS = 0.0_pReal - do k = 1_pInt, res(3); do j = 1_pInt, res(2) - do i = 2_pInt, res1_red -1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice. - err_div_RMS = err_div_RMS & - + 2.0_pReal*(sum (real(math_mul33x3_complex(lambda_fourier(i,j,k,1:3,1:3),& ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2. do not take square root and square again - xi(1:3,i,j,k))*TWOPIIMG)**2.0_pReal)& ! --> sum squared L_2 norm of vector - +sum(aimag(math_mul33x3_complex(lambda_fourier(i,j,k,1:3,1:3),& - xi(1:3,i,j,k))*TWOPIIMG)**2.0_pReal)) - enddo - err_div_RMS = err_div_RMS & ! Those two layers (DC and Nyquist) do not have a conjugate complex counterpart - + sum( real(math_mul33x3_complex(lambda_fourier(1 ,j,k,1:3,1:3),& - xi(1:3,1 ,j,k))*TWOPIIMG)**2.0_pReal)& - + sum(aimag(math_mul33x3_complex(lambda_fourier(1 ,j,k,1:3,1:3),& - xi(1:3,1 ,j,k))*TWOPIIMG)**2.0_pReal)& - + sum( real(math_mul33x3_complex(lambda_fourier(res1_red,j,k,1:3,1:3),& - xi(1:3,res1_red,j,k))*TWOPIIMG)**2.0_pReal)& - + sum(aimag(math_mul33x3_complex(lambda_fourier(res1_red,j,k,1:3,1:3),& - xi(1:3,res1_red,j,k))*TWOPIIMG)**2.0_pReal) - enddo; enddo - - err_div_RMS = sqrt(err_div_RMS)*wgt - err_div = err_div_RMS/pstress_av_L2 -!-------------------------------------------------------------------------------------------------- ! using gamma operator to update F if(memory_efficient) then ! memory saving version, on-the-fly calculation of gamma_hat do k = 1_pInt, res(3); do j = 1_pInt, res(2) ;do i = 1_pInt, res1_red @@ -752,7 +715,7 @@ program DAMASK_spectral_AL gamma_hat(1,1,1, l,u,v,w) = temp33_Real(l,v)*xiDyad(u,w) forall(l = 1_pInt:3_pInt, u = 1_pInt:3_pInt) & temp33_Complex(l,u) = sum(gamma_hat(1,1,1, l,u, 1:3,1:3) *& - lambda_fourier(i,j,k,1:3,1:3)) + tau_fourier(i,j,k,1:3,1:3)) F_fourier(i,j,k,1:3,1:3) = - temp33_Complex endif enddo; enddo; enddo @@ -760,24 +723,39 @@ program DAMASK_spectral_AL do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt,res1_red forall( u = 1_pInt:3_pInt, v = 1_pInt:3_pInt) & temp33_Complex(u,v) = sum(gamma_hat(i,j,k, u,v, 1:3,1:3) *& - lambda_fourier(i,j,k,1:3,1:3)) + tau_fourier(i,j,k,1:3,1:3)) F_fourier(i,j,k, 1:3,1:3) = - temp33_Complex enddo; enddo; enddo endif - F_fourier(1,1,1,1:3,1:3) = cmplx((F_aim_lab - F_star_av)*real(Npoints,pReal),0.0_pReal,pReal) + F_fourier(1,1,1,1:3,1:3) = cmplx((F_aim_lab)*real(Npoints,pReal),0.0_pReal,pReal) !-------------------------------------------------------------------------------------------------- ! doing inverse Fourier transform call fftw_execute_dft_c2r(plan_correction,F_fourier,F_real) ! back transform of fluct deformation gradient - 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_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_av = 0.0_pReal + temp33_real = 0.0_pReal + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + F_star_av = F_star_av + F_star(i,j,k,1:3,1:3) + temp33_real = temp33_real + F_real(i,j,k,1:3,1:3) + enddo; enddo; enddo + + F_star_av = F_star_av *wgt + write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') 'F* =',& + math_transpose33(F_star_av) + temp33_real = temp33_real *wgt + write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') 'F =',& + math_transpose33(temp33_real) !-------------------------------------------------------------------------------------------------- -! - if(callCPFEM) then - write(6,'(a)') '... calling CPFEM to update P(F*) and F*.........................' - F_star_lastIter = F_star + write(6,'(a)') '... calling CPFEM to update P(F*) and F*.........................' + err_nl_point = huge(1.0_pReal) + u = 0_pInt + do while (err_nl_point>1.0e6_pReal .or. u<2_pInt) + u = u + 1_pInt ielem = 0_pInt + err_nl_point = 0.0_pReal + err_nl = 0.0_pReal do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) ielem = ielem + 1_pInt call CPFEM_general(3_pInt,& ! collect cycle @@ -793,83 +771,93 @@ program DAMASK_spectral_AL F_star(i,j,k,1:3,1:3),temperature(i,j,k),timeinc,ielem,1_pInt,& sigma,dsde, P(i,j,k,1:3,1:3), dPdF(i,j,k,1:3,1:3,1:3,1:3)) CPFEM_mode = 2_pInt ! winding forward - temp33_Real = lambda(i,j,k,1:3,1:3) - P(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)) + temp33_Real = P(i,j,k,1:3,1:3) - lambda(i,j,k,1:3,1:3) & + + math_mul3333xx33(C_inc0,F_star(i,j,k,1:3,1:3)-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) + math_mul3333xx33(math_invSym3333(& - C_inc0 + dPdF(i,j,k,1:3,1:3,1:3,1:3)), temp33_Real) + err_nl_point = max(err_nl_point, maxval(abs(temp33_real))) + err_nl = max(err_nl, sqrt(math_mul33xx33(temp33_real,temp33_real))) + + ! temp33_Real = lambda(i,j,k,1:3,1:3) - P(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)) + ! 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(i,j,k,1:3,1:3,1:3,1:3)), temp33_Real) + F_star(i,j,k,1:3,1:3) = F_star(i,j,k,1:3,1:3) - 0.1_pReal * math_mul3333xx33(S_inc0, temp33_Real) enddo; enddo; enddo - else - guesses = guesses +1_pInt - write(6,'(a)')' ... linear approximation for P(F*) and F* ', guesses, ' of ', guessmax + + if(u>1) print*, 'comp wise error of nl eq last it:', err_nl_point/1.0e6_pReal + if(u>1) print*, 'norm error: last it ', err_nl/1.0e6_pReal + + enddo + + + write(6,'(a)') '... update λ..........................' + + err_f = 0.0_pReal + err_f_point = 0.0_pReal + err_p = 0.0_pReal + err_p_point = 0.0_pReal + + F_star_av = 0.0_pReal + lambda_av = 0.0_pReal + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - temp33_Real = lambda(i,j,k,1:3,1:3) - (P(i,j,k,1:3,1:3) + math_mul3333xx33(dPdF(i,j,k,1:3,1:3,1:3,1:3),& - F_star(i,j,k,1:3,1:3) -F_star_lastIter(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)) + 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) & + - F_star(i,j,k,1:3,1:3)) + F_star_av = F_star_av + F_star(i,j,k,1:3,1:3) + lambda_av = lambda_av + lambda(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(i,j,k,1:3,1:3,1:3,1:3)), temp33_Real) + temp33_real = F_star(i,j,k,1:3,1:3) - F_real(i,j,k,1:3,1:3) + err_f_point = max(err_f_point, maxval(abs(temp33_real))) + err_f = max(err_f, sqrt(math_mul33xx33(temp33_real,temp33_real))) + enddo; enddo; enddo - endif - write(6,'(a)') '... update λ..........................' - - err_f = 0.0_pReal - err_f_point = 0.0_pReal - err_p = 0.0_pReal - err_p_point = 0.0_pReal - - F_star_av = 0.0_pReal - P_star_av = 0.0_pReal - lambda_av = 0.0_pReal + ielem = 0_pInt + P_av = 0.0_pReal + C = 0.0_pReal do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - 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) & - - F_star(i,j,k,1:3,1:3)) - F_star_av = F_star_av + F_star(i,j,k,1:3,1:3) - lambda_av = lambda_av + lambda(i,j,k,1:3,1:3) - P_star_av = P_star_av + P(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_point = max(err_f_point, maxval(abs(temp33_real))) - err_f = max(err_f, sqrt(math_mul33xx33(temp33_real,temp33_real))) - - temp33_real = lambda(i,j,k,1:3,1:3) - (P(i,j,k,1:3,1:3) + math_mul3333xx33(dPdF(i,j,k,1:3,1:3,1:3,1:3),& - F_star(i,j,k,1:3,1:3) -F_star_lastIter(i,j,k,1:3,1:3))) + ielem = ielem + 1_pInt + call CPFEM_general(3_pInt,& ! collect cycle + coordinates(i,j,k,1:3), F_lastInc(i,j,k,1:3,1:3),& + F_real(i,j,k,1:3,1:3),temperature(i,j,k),timeinc,ielem,1_pInt,& + sigma,dsde, P(i,j,k,1:3,1:3), dPdF(i,j,k,1:3,1:3,1:3,1:3)) + enddo; enddo; enddo + ielem = 0_pInt + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + ielem = ielem + 1_pInt + call CPFEM_general(2_pInt,& + coordinates(i,j,k,1:3), F_lastInc(i,j,k,1:3,1:3),& + F_real(i,j,k,1:3,1:3),temperature(i,j,k),timeinc,ielem,1_pInt,& + sigma,dsde, P(i,j,k,1:3,1:3), dPdF(i,j,k,1:3,1:3,1:3,1:3)) + P_av = P_av + P(i,j,k,1:3,1:3) + temp33_real = lambda(i,j,k,1:3,1:3) - P(i,j,k,1:3,1:3) err_p_point = max(err_p_point, maxval(abs(temp33_real))) err_p = max(err_p, sqrt(math_mul33xx33(temp33_real,temp33_real))) + C = C + dPdF(i,j,k,1:3,1:3,1:3,1:3) enddo; enddo; enddo F_star_av = F_star_av *wgt write (*,'(a,/,3(3(f12.7,1x)/))',advance='no') 'F* =',& math_transpose33(F_star_av) - P_star_av = P_star_av *wgt - write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'P(F*) / GPa =',& - math_transpose33(P_star_av) /1.e6_pReal + P_av = P_av *wgt + write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'P(F) / GPa =',& + math_transpose33(P_av) /1.e6_pReal lambda_av = lambda_av *wgt write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'λ / GPa =',& math_transpose33(lambda_av) /1.e6_pReal err_f = err_f/sqrt(math_mul33xx33(F_star_av,F_star_av)) - err_p = err_p/sqrt(math_mul33xx33(P_star_av,P_star_av)) + err_p = err_p/sqrt(math_mul33xx33(P_av,P_av)) write(6,'(a,es14.7,es14.7)') 'error F', err_f/1e-4, err_f write(6,'(a,es14.7,es14.7)') 'error P', err_p/1e-3, err_p write(6,'(a,es14.7,es14.7)') 'error stress = ',err_stress/err_stress_tol, err_stress - write(6,'(a,es14.7,es14.7)') 'error divergence = ',err_div/err_div_tol, err_div write(6,*) ' ' - write(6,'(a,es14.7)') 'error divergence FT RMS = ',err_div_RMS write(6,'(a,es14.7)') 'max abs err F', err_f_point write(6,'(a,es14.7)') 'max abs err P', err_p_point - err_crit = max(err_p/1e-3, err_f/1e-4,err_div/err_div_tol,err_stress/err_stress_tol) + err_crit = max(err_p/1e-3, err_f/1e-4,err_stress/err_stress_tol) write(6,'(a,es14.7)') 'critical error', err_crit - if (.not. callCPFEM) then - if(err_crit < 1.0_pReal .or. guesses >= guessmax) callCPFEM = .true. - err_crit =huge(1.0_pReal) - else - if(iter >2 .and. iter< itmax-3) callCPFEM=.false. - guesses = 0_pInt - endif enddo ! end looping when convergency is achieved write(6,'(a)') ' ' @@ -888,15 +876,6 @@ program DAMASK_spectral_AL write(538) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! write result to file endif - if( bc(loadcase)%restartFrequency > 0_pInt .and. & - mod(inc - 1_pInt,bc(loadcase)%restartFrequency) == 0_pInt) then ! at frequency of writing restart information set restart parameter for FEsolving (first call to CPFEM_general will write ToDo: true?) - restartWrite = .true. - write(6,*) 'writing converged results for restart' - call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(F_star)) ! writing deformation gradient field to file - write (777,rec=1) F_star - close (777) - restartInc=totalIncsCounter - endif endif ! end calculation/forwarding enddo ! end looping over incs in current loadcase deallocate(c_reduced) @@ -907,6 +886,6 @@ program DAMASK_spectral_AL write(6,'(i6.6,a,i6.6,a)') notConvergedCounter, ' out of ', & notConvergedCounter + convergedCounter, ' increments did not converge!' close(538) - call fftw_destroy_plan(plan_lambda); call fftw_destroy_plan(plan_correction) + call fftw_destroy_plan(plan_tau); call fftw_destroy_plan(plan_correction) call quit(1_pInt) end program DAMASK_spectral_AL diff --git a/code/spectral_quit.f90 b/code/spectral_quit.f90 index ca8c296ae..5edc2f78a 100644 --- a/code/spectral_quit.f90 +++ b/code/spectral_quit.f90 @@ -53,10 +53,10 @@ subroutine quit(stop_id) write(6,'(a,2(i2.2,a),i2.2)') 'Time: ',dateAndTime(5),':',& dateAndTime(6),':',& dateAndTime(7) - if (stop_id == 1_pInt) stop 1 ! normal termination - if (stop_id <= 0_pInt) then ! trigger regridding - write(6,'(a,i6)') 'restart a', stop_id*(-1_pInt) + if (stop_id == 0_pInt) stop 0 ! normal termination + if (stop_id < 0_pInt) then ! trigger regridding + write(0,'(a,i6)') 'restart a', stop_id*(-1_pInt) stop 2 endif - stop 0 ! error + stop 1 ! error end subroutine