From 724380c624354beb1170327e4ba7017e9471777c Mon Sep 17 00:00:00 2001 From: Pratheek Shanthraj Date: Tue, 12 Nov 2013 00:44:23 +0000 Subject: [PATCH] =?UTF-8?q?Coupled=20phase=20field=20implementation=20for?= =?UTF-8?q?=20the=20basic=20petsc=20solver.=20needs=20clearing=20up?= =?UTF-8?q?=E2=80=A6=20work=20in=20progress?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- code/DAMASK_spectral_driver.f90 | 1 - code/DAMASK_spectral_solverBasicPETSc.f90 | 126 +++++++++++++++++----- code/DAMASK_spectral_utilities.f90 | 117 +++++++++++++++++++- code/crystallite.f90 | 5 +- 4 files changed, 217 insertions(+), 32 deletions(-) diff --git a/code/DAMASK_spectral_driver.f90 b/code/DAMASK_spectral_driver.f90 index ce7745dbd..668da55fd 100644 --- a/code/DAMASK_spectral_driver.f90 +++ b/code/DAMASK_spectral_driver.f90 @@ -523,7 +523,6 @@ program DAMASK_spectral_Driver ' increment ', totalIncsCounter, ' NOT converged' notConvergedCounter = notConvergedCounter + 1_pInt endif; flush(6) - call utilities_temperatureUpdate() if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency write(6,'(1/,a)') ' ... writing results to file ......................................' write(resUnit) materialpoint_results ! write result to file diff --git a/code/DAMASK_spectral_solverBasicPETSc.f90 b/code/DAMASK_spectral_solverBasicPETSc.f90 index 0d821f20a..2af41b462 100644 --- a/code/DAMASK_spectral_solverBasicPETSc.f90 +++ b/code/DAMASK_spectral_solverBasicPETSc.f90 @@ -62,6 +62,10 @@ module DAMASK_spectral_SolverBasicPETSc !-------------------------------------------------------------------------------------------------- ! common pointwise data real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, Fdot, F_lastInc2 + real(pReal), private, dimension(:,:,:), allocatable :: & + heatSource_lastInc, & + phaseField_lastInc, & + phaseFieldDot complex(pReal), private, dimension(:,:,:,:,:), allocatable :: inertiaField_fourier !-------------------------------------------------------------------------------------------------- @@ -78,7 +82,7 @@ module DAMASK_spectral_SolverBasicPETSc C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 stiffness S = 0.0_pReal !< current compliance (filled up with zeros) - real(pReal), private :: err_stress, err_div, err_divPrev, err_divDummy + real(pReal), private :: err_stress, err_div, err_divPrev, err_divDummy, err_phaseField, phaseField_Avg logical, private :: ForwardData integer(pInt), private :: & totalIter = 0_pInt !< total iteration in current increment @@ -144,8 +148,10 @@ subroutine basicPETSc_init(temperature) temperature #include #include +#include real(pReal), dimension(:,:,:,:,:), allocatable :: P - PetscScalar, dimension(:,:,:,:), pointer :: F + PetscScalar, dimension(:,:,:,:), pointer :: xx_psc, F + PetscScalar, dimension(:,:,:), pointer :: phaseField PetscErrorCode :: ierr PetscObject :: dummy real(pReal), dimension(3,3) :: & @@ -167,6 +173,9 @@ subroutine basicPETSc_init(temperature) allocate (F_lastInc2(3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal) allocate (Fdot (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal) allocate (inertiaField_fourier (grid1Red,grid(2),grid(3),3,3),source = cmplx(0.0_pReal,0.0_pReal,pReal)) + allocate (heatSource_lastInc (grid(1),grid(2),grid(3)),source = 0.0_pReal) + allocate (phaseField_lastInc(grid(1),grid(2),grid(3)),source = 0.0_pReal) + allocate (phaseFieldDot (grid(1),grid(2),grid(3)),source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc @@ -174,9 +183,10 @@ subroutine basicPETSc_init(temperature) call DMDACreate3d(PETSC_COMM_WORLD, & DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, & DMDA_STENCIL_BOX,grid(1),grid(2),grid(3),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, & - 9,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr) + 10,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr) CHKERRQ(ierr) call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) + !call DMDASNESSetFunctionLocal(da,INSERT_VALUES,BasicPETSC_formResidual,dummy,ierr); CHKERRQ(ierr) ! needed for newer versions of petsc call DMDASetLocalFunction(da,BasicPETSC_formResidual,ierr); CHKERRQ(ierr) call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) call SNESSetConvergenceTest(snes,BasicPETSC_converged,dummy,PETSC_NULL_FUNCTION,ierr) @@ -188,12 +198,16 @@ subroutine basicPETSc_init(temperature) !-------------------------------------------------------------------------------------------------- ! init fields - call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with - + call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with + F => xx_psc(0:8,:,:,:) + phaseField => xx_psc(9,:,:,:) + if (restartInc == 1_pInt) then ! no deformation (no restart) F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid(3)) ! initialize to identity F = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)]) F_lastInc2 = F_lastInc + phaseField = temperature + phaseField_lastInc = temperature elseif (restartInc > 1_pInt) then ! using old values from file if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & @@ -227,10 +241,10 @@ subroutine basicPETSc_init(temperature) mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(& F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)]) call Utilities_constitutiveResponse(& - reshape(F(0:8,0:grid(1)-1_pInt,0:grid(2)-1_pInt,0:grid(3)-1_pInt),[3,3,grid(1),grid(2),grid(3)]),& - reshape(F(0:8,0:grid(1)-1_pInt,0:grid(2)-1_pInt,0:grid(3)-1_pInt),[3,3,grid(1),grid(2),grid(3)]),& + reshape(F,[3,3,grid(1),grid(2),grid(3)]),& + reshape(F,[3,3,grid(1),grid(2),grid(3)]),& temperature,0.0_pReal,P,C_volAvg,C_minmaxAvg,temp33_Real,.false.,math_I3) - call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! write data back into PETSc + call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! write data back into PETSc if (restartInc == 1_pInt) then ! use initial stiffness as reference stiffness temp3333_Real = C_minMaxAvg endif @@ -267,6 +281,8 @@ type(tSolutionState) function basicPETSc_solution( & use FEsolving, only: & restartWrite, & terminallyIll + use homogenization, only: & + materialpoint_heat implicit none #include @@ -291,12 +307,14 @@ type(tSolutionState) function basicPETSc_solution( & !-------------------------------------------------------------------------------------------------- ! PETSc Data - PetscScalar, pointer :: F(:,:,:,:) + PetscScalar, pointer :: xx_psc(:,:,:,:), F(:,:,:,:), phaseField(:,:,:) PetscErrorCode :: ierr SNESConvergedReason :: reason incInfo = incInfoIn - call DMDAVecGetArrayF90(da,solution_vec,F,ierr) + call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with + F => xx_psc(0:8,:,:,:) + phaseField => xx_psc(9,:,:,:) !-------------------------------------------------------------------------------------------------- ! restart information for spectral solver if (restartWrite) then @@ -327,6 +345,7 @@ type(tSolutionState) function basicPETSc_solution( & F_aim = F_aim_lastInc F = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)]) C_volAvg = C_volAvgLastInc + phaseField = phaseField_lastInc else C_volAvgLastInc = C_volAvg @@ -349,14 +368,18 @@ type(tSolutionState) function basicPETSc_solution( & F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)]) Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,params%rotation_BC), & timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)])) + phaseFieldDot = (phaseField - phaseField_lastInc)/timeinc_old F_lastInc2 = F_lastInc F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid(3)]) + phaseField_lastInc = phaseField + heatSource_lastInc = reshape(materialpoint_heat(1,:),[grid(1),grid(2),grid(3)]) endif F_aim = F_aim + f_aimDot * timeinc F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot,math_rotate_backward33(F_aim, & rotation_BC)),[9,grid(1),grid(2),grid(3)]) - call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) + phaseField = phaseField_lastInc + phaseFieldDot*timeinc + call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) @@ -412,27 +435,51 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) wgt, & field_real, & field_fourier, & + phaseField_real, & + phaseField_fourier, & Utilities_FFTforward, & Utilities_FFTbackward, & + utilities_scalarFFTforward, & + utilities_scalarFFTbackward, & Utilities_fourierConvolution, & Utilities_inverseLaplace, & + Utilities_diffusion, & Utilities_constitutiveResponse, & Utilities_divergenceRMS use IO, only: & IO_intOut + use crystallite, only: & + crystallite_temperature + use homogenization, only: & + materialpoint_heat implicit none DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & in - PetscScalar, dimension(3,3,grid(1),grid(2),grid(3)) :: & - x_scal, & - f_scal + PetscScalar, target, dimension(10, & + XG_RANGE,YG_RANGE,ZG_RANGE) :: & + x_scal + PetscScalar, target, dimension(10, & + X_RANGE,Y_RANGE,Z_RANGE) :: & + f_scal + PetscScalar, pointer, dimension(:,:,:,:) :: & + F, & + residual_F + PetscScalar, pointer, dimension(:,:,:) :: & + phaseField, & + residual_phaseField PetscInt :: & PETScIter, & nfuncs PetscObject :: dummy PetscErrorCode :: ierr + real(pReal) :: mobility, diffusivity + F => x_scal(1:9,1:grid(1),1:grid(2),1:grid(3)) + residual_F => f_scal(1:9,1:grid(1),1:grid(2),1:grid(3)) + phaseField => x_scal(10,1:grid(1),1:grid(2),1:grid(3)) + residual_phaseField => f_scal(10,1:grid(1),1:grid(2),1:grid(3)) + call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) @@ -453,24 +500,26 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! evaluate inertia -if (params%density > 0.0_pReal) then - f_scal = ((x_scal - F_lastInc)/params%timeinc - (F_lastInc - F_lastInc2)/params%timeincOld)/& - ((params%timeinc + params%timeincOld)/2.0_pReal) - f_scal = params%density*product(geomSize/grid)*f_scal + dynamic: if (params%density > 0.0_pReal) then + residual_F = ((F - reshape(F_lastInc,[9,grid(1),grid(2),grid(3)]))/params%timeinc - & + reshape(F_lastInc - F_lastInc2, [9,grid(1),grid(2),grid(3)])/params%timeincOld)/& + ((params%timeinc + params%timeincOld)/2.0_pReal) + residual_F = params%density*product(geomSize/grid)*residual_F field_real = 0.0_pReal - field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3) = reshape(f_scal,[grid(1),grid(2),grid(3),3,3],& + field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3) = reshape(residual_F,[grid(1),grid(2),grid(3),3,3],& order=[4,5,1,2,3]) ! field real has a different order call Utilities_FFTforward() call Utilities_inverseLaplace() inertiaField_fourier = field_fourier - else + else dynamic inertiaField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal) - endif + endif dynamic !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - call Utilities_constitutiveResponse(F_lastInc,x_scal,params%temperature,params%timeinc, & - f_scal,C_volAvg,C_minmaxAvg,P_av,ForwardData,params%rotation_BC) + crystallite_temperature(1,1_pInt:product(grid)) = reshape(phaseField,[product(grid)]) + call Utilities_constitutiveResponse(F_lastInc,F,params%temperature,params%timeinc, & + residual_F,C_volAvg,C_minmaxAvg,P_av,ForwardData,params%rotation_BC) ForwardData = .false. !-------------------------------------------------------------------------------------------------- @@ -482,7 +531,7 @@ if (params%density > 0.0_pReal) then !-------------------------------------------------------------------------------------------------- ! updated deformation gradient using fix point algorithm of basic scheme field_real = 0.0_pReal - field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3) = reshape(f_scal,[grid(1),grid(2),grid(3),3,3],& + field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3) = reshape(residual_F,[grid(1),grid(2),grid(3),3,3],& order=[4,5,1,2,3]) ! field real has a different order call Utilities_FFTforward() field_fourier = field_fourier + inertiaField_fourier @@ -492,7 +541,29 @@ if (params%density > 0.0_pReal) then !-------------------------------------------------------------------------------------------------- ! constructing residual - f_scal = reshape(field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3),shape(x_scal),order=[3,4,5,1,2]) + residual_F = reshape(field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3),& + [9,grid(1),grid(2),grid(3)],order=[2,3,4,1]) + +!-------------------------------------------------------------------------------------------------- +! constructing phase field residual + diffusivity = 400.0; mobility = 4e6 !sample coefficients for copper... need to clear this up with better mode of input + phaseField_real = 0.0_pReal + phaseField_real(1:grid(1),1:grid(2),1:grid(3)) = phaseField_lastInc + call utilities_scalarFFTforward() + call utilities_diffusion(diffusivity/mobility,params%timeinc) + call utilities_scalarFFTbackward() + residual_phaseField = phaseField_real(1:grid(1),1:grid(2),1:grid(3)) + phaseField_real = 0.0_pReal + phaseField_real(1:grid(1),1:grid(2),1:grid(3)) = & + (heatSource_lastInc + reshape(materialpoint_heat(1,1_pInt:product(grid)),[grid(1),grid(2),grid(3)]))* & + params%timeinc/mobility/2.0_pReal + call utilities_scalarFFTforward() + call utilities_diffusion(diffusivity/mobility,params%timeinc/2.0_pReal) + call utilities_scalarFFTbackward() + residual_phaseField = phaseField - & + (residual_phaseField + phaseField_real(1:grid(1),1:grid(2),1:grid(3))) + + err_phaseField = maxval(abs(residual_phaseField)); phaseField_Avg = sum(phaseField)*wgt end subroutine BasicPETSc_formResidual @@ -531,7 +602,8 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du converged: if ((totalIter >= itmin .and. & all([ err_div/divTol, & - err_stress/stressTol ] < 1.0_pReal)) & + err_stress/stressTol, & + err_phaseField/phaseField_Avg/1.0e-3 ] < 1.0_pReal)) & .or. terminallyIll) then reason = 1 elseif (totalIter >= itmax) then converged @@ -547,6 +619,8 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du err_div/divTol, ' (',err_div,' / m, tol =',divTol,')' write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', & err_stress/stressTol, ' (',err_stress, ' Pa, tol =',stressTol,')' + write(6,'(a,f10.2,a,es8.2,a,es9.2,a)') ' error phase field = ', & + err_phaseField/phaseField_Avg/1.0e-3, ' (',err_phaseField/phaseField_Avg, ' Pa, tol =',1.0e-3,')' write(6,'(/,a)') ' ===========================================================================' flush(6) diff --git a/code/DAMASK_spectral_utilities.f90 b/code/DAMASK_spectral_utilities.f90 index 186c787ce..31ecbc59a 100644 --- a/code/DAMASK_spectral_utilities.f90 +++ b/code/DAMASK_spectral_utilities.f90 @@ -47,6 +47,8 @@ module DAMASK_spectral_utilities integer(pInt), public :: grid1Red !< grid(1)/2 real(pReal), public, dimension(:,:,:,:,:), pointer :: field_real !< real representation (some stress or deformation) of field_fourier complex(pReal),public, dimension(:,:,:,:,:), pointer :: field_fourier !< field on which the Fourier transform operates + real(pReal), public, dimension(:,:,:), pointer :: phaseField_real !< real representation (some stress or deformation) of field_fourier + complex(pReal),public, dimension(:,:,:), pointer :: phaseField_fourier !< field on which the Fourier transform operates real(pReal), private, dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat !< gamma operator (field) for spectral method real(pReal), private, dimension(:,:,:,:), allocatable :: xi !< wave vector field for divergence and for gamma operator real(pReal), private, dimension(3,3,3,3) :: C_ref !< reference stiffness @@ -67,6 +69,8 @@ module DAMASK_spectral_utilities type(C_PTR), private :: & planForth, & !< FFTW plan P(x) to P(k) planBack, & !< FFTW plan F(k) to F(x) + planPhaseFieldForth, & !< FFTW plan P(x) to P(k) + planPhaseFieldBack, & !< 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 @@ -101,8 +105,11 @@ module DAMASK_spectral_utilities utilities_updateGamma, & utilities_FFTforward, & utilities_FFTbackward, & + utilities_scalarFFTforward, & + utilities_scalarFFTbackward, & utilities_fourierConvolution, & utilities_inverseLaplace, & + utilities_diffusion, & utilities_divergenceRMS, & utilities_curlRMS, & utilities_maskedCompliance, & @@ -169,6 +176,7 @@ subroutine utilities_init() integer(pInt), parameter :: fileUnit = 228_pInt integer(pInt), dimension(3) :: k_s type(C_PTR) :: & + phaseFieldFFT, & 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) @@ -231,8 +239,11 @@ subroutine utilities_init() ! allocation allocate (xi(3,grid1Red,grid(2),grid(3)),source = 0.0_pReal) ! frequencies, only half the size for first dimension tensorField = fftw_alloc_complex(int(grid1Red*grid(2)*grid(3)*9_pInt,C_SIZE_T)) ! allocate aligned data using a C function, C_SIZE_T is of type integer(8) + phaseFieldFFT = fftw_alloc_complex(int(grid1Red*grid(2)*grid(3),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, [grid(1)+2_pInt-mod(grid(1),2_pInt),grid(2),grid(3),3,3])! place a pointer for a real representation on tensorField call c_f_pointer(tensorField, field_fourier,[grid1Red, grid(2),grid(3),3,3])! place a pointer for a complex representation on tensorField + call c_f_pointer(phaseFieldFFT,phaseField_real,[grid(1)+2_pInt-mod(grid(1),2_pInt),grid(2),grid(3)])! place a pointer for a real representation on tensorField + call c_f_pointer(phaseFieldFFT,phaseField_fourier,[grid1Red,grid(2),grid(3)])! place a pointer for a complex representation on tensorField !-------------------------------------------------------------------------------------------------- ! general initialization of FFTW (see manual on fftw.org for more details) @@ -259,6 +270,21 @@ subroutine utilities_init() 1, grid(3)*grid(2)*(grid(1)+2_pInt-mod(grid(1),2_pInt)), & ! striding, product of physical length in the 3 dimensions fftw_planner_flag) ! planner precision +!-------------------------------------------------------------------------------------------------- +! creating plans for the convolution + planPhaseFieldForth = fftw_plan_many_dft_r2c(3,[grid(3),grid(2) ,grid(1)], 1, & ! dimensions, logical length in each dimension in reversed order, no. of transforms + phaseField_real,[grid(3),grid(2) ,grid(1)+2_pInt-mod(grid(1),2_pInt)], & ! input data, physical length in each dimension in reversed order + 1, grid(3)*grid(2)*(grid(1)+2_pInt-mod(grid(1),2_pInt)), & ! striding, product of physical length in the 3 dimensions + phaseField_fourier,[grid(3),grid(2) ,grid1Red], & ! output data, physical length in each dimension in reversed order + 1, grid(3)*grid(2)* grid1Red, fftw_planner_flag) ! striding, product of physical length in the 3 dimensions, planner precision + + planPhaseFieldBack = fftw_plan_many_dft_c2r(3,[grid(3),grid(2) ,grid(1)], 1, & ! dimensions, logical length in each dimension in reversed order, no. of transforms + phaseField_fourier,[grid(3),grid(2) ,grid1Red], & ! input data, physical length in each dimension in reversed order + 1, grid(3)*grid(2)* grid1Red, & ! striding, product of physical length in the 3 dimensions + phaseField_real,[grid(3),grid(2) ,grid(1)+2_pInt-mod(grid(1),2_pInt)], & ! output data, physical length in each dimension in reversed order + 1, grid(3)*grid(2)*(grid(1)+2_pInt-mod(grid(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 @@ -487,6 +513,59 @@ subroutine utilities_FFTbackward() end subroutine utilities_FFTbackward +!-------------------------------------------------------------------------------------------------- +!> @brief forward FFT of data in field_real to field_fourier with highest freqs. removed +!> @details Does an unweighted FFT transform from real to complex. +!> 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_scalarFFTforward() + use math + + implicit none + integer(pInt), dimension(2:3,2) :: Nyquist ! highest frequencies to be removed (1 if even, 2 if odd) + +!-------------------------------------------------------------------------------------------------- +! doing the FFT + call fftw_execute_dft_r2c(planPhaseFieldForth,phaseField_real,phaseField_fourier) + +!-------------------------------------------------------------------------------------------------- +! removing highest frequencies + Nyquist(2,1:2) = [grid(2)/2_pInt + 1_pInt, grid(2)/2_pInt + 1_pInt + mod(grid(2),2_pInt)] + Nyquist(3,1:2) = [grid(3)/2_pInt + 1_pInt, grid(3)/2_pInt + 1_pInt + mod(grid(3),2_pInt)] + + if(grid(1)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation + phaseField_fourier (grid1Red, 1:grid(2), 1:grid(3)) & + = cmplx(0.0_pReal,0.0_pReal,pReal) + if(grid(2)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation + phaseField_fourier (1:grid1Red,Nyquist(2,1):Nyquist(2,2),1:grid(3)) & + = cmplx(0.0_pReal,0.0_pReal,pReal) + if(grid(3)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation + phaseField_fourier (1:grid1Red,1:grid(2),Nyquist(3,1):Nyquist(3,2)) & + = cmplx(0.0_pReal,0.0_pReal,pReal) +end subroutine utilities_scalarFFTforward + + +!-------------------------------------------------------------------------------------------------- +!> @brief backward FFT of data in field_fourier to field_real +!> @details Does an inverse FFT transform from complex to real +!> 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 +!> results is weighted by number of points stored in wgt +!-------------------------------------------------------------------------------------------------- +subroutine utilities_scalarFFTbackward() + use math !< must use the whole module for use of FFTW + + implicit none + +!-------------------------------------------------------------------------------------------------- +! doing the iFFT + call fftw_execute_dft_c2r(planPhaseFieldBack,phaseField_fourier,phaseField_real) ! back transform of fluct deformation gradient + phaseField_real = phaseField_real * wgt ! normalize the result by number of elements + +end subroutine utilities_scalarFFTbackward + + !-------------------------------------------------------------------------------------------------- !> @brief doing convolution with inverse laplace kernel !-------------------------------------------------------------------------------------------------- @@ -520,6 +599,39 @@ enddo; enddo; enddo field_fourier(1,1,1,1:3,1:3) = cmplx(0.0_pReal,0.0_pReal,pReal) end subroutine utilities_inverseLaplace + + +!-------------------------------------------------------------------------------------------------- +!> @brief doing convolution with inverse laplace kernel +!-------------------------------------------------------------------------------------------------- +subroutine utilities_diffusion(coefficient,timeinc) + use math, only: & + PI + + implicit none + real(pReal),intent(in) :: timeinc, coefficient + integer(pInt) :: i, j, k + integer(pInt), dimension(3) :: k_s + + write(6,'(/,a)') ' ... doing diffusion .......................................................' + flush(6) + +!-------------------------------------------------------------------------------------------------- +! do the actual spectral method calculation (mechanical equilibrium) + do k = 1_pInt, grid(3) + k_s(3) = k - 1_pInt + if(k > grid(3)/2_pInt + 1_pInt) k_s(3) = k_s(3) - grid(3) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1 + do j = 1_pInt, grid(2) + k_s(2) = j - 1_pInt + if(j > grid(2)/2_pInt + 1_pInt) k_s(2) = k_s(2) - grid(2) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1 + do i = 1_pInt, grid1Red + k_s(1) = i - 1_pInt + phaseField_fourier(i,j,k) = phaseField_fourier(i,j,k)* & + cmplx(exp(-sum((2.0_pReal*PI*real(k_s,pReal)/geomSize)*(2.0_pReal*PI*real(k_s,pReal)/geomSize))* & + coefficient*timeinc),0.0_pReal,pReal) ! symmetry, junst running from 0,1,...,N/2,N/2+1 +enddo; enddo; enddo + +end subroutine utilities_diffusion !-------------------------------------------------------------------------------------------------- @@ -1007,19 +1119,20 @@ end function utilities_getFilter !-------------------------------------------------------------------------------------------------- !> @brief calculates filter for fourier convolution depending on type given in numerics.config !-------------------------------------------------------------------------------------------------- -subroutine utilities_temperatureUpdate() +subroutine utilities_temperatureUpdate(timeinc) use crystallite, only: & crystallite_temperature use homogenization, only: & materialpoint_heat implicit none + real(pReal),intent(in) :: timeinc integer :: & x,y,z,e real(pReal) :: & a forall(e=1_pInt:product(grid)) & - crystallite_temperature(1,e) = crystallite_temperature(1,e) + materialpoint_heat(1,e)*0.0001_pReal + crystallite_temperature(1,e) = crystallite_temperature(1,e) + materialpoint_heat(1,e)*timeinc e = 0_pInt z = 0_pInt y = 0_pInt diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 9f6ebfcd0..6fa5fc2b4 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -1303,9 +1303,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains crystallite_heat(g,i,e) = 0.98_pReal* & - math_mul33xx33(crystallite_P(1:3,1:3,g,i,e), & - crystallite_partionedF(1:3,1:3,g,i,e)-& - crystallite_partionedF0(1:3,1:3,g,i,e)) + abs(math_mul33xx33(math_Mandel6to33(crystallite_Tstar_v), & + crystallite_Lp(1:3,1:3,g,i,e))) enddo enddo enddo elementLooping12