diff --git a/code/DAMASK_spectral_driver.f90 b/code/DAMASK_spectral_driver.f90 index 91c3c0132..bfa9accd5 100644 --- a/code/DAMASK_spectral_driver.f90 +++ b/code/DAMASK_spectral_driver.f90 @@ -76,9 +76,7 @@ program DAMASK_spectral_Driver geomSize, & tBoundaryCondition, & tSolutionState, & - phaseFieldDataBin, & - cutBack, & - maxPhaseFields + cutBack use DAMASK_spectral_SolverBasic #ifdef PETSc use DAMASK_spectral_SolverBasicPETSC @@ -91,17 +89,14 @@ program DAMASK_spectral_Driver real(pReal), dimension (3,3) :: rotation = math_I3 !< rotation of BC type(tBoundaryCondition) :: P, & !< stress BC deformation !< deformation BC (Fdot or L) - type(phaseFieldDataBin) :: phaseFieldData(maxPhaseFields) real(pReal) :: time = 0.0_pReal, & !< length of increment - temperature = 300.0_pReal, & !< isothermal starting condition + temperature = 300.0_pReal, & !< isothermal starting conditions density = 0.0_pReal !< density integer(pInt) :: incs = 0_pInt, & !< number of increments outputfrequency = 1_pInt, & !< frequency of result writes restartfrequency = 0_pInt, & !< frequency of restart writes logscale = 0_pInt !< linear/logarithmic time inc flag - logical :: followFormerTrajectory = .true., & !< follow trajectory of former loadcase - thermal_active = .false., & !< activate thermal phase field - fracture_active = .false. !< activate fracture phase field + logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase end type tLoadCase !-------------------------------------------------------------------------------------------------- @@ -117,8 +112,7 @@ program DAMASK_spectral_Driver integer(pInt) :: & N_t = 0_pInt, & !< # of time indicators found in load case file N_n = 0_pInt, & !< # of increment specifiers found in load case file - N_def = 0_pInt, & !< # of rate of deformation specifiers found in load case file - nActivePhaseFields + N_def = 0_pInt !< # of rate of deformation specifiers found in load case file character(len=65536) :: & line @@ -236,22 +230,8 @@ program DAMASK_spectral_Driver loadCases(currentLoadCase)%P%values = math_plain9to33(temp_valueVector) case('t','time','delta') ! increment time loadCases(currentLoadCase)%time = IO_floatValue(line,positions,i+1_pInt) - case('temperature') + case('temp','temperature') ! starting temperature loadCases(currentLoadCase)%temperature = IO_floatValue(line,positions,i+1_pInt) - case('thermal') ! starting temperature, conductivity and mobility - loadCases(:)%phaseFieldData(1)%label = 'thermal' - loadCases(:)%phaseFieldData(1)%active = .true. - loadCases(:)%phaseFieldData(1)%phaseField0 = 300.0_pReal ! initialize to meaningful value if not defined - loadCases(currentLoadCase)%phaseFieldData(1)%phaseField0 = IO_floatValue(line,positions,i+1_pInt) - loadCases(currentLoadCase)%phaseFieldData(1)%diffusion = IO_floatValue(line,positions,i+2_pInt) - loadCases(currentLoadCase)%phaseFieldData(1)%mobility = IO_floatValue(line,positions,i+3_pInt) - case('fracture') ! starting damage, diffusion and mobility - loadCases(:)%phaseFieldData(2)%label = 'fracture' - loadCases(:)%phaseFieldData(2)%active = .true. - loadCases(:)%phaseFieldData(2)%phaseField0 = 1.0_pReal ! initialize to meaningful value if not defined - loadCases(currentLoadCase)%phaseFieldData(2)%phaseField0 = IO_floatValue(line,positions,i+1_pInt) - loadCases(currentLoadCase)%phaseFieldData(2)%diffusion = IO_floatValue(line,positions,i+2_pInt) - loadCases(currentLoadCase)%phaseFieldData(2)%mobility = IO_floatValue(line,positions,i+3_pInt) case('den','density') ! starting density loadCases(currentLoadCase)%density = IO_floatValue(line,positions,i+1_pInt) case('n','incs','increments','steps') ! number of increments @@ -290,15 +270,7 @@ program DAMASK_spectral_Driver loadCases(currentLoadCase)%rotation = math_plain9to33(temp_valueVector) end select enddo; enddo - close(FILEUNIT) - ! reorder phase field data to remove redundant non-active fields - nActivePhaseFields = 0_pInt - do i = 1, maxPhaseFields - if (loadCases(1)%phaseFieldData(i)%active) then - nActivePhaseFields = nActivePhaseFields + 1_pInt - loadCases(:)%phaseFieldData(nActivePhaseFields) = loadCases(:)%phaseFieldData(i) - endif - enddo + close(FILEUNIT) !-------------------------------------------------------------------------------------------------- ! consistency checks and output of load case @@ -364,7 +336,7 @@ program DAMASK_spectral_Driver call basic_init(loadCases(1)%temperature) #ifdef PETSc case (DAMASK_spectral_SolverBasicPETSc_label) - call basicPETSc_init(loadCases(1)%temperature,nActivePhaseFields,loadCases(1)%phaseFieldData(1:nActivePhaseFields)) + call basicPETSc_init(loadCases(1)%temperature) case (DAMASK_spectral_SolverAL_label) if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) & call IO_warning(42_pInt, ext_msg='debug Divergence') @@ -495,9 +467,7 @@ program DAMASK_spectral_Driver F_BC = loadCases(currentLoadCase)%deformation, & temperature_bc = loadCases(currentLoadCase)%temperature, & rotation_BC = loadCases(currentLoadCase)%rotation, & - density = loadCases(currentLoadCase)%density, & - nActivePhaseFields = nActivePhaseFields, & - phaseFieldData = loadCases(1)%phaseFieldData(1:nActivePhaseFields)) + density = loadCases(currentLoadCase)%density) case (DAMASK_spectral_SolverAL_label) solres = AL_solution (& incInfo,guess,timeinc,timeIncOld,remainingLoadCaseTime, & diff --git a/code/DAMASK_spectral_solverBasicPETSc.f90 b/code/DAMASK_spectral_solverBasicPETSc.f90 index 1b7bf6d9f..035b917d0 100644 --- a/code/DAMASK_spectral_solverBasicPETSc.f90 +++ b/code/DAMASK_spectral_solverBasicPETSc.f90 @@ -31,9 +31,7 @@ module DAMASK_spectral_SolverBasicPETSc use math, only: & math_I3 use DAMASK_spectral_Utilities, only: & - tSolutionState, & - phaseFieldDataBin, & - maxPhaseFields + tSolutionState implicit none private @@ -51,8 +49,6 @@ module DAMASK_spectral_SolverBasicPETSc real(pReal) :: timeincOld real(pReal) :: temperature real(pReal) :: density - integer(pInt) :: nActivePhaseFields - type(phaseFieldDataBin) :: phaseFieldData(maxPhaseFields) end type tSolutionParams type(tSolutionParams), private :: params @@ -66,11 +62,6 @@ module DAMASK_spectral_SolverBasicPETSc !-------------------------------------------------------------------------------------------------- ! common pointwise data real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, Fdot, F_lastInc2 - real(pReal), private, dimension(:,:,:,:), allocatable :: & - phaseFieldRHS_lastInc, & - phaseField_lastInc, & - phaseFieldRHS, & - phaseFieldDot complex(pReal), private, dimension(:,:,:,:,:), allocatable :: inertiaField_fourier !-------------------------------------------------------------------------------------------------- @@ -88,7 +79,6 @@ module DAMASK_spectral_SolverBasicPETSc 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, dimension(:), allocatable :: err_phaseField, phaseField_Avg logical, private :: ForwardData integer(pInt), private :: & totalIter = 0_pInt !< total iteration in current increment @@ -121,7 +111,7 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all neccessary fields and fills them with data, potentially from restart info !-------------------------------------------------------------------------------------------------- -subroutine basicPETSc_init(temperature,nActivePhaseFields,phaseFieldData) +subroutine basicPETSc_init(temperature) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) use IO, only: & IO_intOut, & @@ -150,14 +140,13 @@ subroutine basicPETSc_init(temperature,nActivePhaseFields,phaseFieldData) math_invSym3333 implicit none - integer(pInt), intent(in) :: nActivePhaseFields - type(phaseFieldDataBin), intent(in) :: phaseFieldData(nActivePhaseFields) - real(pReal), intent(inOut) :: temperature + real(pReal), intent(inout) :: & + temperature #include #include #include real(pReal), dimension(:,:,:,:,:), allocatable :: P - PetscScalar, dimension(:,:,:,:), pointer :: xx_psc, F + PetscScalar, dimension(:,:,:,:), pointer :: F PetscErrorCode :: ierr PetscObject :: dummy real(pReal), dimension(3,3) :: & @@ -165,27 +154,20 @@ subroutine basicPETSc_init(temperature,nActivePhaseFields,phaseFieldData) real(pReal), dimension(3,3,3,3) :: & temp3333_Real = 0.0_pReal KSP :: ksp - integer(pInt) :: i call Utilities_init() write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasicPETSc init -+>>>' write(6,'(a)') ' $Id$' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - + + allocate (P (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! allocate global fields - allocate (P (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal) allocate (F_lastInc (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal) 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 (phaseFieldRHS_lastInc (nActivePhaseFields,grid(1),grid(2),grid(3)),source = 0.0_pReal) - allocate (phaseField_lastInc (nActivePhaseFields,grid(1),grid(2),grid(3)),source = 0.0_pReal) - allocate (phaseFieldDot (nActivePhaseFields,grid(1),grid(2),grid(3)),source = 0.0_pReal) - allocate (phaseFieldRHS (nActivePhaseFields,grid(1),grid(2),grid(3)),source = 0.0_pReal) - allocate (err_phaseField(nActivePhaseFields), source = 0.0_pReal) - allocate (phaseField_Avg(nActivePhaseFields), source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc @@ -193,7 +175,7 @@ subroutine basicPETSc_init(temperature,nActivePhaseFields,phaseFieldData) 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+nActivePhaseFields,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr) + 9,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) @@ -208,16 +190,12 @@ subroutine basicPETSc_init(temperature,nActivePhaseFields,phaseFieldData) !-------------------------------------------------------------------------------------------------- ! init fields - call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with - F => xx_psc(0:8,:,:,:) - if (restartInc == 1_pInt) then ! no deformation (no restart) + call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with + + 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 - xx_psc(0:8,:,:,:) = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)]) + F = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)]) F_lastInc2 = F_lastInc - do i = 1, nActivePhaseFields - xx_psc(8+i,:,:,:) = phaseFieldData(i)%phaseField0 - phaseField_lastInc(i,:,:,:) = phaseFieldData(i)%phaseField0 - enddo 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)') & @@ -251,10 +229,10 @@ subroutine basicPETSc_init(temperature,nActivePhaseFields,phaseFieldData) 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,[3,3,grid(1),grid(2),grid(3)]),& - reshape(F,[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(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)]),& temperature,0.0_pReal,P,C_volAvg,C_minmaxAvg,temp33_Real,.false.,math_I3) - call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! write data back into PETSc + call DMDAVecRestoreArrayF90(da,solution_vec,F,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,8 +245,7 @@ end subroutine basicPETSc_init !> @brief solution for the Basic PETSC scheme with internal iterations !-------------------------------------------------------------------------------------------------- type(tSolutionState) function basicPETSc_solution( & - incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,temperature_bc,rotation_BC,density, & - nActivePhaseFields,phaseFieldData) + incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,temperature_bc,rotation_BC,density) use numerics, only: & update_gamma, & itmax @@ -292,8 +269,6 @@ type(tSolutionState) function basicPETSc_solution( & use FEsolving, only: & restartWrite, & terminallyIll - use homogenization, only: & - materialpoint_heat implicit none #include @@ -301,8 +276,6 @@ type(tSolutionState) function basicPETSc_solution( & !-------------------------------------------------------------------------------------------------- ! input data for solution - integer(pInt), intent(in) :: nActivePhaseFields - type(phaseFieldDataBin), intent(in) :: phaseFieldData(nActivePhaseFields) real(pReal), intent(in) :: & timeinc, & !< increment in time for current solution timeinc_old, & !< increment in time of last increment @@ -317,17 +290,15 @@ type(tSolutionState) function basicPETSc_solution( & character(len=*), intent(in) :: & incInfoIn real(pReal), dimension(3,3), intent(in) :: rotation_BC - integer(pInt) :: i !-------------------------------------------------------------------------------------------------- ! PETSc Data - PetscScalar, pointer :: xx_psc(:,:,:,:), F(:,:,:,:) + PetscScalar, pointer :: F(:,:,:,:) PetscErrorCode :: ierr SNESConvergedReason :: reason incInfo = incInfoIn - call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with - F => xx_psc(0:8,:,:,:) + call DMDAVecGetArrayF90(da,solution_vec,F,ierr) !-------------------------------------------------------------------------------------------------- ! restart information for spectral solver if (restartWrite) then @@ -356,11 +327,8 @@ type(tSolutionState) function basicPETSc_solution( & F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)]) if (cutBack) then F_aim = F_aim_lastInc - xx_psc(0:8,:,:,:) = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)]) + F = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)]) C_volAvg = C_volAvgLastInc - do i = 1, nActivePhaseFields - xx_psc(8+i,:,:,:) = phaseField_lastInc(i,:,:,:) - enddo else C_volAvgLastInc = C_volAvg @@ -383,22 +351,14 @@ 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)])) - do i = 1, nActivePhaseFields - phaseFieldDot(i,:,:,:) = (xx_psc(8+i,:,:,:) - phaseField_lastInc(i,:,:,:))/timeinc_old - phaseField_lastInc(i,:,:,:) = xx_psc(8+i,:,:,:) - phaseFieldRHS_lastInc(i,:,:,:) = phaseFieldRHS(i,:,:,:) - enddo F_lastInc2 = F_lastInc - F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid(3)]) + F_lastInc = reshape(F,[3,3,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)]) - do i = 1, nActivePhaseFields - xx_psc(8+i,:,:,:) = phaseField_lastInc(i,:,:,:) + phaseFieldDot(i,:,:,:)*timeinc - enddo - call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) + call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) @@ -416,8 +376,6 @@ type(tSolutionState) function basicPETSc_solution( & params%timeincOld = timeinc_old params%temperature = temperature_BC params%density = density - params%nActivePhaseFields = nActivePhaseFields - params%phaseFieldData(1:nActivePhaseFields) = phaseFieldData(1:nActivePhaseFields) call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr); CHKERRQ(ierr) call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr) @@ -456,47 +414,27 @@ 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, & - materialpoint_P implicit none DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & in - PetscScalar, target, dimension(9+params%nActivePhaseFields, & - XG_RANGE,YG_RANGE,ZG_RANGE) :: & - x_scal - PetscScalar, target, dimension(9+params%nActivePhaseFields, & - X_RANGE,Y_RANGE,Z_RANGE) :: & - f_scal - PetscScalar, pointer, dimension(:,:,:,:) :: & - F, & - residual_F + PetscScalar, dimension(3,3,grid(1),grid(2),grid(3)) :: & + x_scal, & + f_scal PetscInt :: & PETScIter, & nfuncs PetscObject :: dummy PetscErrorCode :: ierr - integer(pInt) :: i - 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)) - call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) @@ -517,37 +455,25 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! evaluate inertia - 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 +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 field_real = 0.0_pReal - 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],& + 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],& order=[4,5,1,2,3]) ! field real has a different order call Utilities_FFTforward() call Utilities_inverseLaplace() inertiaField_fourier = field_fourier - else dynamic + else inertiaField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal) - endif dynamic + endif !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - do i = 1, params%nActivePhaseFields - if(params%phaseFieldData(i)%label == 'thermal') & - crystallite_temperature(1,1_pInt:product(grid)) = & - reshape(x_scal(9+i,1:grid(1),1:grid(2),1:grid(3)),[product(grid)]) - enddo - - call Utilities_constitutiveResponse(F_lastInc,F,params%temperature,params%timeinc, & - residual_F,C_volAvg,C_minmaxAvg,P_av,ForwardData,params%rotation_BC) + call Utilities_constitutiveResponse(F_lastInc,x_scal,params%temperature,params%timeinc, & + f_scal,C_volAvg,C_minmaxAvg,P_av,ForwardData,params%rotation_BC) ForwardData = .false. - - do i = 1, params%nActivePhaseFields - if(params%phaseFieldData(i)%label == 'fracture') & - residual_F = residual_F * spread(x_scal(9+i,1:grid(1),1:grid(2),1:grid(3)),dim=1,ncopies=9) - enddo !-------------------------------------------------------------------------------------------------- ! stress BC handling @@ -558,7 +484,7 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! 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(residual_F,[grid(1),grid(2),grid(3),3,3],& + 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],& order=[4,5,1,2,3]) ! field real has a different order call Utilities_FFTforward() field_fourier = field_fourier + inertiaField_fourier @@ -566,79 +492,10 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) call Utilities_fourierConvolution(math_rotate_backward33(F_aim_lastIter-F_aim,params%rotation_BC)) call Utilities_FFTbackward() -!-------------------------------------------------------------------------------------------------- -! constructing phase field residual - do i = 1, params%nActivePhaseFields - select case (params%phaseFieldData(i)%label) - case ('thermal') - phaseField_real = 0.0_pReal - phaseField_real(1:grid(1),1:grid(2),1:grid(3)) = & - phaseField_lastInc(i,1:grid(1),1:grid(2),1:grid(3)) - call utilities_scalarFFTforward() - call utilities_diffusion(params%phaseFieldData(i)%diffusion,params%timeinc) - call utilities_scalarFFTbackward() - f_scal(9+i,1:grid(1),1:grid(2),1:grid(3)) = & - phaseField_real(1:grid(1),1:grid(2),1:grid(3)) - - phaseFieldRHS(i,1:grid(1),1:grid(2),1:grid(3)) = & - reshape(materialpoint_heat(1,1_pInt:product(grid)),[grid(1),grid(2),grid(3)]) - phaseField_real = 0.0_pReal - phaseField_real(1:grid(1),1:grid(2),1:grid(3)) = & - params%timeinc*params%phaseFieldData(i)%mobility* & - (phaseFieldRHS_lastInc(i,1:grid(1),1:grid(2),1:grid(3)) + & - phaseFieldRHS (i,1:grid(1),1:grid(2),1:grid(3)))/2.0_pReal - call utilities_scalarFFTforward() - call utilities_diffusion(params%phaseFieldData(i)%diffusion,params%timeinc/2.0_pReal) - call utilities_scalarFFTbackward() - f_scal(9+i,1:grid(1),1:grid(2),1:grid(3)) = & - x_scal(9+i,1:grid(1),1:grid(2),1:grid(3)) - & - f_scal(9+i,1:grid(1),1:grid(2),1:grid(3)) - & - phaseField_real(1:grid(1),1:grid(2),1:grid(3)) - err_phaseField(i) = maxval(abs(f_scal(9+i,1:grid(1),1:grid(2),1:grid(3)))) - phaseField_Avg(i) = sum(x_scal(9+i,1:grid(1),1:grid(2),1:grid(3)))*wgt - - case ('fracture') - - phaseField_real = 0.0_pReal - phaseField_real(1:grid(1),1:grid(2),1:grid(3)) = & - phaseField_lastInc(i,1:grid(1),1:grid(2),1:grid(3)) - call utilities_scalarFFTforward() - call utilities_diffusion(2.0_pReal*maxval(geomSize/real(grid,pReal))* & - params%phaseFieldData(i)%diffusion,params%timeinc) - call utilities_scalarFFTbackward() - f_scal(9+i,1:grid(1),1:grid(2),1:grid(3)) = & - phaseField_real(1:grid(1),1:grid(2),1:grid(3)) - - phaseFieldRHS(i,1:grid(1),1:grid(2),1:grid(3)) = & - - params%phaseFieldData(i)%mobility* & - sum(residual_F* & - (F-reshape(spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid(3)),[9,grid(1),grid(2),grid(3)])),dim=1) & - - params%phaseFieldData(i)%diffusion*(x_scal(9+i,1:grid(1),1:grid(2),1:grid(3)) - 1.0_pReal)/ & - 8.0_pReal/maxval(geomSize/real(grid,pReal)) - phaseField_real = 0.0_pReal - phaseField_real(1:grid(1),1:grid(2),1:grid(3)) = & - params%timeinc*params%phaseFieldData(i)%mobility* & - (phaseFieldRHS_lastInc(i,1:grid(1),1:grid(2),1:grid(3)) + & - phaseFieldRHS (i,1:grid(1),1:grid(2),1:grid(3)))/2.0_pReal - call utilities_scalarFFTforward() - call utilities_diffusion(2.0_pReal*maxval(geomSize/real(grid,pReal))* & - params%phaseFieldData(i)%diffusion,params%timeinc/2.0_pReal) - call utilities_scalarFFTbackward() - f_scal(9+i,1:grid(1),1:grid(2),1:grid(3)) = & - x_scal(9+i,1:grid(1),1:grid(2),1:grid(3)) - & - f_scal(9+i,1:grid(1),1:grid(2),1:grid(3)) - & - phaseField_real(1:grid(1),1:grid(2),1:grid(3)) - err_phaseField(i) = maxval(abs(f_scal(9+i,1:grid(1),1:grid(2),1:grid(3)))) - phaseField_Avg(i) = sum(x_scal(9+i,1:grid(1),1:grid(2),1:grid(3)))*wgt - - end select - enddo - !-------------------------------------------------------------------------------------------------- ! constructing residual - 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]) - + f_scal = reshape(field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3),shape(f_scal),order=[3,4,5,1,2]) + end subroutine BasicPETSc_formResidual @@ -668,17 +525,15 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du PetscErrorCode :: ierr real(pReal) :: & divTol, & - stressTol, & - phaseField_err = 0.0_pReal + stressTol divTol = max(maxval(abs(P_av))*err_div_tolRel,err_div_tolAbs) stressTol = max(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs) err_divPrev = err_div; err_div = err_divDummy - if (params%nActivePhaseFields .ne. 0_pInt) phaseField_err = maxval(err_phaseField/phaseField_Avg) converged: if ((totalIter >= itmin .and. & - all([ err_div/divTol, err_stress/stressTol] < 1.0_pReal) .and. & - phaseField_err < 1.0e-3_pReal) & + all([ err_div/divTol, & + err_stress/stressTol ] < 1.0_pReal)) & .or. terminallyIll) then reason = 1 elseif (totalIter >= itmax) then converged @@ -694,16 +549,11 @@ 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,')' - if (params%nActivePhaseFields .ne. 0_pInt) then - write(6,'(a,f10.2,a,es8.2,a,es9.2,a)') ' error phase field = ', & - phaseField_err/1.0e-3, ' (',phaseField_err, ' Pa, tol =',1.0e-3,')' - endif write(6,'(/,a)') ' ===========================================================================' flush(6) end subroutine BasicPETSc_converged - !-------------------------------------------------------------------------------------------------- !> @brief convergence check !-------------------------------------------------------------------------------------------------- diff --git a/code/DAMASK_spectral_utilities.f90 b/code/DAMASK_spectral_utilities.f90 index eab515c75..5e3e047d5 100644 --- a/code/DAMASK_spectral_utilities.f90 +++ b/code/DAMASK_spectral_utilities.f90 @@ -48,8 +48,6 @@ 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 @@ -70,8 +68,6 @@ 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 @@ -114,11 +110,8 @@ 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, & @@ -184,7 +177,6 @@ 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) @@ -247,11 +239,8 @@ 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) @@ -278,21 +267,6 @@ 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 @@ -308,7 +282,7 @@ subroutine utilities_init() endif if (debugFFTW) then - scalarField_realC = fftw_alloc_complex(int(product(grid),C_SIZE_T)) ! allllocate data for real representation (no in place transform) + scalarField_realC = fftw_alloc_complex(int(product(grid),C_SIZE_T)) ! allocate data for real representation (no in place transform) scalarField_fourierC = fftw_alloc_complex(int(product(grid),C_SIZE_T)) ! allocate data for fourier representation (no in place transform) call c_f_pointer(scalarField_realC, scalarField_real, grid) ! place a pointer for a real representation call c_f_pointer(scalarField_fourierC, scalarField_fourier, grid) ! place a pointer for a fourier representation @@ -521,59 +495,6 @@ 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 !-------------------------------------------------------------------------------------------------- @@ -607,36 +528,6 @@ 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 - -!-------------------------------------------------------------------------------------------------- -! 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 !--------------------------------------------------------------------------------------------------