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

corrected div calculation regarding dimension and resolution

This commit is contained in:
Martin Diehl 2012-03-19 13:19:15 +00:00
parent 5263366615
commit 632d57cc31
1 changed files with 13 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

27
code/DAMASK_spectral.f90
View File

@ -571,7 +571,9 @@ if (sum(bc(1:N_Loadcases)%incs)>9000_pInt) stop !discuss with Philip, stop code
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around)
if (divergence_correction) then if (divergence_correction) then
virt_dim = 1.0_pReal do i = 1_pInt, 3_pInt
if (i /= minloc(geomdim,1) .and. i /= maxloc(geomdim,1)) virt_dim = geomdim/geomdim(i)
enddo
else else
virt_dim = geomdim virt_dim = geomdim
endif endif
@ -599,7 +601,7 @@ if (sum(bc(1:N_Loadcases)%incs)>9000_pInt) stop !discuss with Philip, stop code
0.0_pReal,ielem,1_pInt,sigma,dsde,P_real(i,j,k,1:3,1:3),dPdF) 0.0_pReal,ielem,1_pInt,sigma,dsde,P_real(i,j,k,1:3,1:3),dPdF)
C = C + dPdF C = C + dPdF
enddo; enddo; enddo enddo; enddo; enddo
C_ref = C * wgt ! linear reference material stiffness C_ref = C * wgt ! linear reference material stiffness
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate the gamma operator ! calculate the gamma operator
@ -775,7 +777,7 @@ 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
CPFEM_mode = 1_pInt ! winding forward CPFEM_mode = 1_pInt ! winding forward
iter = 0_pInt iter = 0_pInt
err_div = 2.0_pReal * err_div_tol ! go into loop err_div = huge(err_div_tol) ! go into loop
!################################################################################################## !##################################################################################################
! convergence loop (looping over iterations) ! convergence loop (looping over iterations)
@ -818,7 +820,6 @@ enddo; enddo; enddo
temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde, & temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde, &
P_real(i,j,k,1:3,1:3),dPdF) P_real(i,j,k,1:3,1:3),dPdF)
CPFEM_mode = 2_pInt CPFEM_mode = 2_pInt
C = C + dPdF C = C + dPdF
enddo; enddo; enddo enddo; enddo; enddo
@ -903,7 +904,7 @@ enddo; enddo; enddo
do i = 2_pInt, res1_red -1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice. do i = 2_pInt, res1_red -1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice.
err_div_RMS = err_div_RMS & err_div_RMS = err_div_RMS &
+ 2.0_pReal*(sum (real(math_mul33x3_complex(P_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 + 2.0_pReal*(sum (real(math_mul33x3_complex(P_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 xi(1:3,i,j,k))*TWOPIIMG)**2.0_pReal)& ! --> sum squared L_2 norm of vector
+sum(aimag(math_mul33x3_complex(P_fourier(i,j,k,1:3,1:3),& +sum(aimag(math_mul33x3_complex(P_fourier(i,j,k,1:3,1:3),&
xi(1:3,i,j,k))*TWOPIIMG)**2.0_pReal)) xi(1:3,i,j,k))*TWOPIIMG)**2.0_pReal))
enddo enddo
@ -918,27 +919,26 @@ enddo; enddo; enddo
xi(1:3,res1_red,j,k))*TWOPIIMG)**2.0_pReal) xi(1:3,res1_red,j,k))*TWOPIIMG)**2.0_pReal)
enddo; enddo enddo; enddo
err_div_RMS = sqrt(err_div_RMS)*wgt ! RMS in real space calculated with Parsevals theorem from Fourier space 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*sqrt(wgt) >err_div& if(err_div_RMS/pstress_av_L2 > err_div&
.and.iter >2_pInt& .and. err_stress < err_stress_tol) then
.and.err_stress < err_stress_tol) then
write(6,'(a)') 'Increasing divergence, stopping iterations' write(6,'(a)') 'Increasing divergence, stopping iterations'
iter = itmax iter = itmax
endif endif
err_div = err_div_RMS/pstress_av_L2*sqrt(wgt) ! criterion to stop iterations err_div = err_div_RMS/pstress_av_L2 ! criterion to stop iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate additional divergence criteria and report ! calculate additional divergence criteria and report
if(debugDivergence) then ! calculate divergence again if(debugDivergence) then ! calculate divergence again
err_div_max = 0.0_pReal err_div_max = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red
temp3_Complex = math_mul33x3_complex(P_fourier(i,j,k,1:3,1:3),& temp3_Complex = math_mul33x3_complex(P_fourier(i,j,k,1:3,1:3)*wgt,& ! weighting P_fourier
xi(1:3,i,j,k))*TWOPIIMG xi(1:3,i,j,k))*TWOPIIMG
err_div_max = max(err_div_max,sqrt(sum(abs(temp3_Complex)**2.0_pReal))) err_div_max = max(err_div_max,sqrt(sum(abs(temp3_Complex)**2.0_pReal)))
divergence_fourier(i,j,k,1:3) = temp3_Complex ! need divergence NOT squared divergence_fourier(i,j,k,1:3) = temp3_Complex ! need divergence NOT squared
enddo; enddo; enddo enddo; enddo; enddo
call fftw_execute_dft_c2r(plan_divergence,divergence_fourier,divergence_real) call fftw_execute_dft_c2r(plan_divergence,divergence_fourier,divergence_real) ! already weighted
divergence_real = divergence_real*wgt
err_real_div_RMS = 0.0_pReal err_real_div_RMS = 0.0_pReal
err_real_div_max = 0.0_pReal err_real_div_max = 0.0_pReal
max_div_error = 0.0_pReal max_div_error = 0.0_pReal
@ -949,11 +949,10 @@ enddo; enddo; enddo
err_real_div_max=max(err_real_div_max,sqrt(sum(divergence_real(i,j,k,1:3)**2.0_pReal)))! maximum of L two norm of div(stress) in real space err_real_div_max=max(err_real_div_max,sqrt(sum(divergence_real(i,j,k,1:3)**2.0_pReal)))! maximum of L two norm of div(stress) in real space
enddo; enddo; enddo enddo; enddo; enddo
err_real_div_RMS = sqrt(wgt*err_real_div_RMS) ! RMS in real space err_real_div_RMS = sqrt(wgt*err_real_div_RMS) ! RMS in real space
err_div_max = err_div_max*sqrt(wgt)
write(6,'(a,es11.4)') 'error divergence FT RMS = ',err_div_RMS write(6,'(a,es11.4)') 'error divergence FT RMS = ',err_div_RMS
write(6,'(a,es11.4)') 'error divergence FT max = ',err_div_max
write(6,'(a,es11.4)') 'error divergence Real RMS = ',err_real_div_RMS write(6,'(a,es11.4)') 'error divergence Real RMS = ',err_real_div_RMS
write(6,'(a,es11.4)') 'error divergence FT max = ',err_div_max
write(6,'(a,es11.4)') 'error divergence Real max = ',err_real_div_max write(6,'(a,es11.4)') 'error divergence Real max = ',err_real_div_max
write(6,'(a,es11.4)') 'max deviat. from postProc = ',max_div_error write(6,'(a,es11.4)') 'max deviat. from postProc = ',max_div_error
endif endif