From f99bf63397ae2b4db083fba670233323707777e9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 10 Aug 2011 17:45:37 +0000 Subject: [PATCH] removed functions added for debugging of divergence calculation to math.f90 corrected calculation of stress BC condition. Depending on given BC, the stiffness matrix is reduced and than inversed. Then it is filled with zeros and used for the calculation of the correct change of deformation gradient. All calculation is done using dP/dF --- code/DAMASK_spectral.f90 | 403 +++++++-------------------------------- code/math.f90 | 58 ------ 2 files changed, 69 insertions(+), 392 deletions(-) diff --git a/code/DAMASK_spectral.f90 b/code/DAMASK_spectral.f90 index 86874f396..ca6f59aab 100644 --- a/code/DAMASK_spectral.f90 +++ b/code/DAMASK_spectral.f90 @@ -119,25 +119,6 @@ program DAMASK_spectral integer(pInt) i, j, k, l, m, n, p integer(pInt) loadcase, ielem, iter, calcmode, CPFEM_mode, ierr, not_converged_counter logical errmatinv - - !real(pReal) temperature ! not used, but needed for call to CPFEM_general - -!!!!!!!!!!!!!!!!!!!!!!!! start divergence debugging - integer*8 plan_div(3) - real(pReal), dimension(:,:,:,:), allocatable :: divergence - complex(pReal), dimension(:,:,:,:), allocatable :: divergence_hat - complex(pReal), dimension(:,:,:,:), allocatable :: divergence_hat_full - complex(pReal), dimension(:,:,:,:), allocatable :: divergence_hat_full2 - complex(pReal), dimension(:,:,:,:,:), allocatable :: pstress_field_hat, pstress_field - real(pReal) ev1, ev2, ev3 - real(pReal), dimension(3,3) :: evb1, evb2, evb3 - real(pReal) p_hat_avg_inf, p_hat_avg_two, p_real_avg_inf, p_real_avg_two, & - err_div_avg_inf, err_div_avg_two, err_div_max_inf, err_div_max_two, & - err_div_avg_inf2, err_div_avg_two2, err_div_max_two2, err_div_max_inf2, & - err_real_div_avg_inf, err_real_div_avg_two, err_real_div_max_inf, err_real_div_max_two, & - err_real_div_avg_inf2, err_real_div_avg_two2, err_real_div_max_inf2, err_real_div_max_two2, & - rho -!!!!!!!!!!!!!!!!!!!!!!!! end divergence debugging !Initializing !$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS @@ -356,17 +337,7 @@ program DAMASK_spectral 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 (temperature( resolution(1),resolution(2),resolution(3))); temperature = bc_temperature(1) ! start out isothermally - -!!!!!!!!!!!!!!!!!!!!!!!! start divergence debugging -!allocate (xi (3,resolution(1)/2+1,resolution(2),resolution(3))); xi = 0.0_pReal - allocate (xi (3,resolution(1),resolution(2),resolution(3))); xi = 0.0_pReal - allocate (divergence (resolution(1) ,resolution(2),resolution(3),3)); divergence = 0.0_pReal - allocate (divergence_hat (resolution(1)/2+1,resolution(2),resolution(3),3)); divergence_hat = 0.0_pReal - allocate (divergence_hat_full(resolution(1),resolution(2),resolution(3),3)); divergence_hat_full = 0.0_pReal - allocate (divergence_hat_full2(resolution(1),resolution(2),resolution(3),3)); divergence_hat_full2 = 0.0_pReal - allocate (pstress_field_hat(resolution(1),resolution(2),resolution(3),3,3)); pstress_field_hat = 0.0_pReal - allocate (pstress_field (resolution(1),resolution(2),resolution(3),3,3)); pstress_field = 0.0_pReal -!!!!!!!!!!!!!!!!!!!!!!!! end divergence debugging + allocate (xi (3,resolution(1)/2+1,resolution(2),resolution(3))); xi =0.0_pReal wgt = 1.0_pReal/real(resolution(1)*resolution(2)*resolution(3), pReal) defgradAim = math_I3 @@ -392,15 +363,10 @@ program DAMASK_spectral if(k > resolution(3)/2+1) k_s(3) = k_s(3)-resolution(3) do j = 1, resolution(2) k_s(2) = j-1 - if(j > resolution(2)/2+1) k_s(2) = k_s(2)-resolution(2) -!!!!!!!!!!!!!!!!!!!!!!!! start divergence debugging - !do i = 1, resolution(1)/2+1 - ! k_s(1) = i-1 - do i = 1, resolution(1) !defining full xi vector field (no conjugate complex symmetry) + if(j > resolution(2)/2+1) k_s(2) = k_s(2)-resolution(2) + do i = 1, resolution(1)/2+1 k_s(1) = i-1 - if(i > resolution(1)/2+1) k_s(1) = k_s(1)-resolution(1) -!!!!!!!!!!!!!!!!!!!!!!!! end divergence debugging - xi(3,i,j,k) = 0.0_pReal ! 2D case + xi(3,i,j,k) = 0.0_pReal ! 2D case if(resolution(3) > 1) xi(3,i,j,k) = real(k_s(3), pReal)/geomdimension(3) ! 3D case xi(2,i,j,k) = real(k_s(2), pReal)/geomdimension(2) xi(1,i,j,k) = real(k_s(1), pReal)/geomdimension(1) @@ -443,21 +409,6 @@ program DAMASK_spectral call dfftw_plan_many_dft_c2r(plan_fft(2),3,(/resolution(1),resolution(2),resolution(3)/),9,& workfft,(/resolution(1)/2+1,resolution(2),resolution(3)/),1,(resolution(1)/2+1)*resolution(2)*resolution(3),& workfft,(/resolution(1) +2,resolution(2),resolution(3)/),1,(resolution(1) +2)*resolution(2)*resolution(3),FFTW_PATIENT) - -!!!!!!!!!!!!!!!!!!!!!!!! start divergence debugging - call dfftw_plan_many_dft(plan_div(1),3,(/resolution(1),resolution(2),resolution(3)/),9,& - pstress_field, (/resolution(1),resolution(2),resolution(3)/),1,(resolution(1)*resolution(2)*resolution(3)),& - pstress_field_hat,(/resolution(1),resolution(2),resolution(3)/),1,(resolution(1)*resolution(2)*resolution(3)),& - FFTW_FORWARD,FFTW_PATIENT) - call dfftw_plan_many_dft_c2r(plan_div(2),3,(/resolution(1),resolution(2),resolution(3)/),3,& - divergence_hat, (/resolution(1)/2+1,resolution(2),resolution(3)/),1,(resolution(1)/2+1)*resolution(2)*resolution(3),& - divergence ,(/resolution(1), resolution(2),resolution(3)/),1, resolution(1)* resolution(2)*resolution(3),& - FFTW_PATIENT) -call dfftw_plan_many_dft(plan_div(3),3,(/resolution(1),resolution(2),resolution(3)/),3,& - divergence_hat_full,(/resolution(1),resolution(2),resolution(3)/),1,resolution(1)*resolution(2)*resolution(3),& - divergence_hat_full2,(/resolution(1),resolution(2),resolution(3)/),1,resolution(1)* resolution(2)*resolution(3),& - FFTW_BACKWARD,FFTW_PATIENT) -!!!!!!!!!!!!!!!!!!!!!!!! end divergence debugging ! write header of output file open(538,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())& @@ -532,7 +483,7 @@ call dfftw_plan_many_dft(plan_div(3),3,(/resolution(1),resolution(2),resolution( ! update local deformation gradient do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) temp33_Real = defgrad(i,j,k,:,:) - if (velGradApplied(loadcase)) & ! using velocity gradient to calculate new deformation gradient (if not guessing) + if (velGradApplied(loadcase)) & ! using velocity gradient to calculate new deformation gradient (if not guessing) deltaF = math_mul33x33(bc_deformation(:,:,loadcase),defgradold(i,j,k,:,:)) defgrad(i,j,k,:,:) = defgrad(i,j,k,:,:) & ! decide if guessing along former trajectory or apply homogeneous addon (addon only for applied deformation) + guessmode * (defgrad(i,j,k,:,:) - defgradold(i,j,k,:,:))& @@ -551,7 +502,6 @@ call dfftw_plan_many_dft(plan_div(3),3,(/resolution(1),resolution(2),resolution( CPFEM_mode = 1_pInt ! winding forward iter = 0_pInt err_div= 2.0_pReal * err_div_tol ! go into loop - defgradAimCorr = 0.0_pReal ! reset damping calculation !************************************************************* ! convergence loop @@ -602,55 +552,55 @@ call dfftw_plan_many_dft(plan_div(3),3,(/resolution(1),resolution(2),resolution( enddo; enddo err_stress = maxval(abs(mask_stress * (pstress_av - bc_stress(:,:,loadcase)))) - err_stress_tol = maxval(abs(pstress_av))*0.8*err_stress_tolrel + err_stress_tol = maxval(abs(pstress_av)) * 0.8 * err_stress_tolrel - call math_invert(9, math_plain3333to99(c_current),s_current99,i,errmatinv) - if(errmatinv) then - print*, 'using symmetric compliance' - pause !somehow not working, maybe we don't need it - !mask_stress6 = math_Plain33to6_logical(bc_mask(:,:,2,loadcase)) - size_reduced = count(mask_stress6) - allocate (c_reduced(size_reduced,size_reduced)); c_reduced = 0.0_pReal - allocate (s_reduced(size_reduced,size_reduced)); s_reduced = 0.0_pReal - c_current66 = math_Plain3333to66(c_current) - k = 0_pInt - do n = 1,6 - if(mask_stress6(n)) then - k = k + 1_pInt - j = 0_pInt - do m = 1,6 - if(mask_stress6(m)) then - j = j + 1_pInt - c_reduced(k,j) = c_current66(n,m) - endif - enddo - endif - enddo - call math_invert(size_reduced, c_reduced, s_reduced, i, errmatinv) - if(errmatinv) call IO_error(800) - s_current66 = 0.0_pReal - k = 0_pInt - do n = 1,6 - if(mask_stress6(n)) then - k = k + 1_pInt - j = 0_pInt - do m = 1,6 - if(mask_stress6(m)) then - j = j + 1_pInt - s_current66(n,m) = s_reduced(k,j) - endif - enddo - endif - enddo - s_current = math_Plain66to3333(s_current66) - else - print*, 'using non-symmetric compliance' - mask_stress9 = reshape(bc_mask(:,:,2,loadcase),(/9/)) - size_reduced = count(mask_stress9) - allocate (c_reduced(size_reduced,size_reduced)); c_reduced = 0.0_pReal - allocate (s_reduced(size_reduced,size_reduced)); s_reduced = 0.0_pReal - c_current99 = math_Plain3333to99(c_current) + call math_invert(9, math_plain3333to99(c_current),s_current99,i,errmatinv) + if(errmatinv) then + print*, 'using symmetric compliance' + pause !maybe we don't need it. Code below is not working + !mask_stress6 = math_Plain33to6_logical(bc_mask(:,:,2,loadcase)) + !if the symmetrized stiffness is not used at all, the allocate terms only needed at the beginning of a load case + size_reduced = count(mask_stress6) + allocate (c_reduced(size_reduced,size_reduced)); c_reduced = 0.0_pReal + allocate (s_reduced(size_reduced,size_reduced)); s_reduced = 0.0_pReal + ! c_current66 = math_Plain3333to66(c_current) k = 0_pInt + do n = 1,6 + if(mask_stress6(n)) then + k = k + 1_pInt + j = 0_pInt + do m = 1,6 + if(mask_stress6(m)) then + j = j + 1_pInt + c_reduced(k,j) = c_current66(n,m) + endif + enddo + endif + enddo + call math_invert(size_reduced, c_reduced, s_reduced, i, errmatinv) + if(errmatinv) call IO_error(800) + s_current66 = 0.0_pReal + k = 0_pInt + do n = 1,6 + if(mask_stress6(n)) then + k = k + 1_pInt + j = 0_pInt + do m = 1,6 + if(mask_stress6(m)) then + j = j + 1_pInt + s_current66(n,m) = s_reduced(k,j) + endif + enddo + endif + enddo + ! s_current = math_Plain66to3333(s_current66) + else + mask_stress9 = reshape(bc_mask(:,:,2,loadcase),(/9/)) + size_reduced = count(mask_stress9) + allocate (c_reduced(size_reduced,size_reduced)); c_reduced = 0.0_pReal + allocate (s_reduced(size_reduced,size_reduced)); s_reduced = 0.0_pReal + c_current99 = math_Plain3333to99(c_current) + k = 0_pInt ! build reduced stiffness do n = 1,9 if(mask_stress9(n)) then k = k + 1_pInt @@ -659,14 +609,12 @@ call dfftw_plan_many_dft(plan_div(3),3,(/resolution(1),resolution(2),resolution( if(mask_stress9(m)) then j = j + 1_pInt c_reduced(k,j) = c_current99(n,m) - endif - enddo - endif - enddo - call math_invert(size_reduced, c_reduced, s_reduced, i, errmatinv) + endif; enddo; endif; enddo + + call math_invert(size_reduced, c_reduced, s_reduced, i, errmatinv) ! invert reduced stiffness if(errmatinv) call IO_error(800) - - s_current99 = 0.0_pReal + + s_current99 = 0.0_pReal ! build full compliance k = 0_pInt do n = 1,9 if(mask_stress9(n)) then @@ -676,35 +624,19 @@ call dfftw_plan_many_dft(plan_div(3),3,(/resolution(1),resolution(2),resolution( if(mask_stress9(m)) then j = j + 1_pInt s_current99(n,m) = s_reduced(k,j) - endif - enddo - endif - enddo + endif; enddo; endif; enddo s_current = math_Plain99to3333(s_current99) endif + deallocate(c_reduced) deallocate(s_reduced) print*, 'Correcting deformation gradient to fullfill BCs' defgradAimCorr = - math_mul3333xx33(s_current, ((pstress_av - bc_stress(:,:,loadcase)))) ! residual on given stress components - - print*, 'change of Stress:' - print '(3(f10.4,x))', -math_transpose3x3(math_mul3333xx33(c_current,defgradAimCorr))/1.e6 - print*, 'defgradAimCorr:' - print '(3(e12.3,x))', math_transpose3x3(defgradAimCorr) - do m=1,3; do n =1,3 ! calculate damper - ! if ( sign(1.0_pReal,defgradAimCorr(m,n))/=sign(1.0_pReal,defgradAimCorrPrev(m,n))) then - if (defgradAimCorr(m,n) * defgradAimCorrPrev(m,n) < -relevantStrain ** 2.0_pReal) then ! insignificant within relevantstrain around zero - damper(m,n) = max(0.01_pReal,damper(m,n)*0.8) - else - damper(m,n) = min(1.0_pReal,damper(m,n) *1.2) - endif - enddo; enddo - defgradAimCorr = damper * defgradAimCorr - defgradAim = defgradAim + defgradAimCorr + defgradAim = defgradAim + defgradAimCorr do m = 1,3; do n = 1,3 - defgrad(:,:,:,m,n) = defgrad(:,:,:,m,n) + (defgradAim(m,n) - defgrad_av(m,n)) ! anticipated target minus current state + defgrad(:,:,:,m,n) = defgrad(:,:,:,m,n) + defgradAim(m,n) - defgrad_av(m,n) ! anticipated target minus current state enddo; enddo err_div = 2.0_pReal * err_div_tol err_defgrad = maxval(abs(mask_defgrad * (defgrad_av - defgradAim))) @@ -727,9 +659,7 @@ call dfftw_plan_many_dft(plan_div(3),3,(/resolution(1),resolution(2),resolution( cstress,dsde, pstress, dPdF) enddo; enddo; enddo ielem = 0_pInt -!!!!!!!!!!!!!!!!!!!!!!!! start divergence debugging - pstress_field = cmplx(0.0_pReal,0.0_pReal) -!!!!!!!!!!!!!!!!!!!!!!!! end divergence debugging + do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) ielem = ielem + 1_pInt call CPFEM_general(CPFEM_mode,& ! first element in first iteration retains CPFEM_mode 1, @@ -739,9 +669,6 @@ call dfftw_plan_many_dft(plan_div(3),3,(/resolution(1),resolution(2),resolution( cstress,dsde, pstress, dPdF) CPFEM_mode = 2_pInt workfft(i,j,k,:,:) = pstress -!!!!!!!!!!!!!!!!!!!!!!!! start divergence debugging - pstress_field(i,j,k,:,:) = pstress -!!!!!!!!!!!!!!!!!!!!!!!! end divergence debugging cstress_av = cstress_av + math_mandel6to33(cstress) enddo; enddo; enddo cstress_av = cstress_av * wgt @@ -752,195 +679,21 @@ call dfftw_plan_many_dft(plan_div(3),3,(/resolution(1),resolution(2),resolution( print *, 'Calculating equilibrium using spectral method' err_div = 0.0_pReal p_hat_avg = 0.0_pReal -!!!!!!!!!!!!!!!!!!!!!!!! start divergence debugging - - p_hat_avg_inf = 0.0_pReal - p_hat_avg_two = 0.0_pReal - p_real_avg_inf = 0.0_pReal - p_real_avg_two = 0.0_pReal - err_div_avg_inf = 0.0_pReal - err_div_avg_inf2 = 0.0_pReal - err_div_avg_two = 0.0_pReal - err_div_avg_two2 = 0.0_pReal - err_div_max_inf = 0.0_pReal - err_div_max_inf2 = 0.0_pReal - err_div_max_two = 0.0_pReal - err_div_max_two2 = 0.0_pReal - err_real_div_avg_inf = 0.0_pReal - err_real_div_avg_two = 0.0_pReal - err_real_div_max_inf = 0.0_pReal - err_real_div_max_two = 0.0_pReal - err_real_div_avg_inf2 = 0.0_pReal - err_real_div_avg_two2 = 0.0_pReal - err_real_div_max_inf2 = 0.0_pReal - err_real_div_max_two2 = 0.0_pReal -!!!!!!!!!!!!!!!!!!!!!!!! end divergence debugging call dfftw_execute_dft_r2c(plan_fft(1),workfft,workfft) ! FFT of pstress do m = 1,3 ! L infinity norm of stress tensor p_hat_avg = max(p_hat_avg, sum(abs(workfft(1,1,1,:,m)))) ! ignore imaginary part as it is always zero (Nyquist freq for real only input) enddo -!!!!!!!!!!!!!!!!!!!!!!!! start divergence debugging - call dfftw_execute_dft(plan_div(1),pstress_field,pstress_field_hat) - p_hat_avg_inf = p_hat_avg ! using L inf norm as criterion - ! L2 matrix norm, NuMI Skript, LNM, TU Muenchen p. 47, again ignore imaginary part - call math_spectral1(math_mul33x33(workfft(1,1,1,:,:),math_transpose3x3(workfft(1,1,1,:,:))),ev1,ev2,ev3,evb1,evb2,evb3) - rho = max (ev1,ev2,ev3) - p_hat_avg_two = sqrt(rho) -!!!!!!!!!!!!!!!!!!!!!!!! end divergence debugging do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)/2+1 - ! print*, workfft(2*i-1,j,k,:,:)+workfft(i*2,j,k,:,:)*cmplx(0.0,1.0)-pstress_field_hat(i,j,k,:,:) err_div = max(err_div, maxval(abs(math_mul33x3_complex(workfft(i*2-1,j,k,:,:)+& ! maximum of L infinity norm of div(stress), Suquet 2001 - workfft(i*2, j,k,:,:)*img,xi(:,i,j,k)*minval(geomdimension))))) -!!!!!!!!!!!!!!!!!!!!!!!! start divergence debugging - !if(abs(sqrt(sum(math_mul33x3_complex(workfft(i*2-1,j,k,:,:)+workfft(i*2,j,k,:,:)*img,xi(:,i,j,k)*minval(geomdimension)))**2.0))>err_div_max_two) print*, i,j,k - err_div_max_two = max(err_div_max_two,abs(sqrt(sum(math_mul33x3_complex(workfft(i*2-1,j,k,:,:)+& ! maximum of L two norm of div(stress), Suquet 2001 - workfft(i*2,j,k,:,:)*img,xi(:,i,j,k)*minval(geomdimension)))**2.0))) - err_div_avg_inf = err_div_avg_inf + (maxval(abs(math_mul33x3_complex(workfft(i*2-1,j,k,:,:)+& ! sum of squared L infinity norm of div(stress), Suquet 1998 - workfft(i*2,j,k,:,:)*img,xi(:,i,j,k)*minval(geomdimension)))))**2.0 - err_div_avg_two = err_div_avg_two + abs(sum((math_mul33x3_complex(workfft(i*2-1,j,k,:,:)+& ! sum of squared L2 norm of div(stress) ((sqrt())**2 missing), Suquet 1998 - workfft(i*2,j,k,:,:)*img,xi(:,i,j,k)*minval(geomdimension)))**2.0)) -!!!!!!!!!!!!!!!!!!!!!!!! end divergence debugging + workfft(i*2, j,k,:,:)*img,xi(:,i,j,k)*minval(geomdimension))))) enddo; enddo; enddo -!!!!!!!!!!!!!!!!!!!!!!!! start divergence debugging - do i = 0, resolution(1)/2-2 ! reconstruct data of conjugated complex (symmetric) part in Fourier space - m = 1 - do k = 1, resolution(3) - n = 1 - do j = 1, resolution(2) - !print*, xi(:,resolution(1)-i,j,k) + xi(:,i+2,n,m), resolution(1)-i,j,k - err_div_avg_inf = err_div_avg_inf + (maxval(abs(math_mul33x3_complex& - (workfft(3+2*i,n,m,:,:)+workfft(4+i*2,n,m,:,:)*cmplx(0.0,-1.0),xi(:,resolution(1)-i,j,k)*minval(geomdimension)))))**2.0 - err_div_avg_two = err_div_avg_two + abs(sum((math_mul33x3_complex(workfft(3+2*i,n,m,:,:)+workfft(4+i*2,n,m,:,:)*cmplx(0.0,-1.0),& - xi(:,resolution(1)-i,j,k)*minval(geomdimension)))**2.0)) -! print*, workfft(3+2*i,n,m,:,:)+workfft(4+i*2,n,m,:,:)*cmplx(0.0,-1.0)-pstress_field_hat(resolution(1)-i,j,k,:,:) + err_div = err_div/p_hat_avg !weigthting of error by average stress (L infinity norm) - ! workfft(resolution(1)-i,j,k,:,:) = conjg(workfft(2+i,n,m,:,:)) original code for complex array, code above is a little bit confusing because compley data is stored in real array - if(n == 1) n = resolution(2) +1 - n = n-1 - enddo - if(m == 1) m = resolution(3) +1 - m = m -1 - enddo; enddo - ! print*, 'new' - do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) !calculating divergence criteria for full field (no complex symmetry) - if (abs(sqrt(sum(math_mul33x3_complex(pstress_field_hat(i,j,k,:,:),xi(:,i,j,k)*minval(geomdimension)))**2.0))>err_div_max_two2) print*, i,j,k - err_div_max_two2 = max(err_div_max_two2,abs(sqrt(sum(math_mul33x3_complex(pstress_field_hat(i,j,k,:,:),xi(:,i,j,k)*& - minval(geomdimension)))**2.0))) - err_div_max_inf2 = max(err_div_max_inf2,(maxval(abs(math_mul33x3_complex(pstress_field_hat(i,j,k,:,:),xi(:,i,j,k)*& - minval(geomdimension)))))) - err_div_avg_inf2 = err_div_avg_inf2 + (maxval(abs(math_mul33x3_complex(pstress_field_hat(i,j,k,:,:),& - xi(:,i,j,k)*minval(geomdimension)))))**2.0 - err_div_avg_two2 = err_div_avg_two2 + abs(sum((math_mul33x3_complex(pstress_field_hat(i,j,k,:,:),& - xi(:,i,j,k)*minval(geomdimension)))**2.0)) - enddo; enddo; enddo - - err_div_max_inf = err_div ! using L inf norm as criterion, others will be just printed on screen - err_div_max_inf = err_div_max_inf/p_hat_avg_inf - err_div_max_inf2 = err_div_max_inf2/p_hat_avg_inf - err_div_max_two = err_div_max_two/p_hat_avg_two - err_div_max_two2 = err_div_max_two2/p_hat_avg_two - err_div_avg_inf = sqrt(err_div_avg_inf*wgt)/p_hat_avg_inf - err_div_avg_two = sqrt(err_div_avg_two*wgt)/p_hat_avg_two - err_div_avg_inf2 = sqrt(err_div_avg_inf2*wgt)/p_hat_avg_inf - err_div_avg_two2 = sqrt(err_div_avg_two2*wgt)/p_hat_avg_two -!!!!!!!!!!!!!!!!!!!!!!!! end divergence debugging - - err_div = err_div/p_hat_avg !weigthting of error by average stress (L infinity norm) - -!!!!!!!!!!!!!!!!!!!!!!!! start divergence debugging -divergence_hat_full=0.0 -!divergence in real space - do k = 1, resolution(3) ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) - do j = 1, resolution(2) - do i = 1, resolution(1)/2+1 - divergence_hat(i,j,k,1) = (workfft(i*2-1,j,k,1,1)+ workfft(i*2,j,k,1,1)*img)*xi(1,i,j,k)*img*pi*2.0& - + (workfft(i*2-1,j,k,2,1)+ workfft(i*2,j,k,2,1)*img)*xi(2,i,j,k)*img*pi*2.0& - + (workfft(i*2-1,j,k,3,1)+ workfft(i*2,j,k,3,1)*img)*xi(3,i,j,k)*img*pi*2.0 - divergence_hat(i,j,k,2) = (workfft(i*2-1,j,k,1,2)+ workfft(i*2,j,k,1,2)*img)*xi(1,i,j,k)*img*pi*2.0& - + (workfft(i*2-1,j,k,2,2)+ workfft(i*2,j,k,2,2)*img)*xi(2,i,j,k)*img*pi*2.0& - + (workfft(i*2-1,j,k,3,2)+ workfft(i*2,j,k,3,2)*img)*xi(3,i,j,k)*img*pi*2.0 - divergence_hat(i,j,k,3) = (workfft(i*2-1,j,k,1,3)+ workfft(i*2,j,k,1,3)*img)*xi(1,i,j,k)*img*pi*2.0& - + (workfft(i*2-1,j,k,2,3)+ workfft(i*2,j,k,2,3)*img)*xi(2,i,j,k)*img*pi*2.0& - + (workfft(i*2-1,j,k,3,3)+ workfft(i*2,j,k,3,3)*img)*xi(3,i,j,k)*img*pi*2.0 - enddo; enddo; enddo - do k = 1, resolution(3) ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) - do j = 1, resolution(2) - do i = 1, resolution(1)/2+1 - divergence_hat_full(i,j,k,1) = pstress_field_hat(i,j,k,1,1)*xi(1,i,j,k)*img*pi*2.0& - + pstress_field_hat(i,j,k,2,1)*xi(2,i,j,k)*img*pi*2.0& - + pstress_field_hat(i,j,k,3,1)*xi(3,i,j,k)*img*pi*2.0 - divergence_hat_full(i,j,k,2) = pstress_field_hat(i,j,k,1,2)*xi(1,i,j,k)*img*pi*2.0& - + pstress_field_hat(i,j,k,2,2)*xi(2,i,j,k)*img*pi*2.0& - + pstress_field_hat(i,j,k,3,2)*xi(3,i,j,k)*img*pi*2.0 - divergence_hat_full(i,j,k,3) = pstress_field_hat(i,j,k,1,3)*xi(1,i,j,k)*img*pi*2.0& - + pstress_field_hat(i,j,k,2,3)*xi(2,i,j,k)*img*pi*2.0& - + pstress_field_hat(i,j,k,3,3)*xi(3,i,j,k)*img*pi*2.0 - enddo; enddo; enddo - - ! do i = 0, resolution(1)/2-2 ! reconstruct data of conjugated complex (symmetric) part in Fourier spaced - ! m = 1 - ! do k = 1, resolution(3) - ! n = 1 - ! do j = 1, resolution(2) - ! divergence_hat_full(resolution(1)-i,j,k,1) = pstress_field_hat(resolution(1)-i,j,k,1,1)*xi(1,2+i,n,m)*img*pi*2.0& - ! + pstress_field_hat(resolution(1)-i,j,k,2,1)*xi(2,2+i,n,m)*img*pi*2.0& - ! + pstress_field_hat(resolution(1)-i,j,k,3,1)*xi(3,2+i,n,m)*img*pi*2.0 - ! divergence_hat_full(resolution(1)-i,j,k,2) = pstress_field_hat(resolution(1)-i,j,k,1,2)*xi(1,2+i,n,m)*img*pi*2.0& - ! + pstress_field_hat(resolution(1)-i,j,k,2,2)*xi(2,2+i,n,m)*img*pi*2.0& - ! + pstress_field_hat(resolution(1)-i,j,k,3,2)*xi(3,2+i,n,m)*img*pi*2.0 - ! divergence_hat_full(resolution(1)-i,j,k,3) = pstress_field_hat(resolution(1)-i,j,k,1,3)*xi(1,2+i,n,m)*img*pi*2.0& - ! + pstress_field_hat(resolution(1)-i,j,k,2,3)*xi(2,2+i,n,m)*img*pi*2.0& - ! + pstress_field_hat(resolution(1)-i,j,k,3,3)*xi(3,2+i,n,m)*img*pi*2.0 - ! if(n == 1) n = resolution(2) +1 - ! n = n-1 - ! enddo - ! if(m == 1) m = resolution(3) +1 - ! m = m -1 - ! enddo; enddo - print*, divergence_hat_full - - pause - call dfftw_execute_dft_c2r(plan_div(2), divergence_hat, divergence) - call dfftw_execute_dft(plan_div(3), divergence_hat_full, divergence_hat_full2) - - divergence = divergence*wgt - divergence_hat_full2 = divergence_hat_full2*wgt - print*, divergence_hat_full2 - - do m = 1,3 ! L infinity norm of stress tensor - p_real_avg_inf = max(p_real_avg_inf, sum(abs(pstress_av(:,m)))) - enddo - - call math_spectral1(math_mul33x33(pstress_av,math_transpose3x3(pstress_av)),ev1,ev2,ev3,evb1,evb2,evb3) - rho = max (ev1,ev2,ev3) - p_real_avg_two = sqrt(rho) - - do k = 1, resolution(3); do j = 1, resolution(2) ;do i = 1, resolution(1) - err_real_div_max_inf = max(err_real_div_max_inf, maxval(divergence(i,j,k,:))) - err_real_div_max_inf2 = max(err_real_div_max_inf2, maxval(real(divergence_hat_full2(i,j,k,:)))) - err_real_div_max_two = max(err_real_div_max_two, sqrt(sum(divergence(i,j,k,:)**2.0))) - err_real_div_max_two2 = max(err_real_div_max_two2, sqrt(sum(real(divergence_hat_full2(i,j,k,:))**2.0))) - err_real_div_avg_inf = err_real_div_avg_inf + (maxval(divergence(i,j,k,:)))**2.0 - err_real_div_avg_inf2 = err_real_div_avg_inf2 + (maxval(real(divergence_hat_full2(i,j,k,:))))**2.0 - err_real_div_avg_two = err_real_div_avg_two + sum(divergence(i,j,k,:)**2.0) ! don't take square root just to square it again - err_real_div_avg_two2 = err_real_div_avg_two2 + sum(real(divergence_hat_full2(i,j,k,:))**2.0) ! don't take square root just to square it again - enddo; enddo; enddo - - err_real_div_max_inf = err_real_div_max_inf/p_real_avg_inf - err_real_div_max_inf2 = err_real_div_max_inf2/p_real_avg_inf - err_real_div_max_two = err_real_div_max_two/p_real_avg_two - err_real_div_max_two2 = err_real_div_max_two2/p_real_avg_two - err_real_div_avg_inf = sqrt(err_real_div_avg_inf*wgt)/p_real_avg_inf - err_real_div_avg_inf2 = sqrt(err_real_div_avg_inf2*wgt)/p_real_avg_inf - err_real_div_avg_two = sqrt(err_real_div_avg_two*wgt)/p_real_avg_two - err_real_div_avg_two2 = sqrt(err_real_div_avg_two2*wgt)/p_real_avg_two -!!!!!!!!!!!!!!!!!!!!!!!! end divergence debugging - - if(memory_efficient) then ! memory saving version, on-the-fly calculation of gamma_hat + if(memory_efficient) then ! memory saving version, on-the-fly calculation of gamma_hat do k = 1, resolution(3); do j = 1, resolution(2) ;do i = 1, resolution(1)/2+1 if (any(xi(:,i,j,k) /= 0.0_pReal)) then do l = 1,3; do m = 1,3 @@ -980,34 +733,16 @@ divergence_hat_full=0.0 defgrad = defgrad + workfft(1:resolution(1),:,:,:,:)*wgt do m = 1,3; do n = 1,3 defgrad_av(m,n) = sum(defgrad(:,:,:,m,n))*wgt - defgrad(:,:,:,m,n) = defgrad(:,:,:,m,n) + mask_defgrad(m,n)*(defgradAim(m,n) - defgrad_av(m,n)) ! anticipated target minus current state on components with prescribed deformation + defgrad(:,:,:,m,n) = defgrad(:,:,:,m,n) + (defgradAim(m,n) - defgrad_av(m,n)) ! anticipated target minus current state enddo; enddo err_stress = maxval(abs(mask_stress * (pstress_av - bc_stress(:,:,loadcase)))) err_stress_tol = maxval(abs(pstress_av))*err_stress_tolrel ! accecpt relative error specified err_defgrad = maxval(abs(mask_defgrad * (defgrad_av - defgradAim))) - print '(2(a,E8.2))', ' error divergence: ',err_div, ' Tol. = ', err_div_tol -!!!!!!!!!!!!!!!!!!!!!!!! start divergence debugging - print '((a,E16.11))', ' error divergence FT (max,inf): ',err_div_max_inf - print '((a,E16.11))', ' error divergence FT (max,inf2): ',err_div_max_inf2 - print '((a,E16.11))', ' error divergence FT (max,two): ',err_div_max_two - print '((a,E16.11))', ' error divergence FT (max,two2): ',err_div_max_two2 - print '((a,E16.11))', ' error divergence FT (avg,inf): ',err_div_avg_inf - print '((a,E16.11))', ' error divergence FT (avg,inf2): ',err_div_avg_inf2 - print '((a,E16.11))', ' error divergence FT (avg,two): ',err_div_avg_two - print '((a,E16.11))', ' error divergence FT (avg,two2): ',err_div_avg_two2 - print '((a,E8.2))', ' error divergence Real (max,inf): ',err_real_div_max_inf - print '((a,E8.2))', ' error divergence Real (max,inf2): ',err_real_div_max_inf2 - print '((a,E8.2))', ' error divergence Real (max,two): ',err_real_div_max_two - print '((a,E8.2))', ' error divergence Real (max,two2): ',err_real_div_max_two2 - print '((a,E8.2))', ' error divergence Real (avg,inf): ',err_real_div_avg_inf - print '((a,E8.2))', ' error divergence Real (avg,inf2): ',err_real_div_avg_inf2 - print '((a,E8.2))', ' error divergence Real (avg,two): ',err_real_div_avg_two - print '((a,E8.2))', ' error divergence Real (avg,two2): ',err_real_div_avg_two2 -!!!!!!!!!!!!!!!!!!!!!!!! end divergence debugging - print '(2(a,E8.2))', ' error stress: ',err_stress, ' Tol. = ', err_stress_tol - print '(2(a,E8.2))', ' error deformation gradient: ',err_defgrad,' Tol. = ', err_defgrad_tol + print '(2(a,E13.8))', ' error divergence: ',err_div, ' Tol. = ', err_div_tol + print '(2(a,E13.8))', ' error stress: ',err_stress, ' Tol. = ', err_stress_tol + print '(2(a,E13.8))', ' error deformation gradient: ',err_defgrad,' Tol. = ', err_defgrad_tol !ToDo: usefull .and. for err_div? if((err_stress > err_stress_tol .or. err_defgrad > err_defgrad_tol) .and. err_div < err_div_tol) then ! change to calculation of BCs, reset damper etc. calcmode = 0_pInt diff --git a/code/math.f90 b/code/math.f90 index 46109cc66..eb84d99f3 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -1165,25 +1165,6 @@ pure function math_transpose3x3(A) forall (i=1:9) math_Plain33to9(i) = m33(mapPlain(1,i),mapPlain(2,i)) endfunction math_Plain33to9 - - -!******************************************************************** -! convert symmetric 3x3 matrix into plain vector 6x1 -!******************************************************************** - pure function math_Plain33to6(m33) - - use prec, only: pReal,pInt - implicit none - - real(pReal), dimension(3,3), intent(in) :: m33 - real(pReal), dimension(6) :: math_Plain33to6 - integer(pInt) i - - forall (i=1:6) math_Plain33to6(i) = m33(mapMandel(1,i),mapMandel(2,i)) - - endfunction math_Plain33to6 - - !******************************************************************** ! convert Plain 9x1 back to 3x3 matrix !******************************************************************** @@ -1330,45 +1311,6 @@ pure function math_transpose3x3(A) endfunction math_Mandel3333to66 -!******************************************************************** -! convert symmetric 3x3x3x3 tensor into Plain matrix 6x6 -!******************************************************************** - pure function math_Plain3333to66(m3333) - - use prec, only: pReal,pInt - implicit none - - real(pReal), dimension(3,3,3,3), intent(in) :: m3333 - real(pReal), dimension(6,6) :: math_Plain3333to66 - integer(pInt) i,j - - forall (i=1:6,j=1:6) math_Plain3333to66(i,j) = & - m3333(mapMandel(1,i),mapMandel(2,i),mapMandel(1,j),mapMandel(2,j)) - - endfunction math_Plain3333to66 - -!******************************************************************** -! convert Plain matrix 6x6 back to symmetric 3x3x3x3 tensor -!******************************************************************** - pure function math_Plain66to3333(m66) - - use prec, only: pReal,pInt - implicit none - - real(pReal), dimension(6,6), intent(in) :: m66 - real(pReal), dimension(3,3,3,3) :: math_Plain66to3333 - integer(pInt) i,j - - forall (i=1:6,j=1:6) - math_Plain66to3333(mapMandel(1,i),mapMandel(2,i),mapMandel(1,j),mapMandel(2,j)) = m66(i,j) - math_Plain66to3333(mapMandel(2,i),mapMandel(1,i),mapMandel(1,j),mapMandel(2,j)) = m66(i,j) - math_Plain66to3333(mapMandel(1,i),mapMandel(2,i),mapMandel(2,j),mapMandel(1,j)) = m66(i,j) - math_Plain66to3333(mapMandel(2,i),mapMandel(1,i),mapMandel(2,j),mapMandel(1,j)) = m66(i,j) - end forall - - endfunction math_Plain66to3333 - - !******************************************************************** ! convert Mandel matrix 6x6 back to symmetric 3x3x3x3 tensor !********************************************************************