diff --git a/code/DAMASK_spectral_solverAL.f90 b/code/DAMASK_spectral_solverAL.f90 index 873b0c4fe..aad6b9632 100644 --- a/code/DAMASK_spectral_solverAL.f90 +++ b/code/DAMASK_spectral_solverAL.f90 @@ -527,8 +527,8 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! stress BC handling - F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%P_BC))) ! S = 0.0 for no bc - err_stress = maxval(abs(mask_stress * (P_av - params%P_BC))) ! mask = 0.0 for no bc + F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%P_BC))) ! S = 0.0 for no bc + err_stress = maxval(abs(mask_stress * (P_av - params%P_BC))) ! mask = 0.0 for no bc !-------------------------------------------------------------------------------------------------- ! constructing residual diff --git a/code/DAMASK_spectral_solverBasic.f90 b/code/DAMASK_spectral_solverBasic.f90 index d80aa196f..44ec7f164 100644 --- a/code/DAMASK_spectral_solverBasic.f90 +++ b/code/DAMASK_spectral_solverBasic.f90 @@ -228,7 +228,7 @@ type(tSolutionState) function & !-------------------------------------------------------------------------------------------------- ! loop variables, convergence etc. real(pReal) :: err_div, err_stress - integer(pInt) :: iter, row, column + integer(pInt) :: iter logical :: ForwardData !-------------------------------------------------------------------------------------------------- @@ -375,9 +375,9 @@ logical function basic_Converged(err_div,pAvgDiv,err_stress,pAvgStress) err_stress/err_stress_tol ] < 1.0_pReal) write(6,'(/,a,f10.2,a,es11.5,a,es11.4,a)') ' error divergence = ', & - err_div/pAvgDivL2/err_div_tol, ' (',err_div/pAvgDivL2,' / m, tol =',err_div_tol,')' - write(6,'(a,f8.2,a,es11.5,a,es11.4,a)') ' error stress BC = ', & - err_stress/err_stress_tol, ' (',err_stress, ' Pa , tol =',err_stress_tol,')' + err_div/pAvgDivL2/err_div_tol, ' (',err_div/pAvgDivL2,' / m, tol =',err_div_tol,')' + write(6, '(a,f10.2,a,es11.5,a,es11.4,a)') ' error stress BC = ', & + err_stress/err_stress_tol, ' (',err_stress, ' Pa, tol =',err_stress_tol,')' flush(6) end function basic_Converged diff --git a/code/DAMASK_spectral_utilities.f90 b/code/DAMASK_spectral_utilities.f90 index 79a0ea2b0..2f9da6ce5 100644 --- a/code/DAMASK_spectral_utilities.f90 +++ b/code/DAMASK_spectral_utilities.f90 @@ -52,14 +52,17 @@ module DAMASK_spectral_utilities !-------------------------------------------------------------------------------------------------- ! debug divergence - real(pReal), private, dimension(:,:,:,:), pointer :: divergence_real !< scalar field real representation for debugging divergence calculation - complex(pReal),private, dimension(:,:,:,:), pointer :: divergence_fourier !< scalar field real representation for debugging divergence calculation + real(pReal), private, dimension(:,:,:,:), pointer :: divReal !< scalar field real representation for debugging divergence calculation + complex(pReal),private, dimension(:,:,:,:), pointer :: divFourier !< scalar field real representation for debugging divergence calculation !-------------------------------------------------------------------------------------------------- ! plans for FFTW - type(C_PTR), private :: plan_scalarField_forth, plan_scalarField_back !< plans for FFTW in case of debugging the Fourier transform - type(C_PTR), private :: plan_forward, plan_backward !< plans for FFTW - type(C_PTR), private :: plan_divergence !< plan for FFTW in case of debugging divergence calculation + type(C_PTR), private :: & + planForth, & !< FFTW plan P(x) to P(k) + planBack, & !< FFTW plan F(k) to F(x) + planDebugForth, & !< FFTW plan for scalar field (proof that order of usual transform is correct) + planDebugBack, & !< FFTW plan for scalar field inverse (proof that order of usual transform is correct) + planDiv !< plan for FFTW in case of debugging divergence calculation !-------------------------------------------------------------------------------------------------- ! variables controlling debugging @@ -73,14 +76,14 @@ module DAMASK_spectral_utilities !-------------------------------------------------------------------------------------------------- ! derived types - type, public :: tSolutionState !< return type of solution from spectral solver variants + type, public :: tSolutionState !< return type of solution from spectral solver variants logical :: converged = .true. logical :: regrid = .false. logical :: termIll = .false. integer(pInt) :: iterationsNeeded = 0_pInt end type tSolutionState - type, public :: tBoundaryCondition !< set of parameters defining a boundary condition + type, public :: tBoundaryCondition !< set of parameters defining a boundary condition real(pReal), dimension(3,3) :: values = 0.0_pReal real(pReal), dimension(3,3) :: maskFloat = 0.0_pReal logical, dimension(3,3) :: maskLogical = .false. @@ -159,7 +162,7 @@ subroutine utilities_init() tensorField, & !< field cotaining data for FFTW in real and fourier space (in place) scalarField_realC, & !< field cotaining data for FFTW in real space when debugging FFTW (no in place) scalarField_fourierC, & !< field cotaining data for FFTW in fourier space when debugging FFTW (no in place) - divergence !< field cotaining data for FFTW in real and fourier space when debugging divergence (in place) + div !< field cotaining data for FFTW in real and fourier space when debugging divergence (in place) write(6,'(/,a)') ' <<<+- DAMASK_spectral_utilities init -+>>>' write(6,'(a)') ' $Id$' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() @@ -188,8 +191,8 @@ subroutine utilities_init() ! allocation allocate (xi(3,res1_red,res(2),res(3)),source = 0.0_pReal) ! frequencies, only half the size for first dimension tensorField = fftw_alloc_complex(int(res1_red*res(2)*res(3)*9_pInt,C_SIZE_T)) ! allocate aligned data using a C function, C_SIZE_T is of type integer(8) - call c_f_pointer(tensorField, field_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, field_fourier, [ res1_red, res(2),res(3),3,3]) ! place a pointer for a complex representation on tensorField + call c_f_pointer(tensorField, field_real, [res(1)+2_pInt-mod(res(1),2_pInt),res(2),res(3),3,3]) ! place a pointer for a real representation on tensorField + call c_f_pointer(tensorField, field_fourier,[res1_red, res(2),res(3),3,3]) ! place a pointer for a complex representation on tensorField !-------------------------------------------------------------------------------------------------- ! general initialization of FFTW (see manual on fftw.org for more details) @@ -203,29 +206,31 @@ subroutine utilities_init() !-------------------------------------------------------------------------------------------------- ! creating plans for the convolution - plan_forward = fftw_plan_many_dft_r2c(3, [res(3),res(2) ,res(1)], 9,& ! dimensions, logical length in each dimension in reversed order, no. of transforms - field_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 length in the 3 dimensions - field_fourier, [res(3),res(2) ,res1_red],& ! output data, physical length in each dimension in reversed order - 1, res(3)*res(2)* res1_red, fftw_planner_flag) ! striding, product of physical length in the 3 dimensions, planner precision + planForth = fftw_plan_many_dft_r2c(3,[res(3),res(2) ,res(1)], 9, & ! dimensions, logical length in each dimension in reversed order, no. of transforms + field_real,[res(3),res(2) ,res(1)+2_pInt-mod(res(1),2_pInt)], & ! input data, physical length in each dimension in reversed order + 1, res(3)*res(2)*(res(1)+2_pInt-mod(res(1),2_pInt)), & ! striding, product of physical length in the 3 dimensions + field_fourier,[res(3),res(2) ,res1_red], & ! output data, physical length in each dimension in reversed order + 1, res(3)*res(2)* res1_red, fftw_planner_flag) ! striding, product of physical length in the 3 dimensions, planner precision - plan_backward = fftw_plan_many_dft_c2r(3, [res(3),res(2) ,res(1)], 9,& ! dimensions, logical length in each dimension in reversed order, no. of transforms - field_fourier, [res(3),res(2) ,res1_red],& ! input data, physical length in each dimension in reversed order - 1, res(3)*res(2)* res1_red,& ! striding, product of physical length in the 3 dimensions - field_real, [res(3),res(2) ,res(1)+2_pInt],& ! output data, physical length in each dimension in reversed order - 1, res(3)*res(2)*(res(1)+2_pInt), fftw_planner_flag) ! striding, product of physical length in the 3 dimensions, planner precision + planBack = fftw_plan_many_dft_c2r(3,[res(3),res(2) ,res(1)], 9, & ! dimensions, logical length in each dimension in reversed order, no. of transforms + field_fourier,[res(3),res(2) ,res1_red], & ! input data, physical length in each dimension in reversed order + 1, res(3)*res(2)* res1_red, & ! striding, product of physical length in the 3 dimensions + field_real,[res(3),res(2) ,res(1)+2_pInt-mod(res(1),2_pInt)], & ! output data, physical length in each dimension in reversed order + 1, res(3)*res(2)*(res(1)+2_pInt-mod(res(1),2_pInt)), & ! striding, product of physical length in the 3 dimensions + fftw_planner_flag) ! planner precision !-------------------------------------------------------------------------------------------------- ! depending on debug options, allocate more memory and create additional plans if (debugDivergence) then - divergence = fftw_alloc_complex(int(res1_red*res(2)*res(3)*3_pInt,C_SIZE_T)) - call c_f_pointer(divergence, divergence_real, [ res(1)+2_pInt,res(2),res(3),3]) - call c_f_pointer(divergence, divergence_fourier, [ res1_red, res(2),res(3),3]) - plan_divergence = fftw_plan_many_dft_c2r(3,[ res(3),res(2) ,res(1)],3,& - divergence_fourier,[ res(3),res(2) ,res1_red],& - 1, res(3)*res(2)* res1_red,& - divergence_real,[ res(3),res(2) ,res(1)+2_pInt],& - 1, res(3)*res(2)*(res(1)+2_pInt),fftw_planner_flag) + div = fftw_alloc_complex(int(res1_red*res(2)*res(3)*3_pInt,C_SIZE_T)) + call c_f_pointer(div,divReal, [res(1)+2_pInt-mod(res(1),2_pInt),res(2),res(3),3]) + call c_f_pointer(div,divFourier,[res1_red, res(2),res(3),3]) + planDiv = fftw_plan_many_dft_c2r(3,[res(3),res(2) ,res(1)],3,& + divFourier,[res(3),res(2) ,res1_red],& + 1, res(3)*res(2)* res1_red,& + divReal,[res(3),res(2) ,res(1)+2_pInt-mod(res(1),2_pInt)], & + 1, res(3)*res(2)*(res(1)+2_pInt-mod(res(1),2_pInt)), & + fftw_planner_flag) endif if (debugFFTW) then @@ -233,9 +238,9 @@ subroutine utilities_init() scalarField_fourierC = fftw_alloc_complex(int(res(1)*res(2)*res(3),C_SIZE_T)) ! allocate data for fourier representation (no in place transform) call c_f_pointer(scalarField_realC, scalarField_real, [res(1),res(2),res(3)]) ! place a pointer for a real representation call c_f_pointer(scalarField_fourierC, scalarField_fourier, [res(1),res(2),res(3)]) ! place a pointer for a fourier representation - plan_scalarField_forth = fftw_plan_dft_3d(res(3),res(2),res(1),& ! reversed order (C style) + planDebugForth = fftw_plan_dft_3d(res(3),res(2),res(1),& ! reversed order (C style) scalarField_real,scalarField_fourier,-1,fftw_planner_flag) ! input, output, forward FFT(-1), planner precision - plan_scalarField_back = fftw_plan_dft_3d(res(3),res(2),res(1),& ! reversed order (C style) + planDebugBack = fftw_plan_dft_3d(res(3),res(2),res(1),& ! reversed order (C style) scalarField_fourier,scalarField_real,+1,fftw_planner_flag) ! input, output, backward (1), planner precision endif @@ -326,7 +331,7 @@ end subroutine utilities_updateGamma !> In case of debugging the FFT, also one component of the tensor (specified by row and column) !> is independetly transformed complex to complex and compared to the whole tensor transform !-------------------------------------------------------------------------------------------------- -subroutine utilities_FFTforward(row,column) +subroutine utilities_FFTforward() !< @ToDo make row and column between randomly between 1 and 3 use math use mesh, only : & scaledDim, & @@ -334,42 +339,55 @@ subroutine utilities_FFTforward(row,column) res1_red implicit none - integer(pInt), intent(in), optional :: row, column !< if debug FFTW, compare 3D array field of row and column - + integer(pInt) :: row, column ! if debug FFTW, compare 3D array field of row and column + integer(pInt), dimension(2:3,2) :: Nyquist ! highest frequencies to be removed (1 if even, 2 if odd) + real(pReal), dimension(2) :: myRand ! random numbers + !-------------------------------------------------------------------------------------------------- ! copy one component of the stress field to to a single FT and check for mismatch - if (debugFFTW .and. present(row) .and. present(column)) & - scalarField_real(1:res(1),1:res(2),1:res(3)) =& ! store the selected component - cmplx(field_real(1:res(1),1:res(2),1:res(3),row,column),0.0_pReal,pReal) + if (debugFFTW) then + call random_number(myRand) ! two numbers: 0 <= x < 1 + row = nint(myRand(1)*2_pReal + 1_pReal,pInt) + column = nint(myRand(2)*2_pReal + 1_pReal,pInt) + scalarField_real = cmplx(field_real(1:res(1),1:res(2),1:res(3),row,column),0.0_pReal,pReal) ! store the selected component + endif !-------------------------------------------------------------------------------------------------- ! doing the FFT - call fftw_execute_dft_r2c(plan_forward,field_real,field_fourier) + call fftw_execute_dft_r2c(planForth,field_real,field_fourier) !-------------------------------------------------------------------------------------------------- ! comparing 1 and 3x3 FT results - if (debugFFTW .and. present(row) .and. present(column)) then - call fftw_execute_dft(plan_scalarField_forth,scalarField_real,scalarField_fourier) - write(6,'(/,a,i1,1x,i1,a)') ' .. checking FT results of compontent ', row, column, ' ..' - flush(6) - write(6,'(/,a,2(es11.4,1x))') ' max FT relative error = ',& ! print real and imaginary part seperately - maxval( real((scalarField_fourier(1:res1_red,1:res(2),1:res(3))-& - field_fourier(1:res1_red,1:res(2),1:res(3),row,column))/& - scalarField_fourier(1:res1_red,1:res(2),1:res(3)))), & - maxval(aimag((scalarField_fourier(1:res1_red,1:res(2),1:res(3))-& - field_fourier(1:res1_red,1:res(2),1:res(3),row,column))/& - scalarField_fourier(1:res1_red,1:res(2),1:res(3)))) - flush(6) - endif + if (debugFFTW) then + call fftw_execute_dft(planDebugForth,scalarField_real,scalarField_fourier) + where(abs(scalarField_fourier(1:res1_red,1:res(2),1:res(3))) > tiny(1.0_pReal)) ! avoid division by zero + scalarField_fourier(1:res1_red,1:res(2),1:res(3)) = & + (scalarField_fourier(1:res1_red,1:res(2),1:res(3))-& + field_fourier(1:res1_red,1:res(2),1:res(3),row,column))/& + scalarField_fourier(1:res1_red,1:res(2),1:res(3)) + else where + scalarField_real = cmplx(0.0,0.0,pReal) + end where + write(6,'(/,a,i1,1x,i1,a)') ' .. checking FT results of compontent ', row, column, ' ..' + write(6,'(/,a,2(es11.4,1x))') ' max FT relative error = ',& ! print real and imaginary part seperately + maxval(real (scalarField_fourier(1:res1_red,1:res(2),1:res(3)))),& + maxval(aimag(scalarField_fourier(1:res1_red,1:res(2),1:res(3)))) + flush(6) + endif !-------------------------------------------------------------------------------------------------- ! removing highest frequencies - field_fourier ( res1_red,1:res(2) , 1:res(3) ,1:3,1:3)& + Nyquist(2,1:2) = [res(2)/2_pInt + 1_pInt, res(2)/2_pInt + 1_pInt + mod(res(2),2_pInt)] + Nyquist(3,1:2) = [res(3)/2_pInt + 1_pInt, res(3)/2_pInt + 1_pInt + mod(res(3),2_pInt)] + + if(res(1)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation + field_fourier (res1_red, 1:res(2), 1:res(3), 1:3,1:3) & = cmplx(0.0_pReal,0.0_pReal,pReal) - field_fourier (1:res1_red, res(2)/2_pInt+1_pInt,1:res(3) ,1:3,1:3)& + if(res(2)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation + field_fourier (1:res1_red,Nyquist(2,1):Nyquist(2,2),1:res(3), 1:3,1:3) & = cmplx(0.0_pReal,0.0_pReal,pReal) - if(res(3)>1_pInt) & ! do not delete the whole slice in case of 2D calculation - field_fourier (1:res1_red,1:res(2), res(3)/2_pInt+1_pInt,1:3,1:3)& + if(res(3)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation + field_fourier (1:res1_red,1:res(2), Nyquist(3,1):Nyquist(3,2),1:3,1:3) & = cmplx(0.0_pReal,0.0_pReal,pReal) end subroutine utilities_FFTforward @@ -381,7 +399,7 @@ end subroutine utilities_FFTforward !> is independetly transformed complex to complex and compared to the whole tensor transform !> results is weighted by number of points stored in wgt !-------------------------------------------------------------------------------------------------- -subroutine utilities_FFTbackward(row,column) +subroutine utilities_FFTbackward() use math !< must use the whole module for use of FFTW use mesh, only: & wgt, & @@ -389,14 +407,19 @@ subroutine utilities_FFTbackward(row,column) res1_red implicit none - integer(pInt), intent(in), optional :: row, column !< if debug FFTW, compare 3D array field of row and column + integer(pInt) :: row, column !< if debug FFTW, compare 3D array field of row and column integer(pInt) :: i, j, k, m, n - + real(pReal), dimension(2) :: myRand + !-------------------------------------------------------------------------------------------------- ! unpack FFT data for conj complex symmetric part. This data is not transformed when using c2r - if (debugFFTW .and. present(row) .and. present(column)) then - scalarField_fourier = field_fourier(1:res1_red,1:res(2),1:res(3),row,column) - do i = 0_pInt, res(1)/2_pInt-2_pInt + if (debugFFTW) then + call random_number(myRand) ! two numbers: 0 <= x < 1 + row = nint(myRand(1)*2_pReal + 1_pReal,pInt) + column = nint(myRand(2)*2_pReal + 1_pReal,pInt) + scalarField_fourier(1:res1_red,1:res(2),1:res(3)) & + = field_fourier(1:res1_red,1:res(2),1:res(3),row,column) + do i = 0_pInt, res(1)/2_pInt-2_pInt + mod(res(1),2_pInt) m = 1_pInt do k = 1_pInt, res(3) n = 1_pInt @@ -412,22 +435,25 @@ subroutine utilities_FFTbackward(row,column) !-------------------------------------------------------------------------------------------------- ! doing the iFFT - call fftw_execute_dft_c2r(plan_backward,field_fourier,field_real) ! back transform of fluct deformation gradient + call fftw_execute_dft_c2r(planBack,field_fourier,field_real) ! back transform of fluct deformation gradient !-------------------------------------------------------------------------------------------------- ! comparing 1 and 3x3 inverse FT results - if (debugFFTW .and. present(row) .and. present(column)) then - write(6,'(/,a,i1,1x,i1,a)') ' ... checking iFT results of compontent ', row, column, ' ..' - flush(6) - call fftw_execute_dft(plan_scalarField_back,scalarField_fourier,scalarField_real) - write(6,'(/,a,es11.4)') ' max iFT relative error = ',& - maxval((real(scalarField_real(1:res(1),1:res(2),1:res(3)))-& - field_real(1:res(1),1:res(2),1:res(3),row,column))/& - real(scalarField_real(1:res(1),1:res(2),1:res(3)))) - flush(6) - endif - - field_real = field_real * wgt ! normalize the result by number of elements + if (debugFFTW) then + call fftw_execute_dft(planDebugBack,scalarField_fourier,scalarField_real) + where(abs(real(scalarField_real,pReal)) > tiny(1.0_pReal)) ! avoid division by zero + scalarField_real = (scalarField_real & + - cmplx(field_real(1:res(1),1:res(2),1:res(3),row,column), 0.0, pReal))/ & + scalarField_real + else where + scalarField_real = cmplx(0.0,0.0,pReal) + end where + write(6,'(/,a,i1,1x,i1,a)') ' ... checking iFT results of compontent ', row, column, ' ..' + write(6,'(/,a,es11.4)') ' max iFT relative error = ', maxval(real(scalarField_real,pReal)) + flush(6) + endif + + field_real = field_real * wgt ! normalize the result by number of elements end subroutine utilities_FFTbackward @@ -507,46 +533,47 @@ real(pReal) function utilities_divergenceRMS() write(6,'(/,a)') ' ... calculating divergence ................................................' flush(6) + !-------------------------------------------------------------------------------------------------- ! calculating RMS divergence criterion in Fourier space utilities_divergenceRMS = 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. + do i = 2_pInt, res1_red -1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice. utilities_divergenceRMS = utilities_divergenceRMS & - + 2.0_pReal*(sum (real(math_mul33x3_complex(field_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 + + 2.0_pReal*(sum (real(math_mul33x3_complex(field_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(field_fourier(i,j,k,1:3,1:3),& xi(1:3,i,j,k))*TWOPIIMG)**2.0_pReal)) enddo - utilities_divergenceRMS = utilities_divergenceRMS & ! these two layers (DC and Nyquist) do not have a conjugate complex counterpart - + sum( real(math_mul33x3_complex(field_fourier(1 ,j,k,1:3,1:3),& - xi(1:3,1 ,j,k))*TWOPIIMG)**2.0_pReal)& - + sum(aimag(math_mul33x3_complex(field_fourier(1 ,j,k,1:3,1:3),& - xi(1:3,1 ,j,k))*TWOPIIMG)**2.0_pReal)& - + sum( real(math_mul33x3_complex(field_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(field_fourier(res1_red,j,k,1:3,1:3),& - xi(1:3,res1_red,j,k))*TWOPIIMG)**2.0_pReal) + utilities_divergenceRMS = utilities_divergenceRMS & ! these two layers (DC and Nyquist) do not have a conjugate complex counterpart (if res(1) /= 1) + + sum( real(math_mul33x3_complex(field_fourier(1 ,j,k,1:3,1:3),& + xi(1:3,1 ,j,k))*TWOPIIMG)**2.0_pReal)& + + sum(aimag(math_mul33x3_complex(field_fourier(1 ,j,k,1:3,1:3),& + xi(1:3,1 ,j,k))*TWOPIIMG)**2.0_pReal)& + + sum( real(math_mul33x3_complex(field_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(field_fourier(res1_red,j,k,1:3,1:3),& + xi(1:3,res1_red,j,k))*TWOPIIMG)**2.0_pReal) enddo; enddo + if(res(1) == 1_pInt) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of res(1) == 1 + utilities_divergenceRMS = sqrt(utilities_divergenceRMS) * wgt ! RMS in real space calculated with Parsevals theorem from Fourier space - utilities_divergenceRMS = sqrt(utilities_divergenceRMS) *wgt ! RMS in real space calculated with Parsevals theorem from Fourier space - !-------------------------------------------------------------------------------------------------- ! calculate additional divergence criteria and report - if (debugDivergence) then ! calculate divergence again + if (debugDivergence) then ! calculate divergence again err_div_max = 0.0_pReal do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red - temp3_Complex = math_mul33x3_complex(field_fourier(i,j,k,1:3,1:3)*wgt,& ! weighting P_fourier + temp3_Complex = math_mul33x3_complex(field_fourier(i,j,k,1:3,1:3)*wgt,& ! weighting P_fourier xi(1:3,i,j,k))*TWOPIIMG err_div_max = max(err_div_max,sum(abs(temp3_Complex)**2.0_pReal)) - divergence_fourier(i,j,k,1:3) = temp3_Complex ! need divergence NOT squared + divFourier(i,j,k,1:3) = temp3_Complex ! need divergence NOT squared enddo; enddo; enddo - call fftw_execute_dft_c2r(plan_divergence,divergence_fourier,divergence_real) ! already weighted + call fftw_execute_dft_c2r(planDiv,divFourier,divReal) ! already weighted - err_real_div_RMS = sqrt(wgt*sum(divergence_real**2.0_pReal)) ! RMS in real space - err_real_div_max = sqrt(maxval(sum(divergence_real**2.0_pReal,dim=4))) ! max in real space - err_div_max = sqrt( err_div_max) ! max in Fourier space + err_real_div_RMS = sqrt(wgt*sum(divReal**2.0_pReal)) ! RMS in real space + err_real_div_max = sqrt(maxval(sum(divReal**2.0_pReal,dim=4))) ! max in real space + err_div_max = sqrt( err_div_max) ! max in Fourier space write(6,'(/,1x,a,es11.4)') 'error divergence FT RMS = ',utilities_divergenceRMS write(6,'(1x,a,es11.4)') 'error divergence Real RMS = ',err_real_div_RMS @@ -763,10 +790,10 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,& integer(pInt) :: & calcMode, & !< CPFEM mode for calculation - collectMode !< CPFEM mode for collection + collectMode, & !< CPFEM mode for collection + j,k real(pReal), dimension(3,3,3,3) :: max_dPdF, min_dPdF - real(pReal) :: max_dPdF_norm, min_dPdF_norm, defgradDetMin, defgradDetMax, defgradDet - integer(pInt) :: i,j,k + real(pReal) :: max_dPdF_norm, min_dPdF_norm, defgradDetMin, defgradDetMax, defgradDet write(6,'(/,a)') ' ... evaluating constitutive response ......................................' calcMode = CPFEM_CALCRESULTS @@ -939,11 +966,11 @@ subroutine utilities_destroy() implicit none - if (debugDivergence) call fftw_destroy_plan(plan_divergence) - if (debugFFTW) call fftw_destroy_plan(plan_scalarField_forth) - if (debugFFTW) call fftw_destroy_plan(plan_scalarField_back) - call fftw_destroy_plan(plan_forward) - call fftw_destroy_plan(plan_backward) + if (debugDivergence) call fftw_destroy_plan(planDiv) + if (debugFFTW) call fftw_destroy_plan(planDebugForth) + if (debugFFTW) call fftw_destroy_plan(planDebugBack) + call fftw_destroy_plan(planForth) + call fftw_destroy_plan(planBack) end subroutine utilities_destroy diff --git a/code/config/numerics.config b/code/config/numerics.config index efeb5ff03..0d71ed875 100644 --- a/code/config/numerics.config +++ b/code/config/numerics.config @@ -66,7 +66,7 @@ itmin 2 # Minimum iteration number maxCutBack 3 # maximum cut back level (0: 1, 1: 0.5, 2: 0.25, etc) memory_efficient 1 # Precalculate Gamma-operator (81 double per point) update_gamma 0 # Update Gamma-operator with current dPdF (not possible if memory_efficient=1) -divergence_correction 2 # Use dimension-independent divergence criterion +divergence_correction 2 # Use size-independent divergence criterion myspectralsolver basic # Type of spectral solver (basic: basic, basicPETSc: basic with PETSc, AL: augmented Lagrange) myfilter none # Type of filtering method to mitigate Gibb's phenomenon (none, cosine, ...) petsc_options -snes_type ngmres -snes_ngmres_anderson # PetSc solver options diff --git a/code/mesh.f90 b/code/mesh.f90 index 3aaae97ae..c68d41680 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -444,9 +444,9 @@ subroutine mesh_init(ip,el) integer(pInt), intent(in) :: el, ip integer(pInt) :: j - write(6,'(/,a)') ' <<<+- mesh init -+>>>' - write(6,'(a)') ' $Id$' - write(6,'(a16,a)') ' Current time : ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- mesh init -+>>>' + write(6,'(a)') ' $Id$' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" if (allocated(mesh_mapFEtoCPelem)) deallocate(mesh_mapFEtoCPelem) @@ -486,17 +486,12 @@ subroutine mesh_init(ip,el) do j = 1_pInt, 3_pInt if (j /= minloc(geomdim/res,1) .and. j /= maxloc(geomdim/res,1)) scaledDim = geomdim/geomdim(j)*res(j) enddo - elseif (divergence_correction == 3_pInt) then - do j = 1_pInt, 3_pInt - if (j/=minloc(geomdim/sqrt(real(res,pReal)),1) .and. j/=maxloc(geomdim/sqrt(real(res,pReal)),1))& - scaledDim=geomdim/geomdim(j)*sqrt(real(res(j),pReal)) - enddo else scaledDim = geomdim endif - write(6,'(a,3(i12 ))') ' resolution a b c:', res - write(6,'(a,3(f12.5))') ' dimension x y z:', geomdim - write(6,'(a,i5,/)') ' homogenization: ', homog + write(6,'(a,3(i12 ))') ' grid a b c: ', res + write(6,'(a,3(f12.5))') ' size x y z: ', geomdim + write(6,'(a,i5,/)') ' homogenization: ', homog call mesh_spectral_count_nodesAndElements call mesh_spectral_count_cpElements call mesh_spectral_map_elements @@ -566,8 +561,8 @@ end subroutine mesh_init !! valid questions (what) are 'elem', 'node' !-------------------------------------------------------------------------------------------------- integer(pInt) function mesh_FEasCP(what,myID) - - use IO, only: IO_lc + use IO, only: & + IO_lc implicit none character(len=*), intent(in) :: what @@ -654,11 +649,11 @@ end subroutine mesh_build_subNodeCoords !> @brief Calculates IP volume. Allocates global array 'mesh_ipVolume' !-------------------------------------------------------------------------------------------------- subroutine mesh_build_ipVolumes + use math, only: & + math_volTetrahedron, & + math_areaTriangle - use math, only: math_volTetrahedron, & - math_areaTriangle implicit none - integer(pInt) :: e,f,t,i,j,n, & NmySubNodes, & !< number of subnodes already found for this (2d) ip cell Ntriangles !< each interface is made up of this many triangles @@ -740,8 +735,8 @@ end subroutine mesh_build_ipVolumes ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !-------------------------------------------------------------------------------------------------- subroutine mesh_build_ipCoordinates - - use prec, only: tol_gravityNodePos + use prec, only: & + tol_gravityNodePos implicit none integer(pInt) :: e,f,t,i,j,k,n @@ -786,11 +781,10 @@ end subroutine mesh_build_ipCoordinates !> @brief Calculates cell center coordinates. !-------------------------------------------------------------------------------------------------- pure function mesh_cellCenterCoordinates(i,e) - -use prec, only: tol_gravityNodePos +use prec, only: & + tol_gravityNodePos implicit none - !*** input variables integer(pInt), intent(in) :: e, & ! element number i ! integration point number @@ -900,9 +894,7 @@ function mesh_spectral_getGrid(fileUnit) if (.not. gotGrid) & call IO_error(error_ID = 845_pInt, ext_msg='grid') - if(mesh_spectral_getGrid(1) < 2_pInt .or. & ! must be at least 2 - mesh_spectral_getGrid(2) < 2_pInt .or. & ! -"- - mesh_spectral_getGrid(3) < 1_pInt) & ! must be at least 1 + if(any(mesh_spectral_getGrid < 1_pInt)) & call IO_error(error_ID = 843_pInt, ext_msg='mesh_spectral_getGrid') end function mesh_spectral_getGrid @@ -1116,7 +1108,7 @@ subroutine mesh_spectral_count_cpSizes implicit none integer(pInt) :: t - t = FE_geomtype(FE_mapElemtype('C3D8R')) ! fake 3D hexahedral 8 node 1 IP element + t = FE_geomtype(FE_mapElemtype('C3D8R')) ! fake 3D hexahedral 8 node 1 IP element mesh_maxNnodes = FE_Nnodes(t) mesh_maxNips = FE_Nips(t) @@ -1131,8 +1123,8 @@ end subroutine mesh_spectral_count_cpSizes !! Allocates global arrays 'mesh_node0' and 'mesh_node' !-------------------------------------------------------------------------------------------------- subroutine mesh_spectral_build_nodes() - - use numerics, only: numerics_unitlength + use numerics, & + only: numerics_unitlength implicit none integer(pInt) :: n @@ -1152,7 +1144,7 @@ subroutine mesh_spectral_build_nodes() / real(res(3),pReal) end forall - mesh_node = mesh_node0 !why? + mesh_node = mesh_node0 ! why? end subroutine mesh_spectral_build_nodes @@ -1212,8 +1204,8 @@ subroutine mesh_spectral_build_elements(myUnit) read(myUnit,'(a65536)') line enddo - allocate (mesh_element (4_pInt+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt - allocate (microstructures (1_pInt+maxIntCount)) ; microstructures = 2_pInt + allocate (mesh_element (4_pInt+mesh_maxNnodes,mesh_NcpElems)); mesh_element = 0_pInt + allocate (microstructures (1_pInt+maxIntCount)); microstructures = 2_pInt e = 0_pInt do while (e < mesh_NcpElems .and. microstructures(1) > 0_pInt) ! fill expected number of elements, stop at end of data (or blank line!) @@ -1239,7 +1231,7 @@ subroutine mesh_spectral_build_elements(myUnit) enddo deallocate(microstructures) - if (e /= mesh_NcpElems) call IO_error(880_pInt,e) + if (e /= mesh_NcpElems) call IO_error(880_pInt,e) !@ToDo does that make sense? end subroutine mesh_spectral_build_elements @@ -1857,8 +1849,8 @@ function mesh_deformedCoordsLinear(gDim,F,FavgIn) result(coords) ! report if (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) then write(6,'(a)') ' Restore geometry using linear integration' - write(6,'(a,3(e12.5))') ' Dimension: ', gDim - write(6,'(a,3(i5))') ' Resolution:', iRes + write(6,'(a,3(i12 ))') ' grid a b c: ', iRes + write(6,'(a,3(f12.5))') ' size x y z: ', gDim endif !-------------------------------------------------------------------------------------------------- @@ -1961,20 +1953,20 @@ function mesh_deformedCoordsFFT(gDim,F,FavgIn,scalingIn) result(coords) real(pReal), intent(in), dimension(3) :: gDim real(pReal), intent(in), dimension(3,3), optional :: FavgIn real(pReal), intent(in), dimension(3), optional :: scalingIn -! function ! allocatable arrays for fftw c routines - type(C_PTR) :: fftw_forth, fftw_back + type(C_PTR) :: planForth, planBack type(C_PTR) :: coords_fftw, defgrad_fftw real(pReal), dimension(:,:,:,:,:), pointer :: F_real complex(pReal), dimension(:,:,:,:,:), pointer :: F_fourier real(pReal), dimension(:,:,:,:), pointer :: coords_real complex(pReal), dimension(:,:,:,:), pointer :: coords_fourier ! other variables - integer(pInt) :: i, j, k, m, res1_red + integer(pInt) :: i, j, k, m, res1Red integer(pInt), dimension(3) :: k_s, iRes real(pReal), dimension(3) :: scaling, step, offset_coords, integrator real(pReal), dimension(3,3) :: Favg + integer(pInt), dimension(2:3,2) :: Nyquist ! highest frequencies to be removed (1 if even, 2 if odd) if (present(scalingIn)) then where (scalingIn < 0.0_pReal) ! the f2py way to tell it is not present @@ -1988,53 +1980,56 @@ function mesh_deformedCoordsFFT(gDim,F,FavgIn,scalingIn) result(coords) iRes = [size(F,3),size(F,4),size(F,5)] integrator = gDim / 2.0_pReal / PI ! see notes where it is used - res1_red = iRes(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c) + res1Red = iRes(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c) step = gDim/real(iRes, pReal) + Nyquist(2,1:2) = [res(2)/2_pInt + 1_pInt, res(2)/2_pInt + 1_pInt + mod(res(2),2_pInt)] + Nyquist(3,1:2) = [res(3)/2_pInt + 1_pInt, res(3)/2_pInt + 1_pInt + mod(res(3),2_pInt)] !-------------------------------------------------------------------------------------------------- ! report if (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) then write(6,'(a)') ' Restore geometry using FFT-based integration' - write(6,'(a,3(e12.5))') ' Dimension: ', gDim - write(6,'(a,3(i5))') ' Resolution:', iRes + write(6,'(a,3(i12 ))') ' grid a b c: ', iRes + write(6,'(a,3(f12.5))') ' size x y z: ', gDim endif !-------------------------------------------------------------------------------------------------- -! sanity checks - - if ((mod(iRes(3),2_pInt)/=0_pInt .and. iRes(3) /= 1_pInt) .or. & - mod(iRes(2),2_pInt)/=0_pInt .or. & - mod(iRes(1),2_pInt)/=0_pInt) & - call IO_error(0_pInt,ext_msg='Resolution in mesh_deformedCoordsFFT') +! sanity check if (pReal /= C_DOUBLE .or. pInt /= C_INT) & call IO_error(0_pInt,ext_msg='Fortran to C in mesh_deformedCoordsFFT') !-------------------------------------------------------------------------------------------------- ! allocation and FFTW initialization - call fftw_set_timelimit(fftw_timelimit) - defgrad_fftw = fftw_alloc_complex(int(res1_red *iRes(2)*iRes(3)*9_pInt,C_SIZE_T)) ! C_SIZE_T is of type integer(8) - call c_f_pointer(defgrad_fftw, F_real, [iRes(1)+2_pInt,iRes(2),iRes(3),3_pInt,3_pInt]) - call c_f_pointer(defgrad_fftw, F_fourier,[res1_red ,iRes(2),iRes(3),3_pInt,3_pInt]) - coords_fftw = fftw_alloc_complex(int(res1_red *iRes(2)*iRes(3)*3_pInt,C_SIZE_T)) ! C_SIZE_T is of type integer(8) - call c_f_pointer(coords_fftw, coords_real, [iRes(1)+2_pInt,iRes(2),iRes(3),3_pInt]) - call c_f_pointer(coords_fftw, coords_fourier, [res1_red ,iRes(2),iRes(3),3_pInt]) - fftw_forth = fftw_plan_many_dft_r2c(3_pInt,[iRes(3),iRes(2) ,iRes(1)],9_pInt,& ! dimensions , length in each dimension in reversed order - F_real,[iRes(3),iRes(2) ,iRes(1)+2_pInt],& ! input data , physical length in each dimension in reversed order - 1_pInt, iRes(3)*iRes(2)*(iRes(1)+2_pInt),& ! striding , product of physical lenght in the 3 dimensions - F_fourier,[iRes(3),iRes(2) ,res1_red],& - 1_pInt, iRes(3)*iRes(2)* res1_red,fftw_planner_flag) + defgrad_fftw = fftw_alloc_complex(int(res1Red *iRes(2)*iRes(3)*9_pInt,C_SIZE_T)) ! C_SIZE_T is of type integer(8) + coords_fftw = fftw_alloc_complex(int(res1Red *iRes(2)*iRes(3)*3_pInt,C_SIZE_T)) ! C_SIZE_T is of type integer(8) + call c_f_pointer(defgrad_fftw, F_real, & + [iRes(1)+2_pInt-mod(iRes(1),2_pInt),iRes(2),iRes(3),3_pInt,3_pInt]) + call c_f_pointer(defgrad_fftw, F_fourier, & + [res1Red, iRes(2),iRes(3),3_pInt,3_pInt]) + call c_f_pointer(coords_fftw, coords_real, & + [iRes(1)+2_pInt-mod(iRes(1),2_pInt),iRes(2),iRes(3),3_pInt]) + call c_f_pointer(coords_fftw, coords_fourier, & + [res1Red, iRes(2),iRes(3),3_pInt]) - fftw_back = fftw_plan_many_dft_c2r(3_pInt,[iRes(3),iRes(2) ,iRes(1)],3_pInt,& - coords_fourier,[iRes(3),iRes(2) ,res1_red],& - 1_pInt, iRes(3)*iRes(2)* res1_red,& - coords_real,[iRes(3),iRes(2) ,iRes(1)+2_pInt],& - 1_pInt, iRes(3)*iRes(2)*(iRes(1)+2_pInt),fftw_planner_flag) - forall(k = 1_pInt:iRes(3), j = 1_pInt:iRes(2), i = 1_pInt:iRes(1)) & - F_real(i,j,k,1:3,1:3) = F(1:3,1:3,i,j,k) ! F_real is overwritten during plan creation and is larger (padding) + call fftw_set_timelimit(fftw_timelimit) + planForth = fftw_plan_many_dft_r2c(3_pInt,[iRes(3),iRes(2) ,iRes(1)],9_pInt,& ! dimensions , length in each dimension in reversed order + F_real,[iRes(3),iRes(2) ,iRes(1)+2_pInt-mod(iRes(1),2_pInt)],& ! input data , physical length in each dimension in reversed order + 1_pInt, iRes(3)*iRes(2)*(iRes(1)+2_pInt-mod(iRes(1),2_pInt)),& ! striding , product of physical lenght in the 3 dimensions + F_fourier,[iRes(3),iRes(2) ,res1Red],& + 1_pInt, iRes(3)*iRes(2)* res1Red,fftw_planner_flag) + + planBack = fftw_plan_many_dft_c2r(3_pInt,[iRes(3),iRes(2) ,iRes(1)],3_pInt,& + coords_fourier,[iRes(3),iRes(2) ,res1Red],& + 1_pInt, iRes(3)*iRes(2)* res1Red,& + coords_real,[iRes(3),iRes(2) ,iRes(1)+2_pInt-mod(iRes(1),2_pInt)],& + 1_pInt, iRes(3)*iRes(2)*(iRes(1)+2_pInt-mod(iRes(1),2_pInt)),& + fftw_planner_flag) + F_real(1:iRes(1),1:iRes(2),1:iRes(3),1:3,1:3) = & + reshape(F,[res(1),res(2),res(3),3,3], order = [4,5,1,2,3]) ! F_real is overwritten during plan creatio, is larger (padding) and has different order !-------------------------------------------------------------------------------------------------- ! FFT - call fftw_execute_dft_r2c(fftw_forth, F_real, F_fourier) + call fftw_execute_dft_r2c(planForth, F_real, F_fourier) !-------------------------------------------------------------------------------------------------- ! if no average F is given, compute it in Fourier space @@ -2050,13 +2045,16 @@ function mesh_deformedCoordsFFT(gDim,F,FavgIn,scalingIn) result(coords) !-------------------------------------------------------------------------------------------------- ! remove highest frequency in each direction, in third direction only if not 2D - F_fourier( iRes(1)/2_pInt+1_pInt,1:iRes(2) ,1:iRes(3) ,1:3,1:3) & - = cmplx(0.0_pReal,0.0_pReal,pReal) - F_fourier( 1:res1_red ,iRes(2)/2_pInt+1_pInt,1:iRes(3) ,1:3,1:3) & - = cmplx(0.0_pReal,0.0_pReal,pReal) - if(iRes(3)>1_pInt) & - F_fourier(1:res1_red ,1:iRes(2) ,iRes(3)/2_pInt+1_pInt,1:3,1:3) & - = cmplx(0.0_pReal,0.0_pReal,pReal) + + if(iRes(1)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation + F_fourier (res1Red, 1:iRes(2), 1:iRes(3), 1:3,1:3) & + = cmplx(0.0_pReal,0.0_pReal,pReal) + if(iRes(2)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation + F_fourier (1:res1Red,Nyquist(2,1):Nyquist(2,2),1:iRes(3), 1:3,1:3) & + = cmplx(0.0_pReal,0.0_pReal,pReal) + if(iRes(3)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation + F_fourier (1:res1Red,1:iRes(2), Nyquist(3,1):Nyquist(3,2),1:3,1:3) & + = cmplx(0.0_pReal,0.0_pReal,pReal) !-------------------------------------------------------------------------------------------------- ! integration in Fourier space @@ -2067,7 +2065,7 @@ function mesh_deformedCoordsFFT(gDim,F,FavgIn,scalingIn) result(coords) do j = 1_pInt, iRes(2) k_s(2) = j-1_pInt if(j > iRes(2)/2_pInt+1_pInt) k_s(2) = k_s(2)-iRes(2) - do i = 1_pInt, res1_red + do i = 1_pInt, res1Red k_s(1) = i-1_pInt do m = 1_pInt,3_pInt coords_fourier(i,j,k,m) = sum(F_fourier(i,j,k,m,1:3)*& @@ -2079,12 +2077,11 @@ function mesh_deformedCoordsFFT(gDim,F,FavgIn,scalingIn) result(coords) !-------------------------------------------------------------------------------------------------- ! iFFT and freeing memory - call fftw_execute_dft_c2r(fftw_back,coords_fourier,coords_real) - coords_real = coords_real/real(product(iRes),pReal) - forall(k = 1_pInt:iRes(3), j = 1_pInt:iRes(2), i = 1_pInt:iRes(1)) & - coords(1:3,i,j,k) = coords_real(i,j,k,1:3) - call fftw_destroy_plan(fftw_forth) - call fftw_destroy_plan(fftw_back) + call fftw_execute_dft_c2r(planBack,coords_fourier,coords_real) + coords = reshape(coords_real(1:iRes(1),1:iRes(2),1:iRes(3),1:3), [3,iRes(1),iRes(2),iRes(3)], & + order = [2,3,4,1])/real(product(iRes),pReal) ! weight and change order + call fftw_destroy_plan(planForth) + call fftw_destroy_plan(planBack) call fftw_free(defgrad_fftw) call fftw_free(coords_fftw) @@ -2137,8 +2134,8 @@ function mesh_volumeMismatch(gDim,F,nodes) result(vMismatch) ! report and check if (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) then write(6,'(a)') ' Calculating volume mismatch' - write(6,'(a,3(e12.5))') ' Dimension: ', gDim - write(6,'(a,3(i5))') ' Resolution:', iRes + write(6,'(a,3(i12 ))') ' grid a b c: ', iRes + write(6,'(a,3(f12.5))') ' size x y z: ', gDim endif if (any([iRes/=size(nodes,2)-1_pInt,iRes/=size(nodes,3)-1_pInt,iRes/=size(nodes,4)-1_pInt]))& @@ -2209,8 +2206,8 @@ function mesh_shapeMismatch(gDim,F,nodes,centres) result(sMismatch) ! report and check if (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) then write(6,'(a)') ' Calculating shape mismatch' - write(6,'(a,3(e12.5))') ' Dimension: ', gDim - write(6,'(a,3(i5))') ' Resolution:', iRes + write(6,'(a,3(i12 ))') ' grid a b c: ', iRes + write(6,'(a,3(f12.5))') ' size x y z: ', gDim endif if(any([iRes/=size(nodes,2)-1_pInt,iRes/=size(nodes,3)-1_pInt,iRes/=size(nodes,4)-1_pInt]) .or.& @@ -3553,8 +3550,8 @@ subroutine mesh_build_ipAreas enddo enddo - end select - enddo + end select + enddo end subroutine mesh_build_ipAreas @@ -4610,6 +4607,12 @@ subroutine mesh_build_FEdata ! indicates which subnodes make up the interfaces enclosing the IP volume. ! The sorting convention is such that the outward pointing normal ! follows from a right-handed traversal of the face node list. + ! For two-dimensional elements, which only have lines as "interface" + ! one nevertheless has to specify each interface by a closed path, + ! e.g., 1,2, 2,1, assuming the line connects nodes 1 and 2. + ! This will result in zero ipVolume and interfaceArea, but is not + ! detrimental at the moment since non-local constitutive laws are + ! currently not foreseen in 2D cases. me = 0_pInt me = me + 1_pInt diff --git a/code/numerics.f90 b/code/numerics.f90 index 66be4d3e0..701522ee4 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -107,8 +107,8 @@ module numerics itmin = 2_pInt, & !< minimum number of iterations maxCutBack = 3_pInt, & !< max number of cut backs regridMode = 0_pInt, & !< 0: no regrid; 1: regrid if DAMASK doesn't converge; 2: regrid if DAMASK or BVP Solver doesn't converge - divergence_correction = 2_pInt !< correct divergence calculation in fourier space 0: no correction, 1: dimension scaled to 1, 2: dimension scaled to Npoints, 3: dimension scaled to sqrt(Npoints) - logical, protected , public :: & + divergence_correction = 2_pInt !< correct divergence calculation in fourier space 0: no correction, 1: size scaled to 1, 2: size scaled to Npoints + logical, protected, public :: & memory_efficient = .true., & !< for fast execution (pre calculation of gamma_hat), Default .true.: do not precalculate update_gamma = .false. !< update gamma operator with current stiffness, Default .false.: use initial stiffness #endif @@ -484,7 +484,7 @@ subroutine numerics_init if (itmax <= 1_pInt) call IO_error(301_pInt,ext_msg='itmax') if (itmin > itmax .or. itmin < 1_pInt) call IO_error(301_pInt,ext_msg='itmin') if (divergence_correction < 0_pInt .or. & - divergence_correction > 3_pInt) call IO_error(301_pInt,ext_msg='divergence_correction') + divergence_correction > 2_pInt) call IO_error(301_pInt,ext_msg='divergence_correction') if (maxCutBack < 0_pInt) call IO_error(301_pInt,ext_msg='maxCutBack') if (update_gamma .and. & .not. memory_efficient) call IO_error(error_ID = 847_pInt)