diff --git a/code/DAMASK_spectral_solverAL.f90 b/code/DAMASK_spectral_solverAL.f90 index 739b4fd8a..43d9d51ae 100644 --- a/code/DAMASK_spectral_solverAL.f90 +++ b/code/DAMASK_spectral_solverAL.f90 @@ -1,20 +1,18 @@ !-------------------------------------------------------------------------------------------------- -! $Id: DAMASK_spectral_SolverAL.f90 1654 2012-08-03 09:25:48Z MPIE\m.diehl $ +! $Id: DAMASK_spectral_solverAL.f90 1654 2012-08-03 09:25:48Z MPIE\m.diehl $ !-------------------------------------------------------------------------------------------------- !> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @brief AL scheme solver !-------------------------------------------------------------------------------------------------- -module DAMASK_spectral_SolverAL +module DAMASK_spectral_solverAL use prec, only: & pInt, & pReal - use math, only: & math_I3 - - use DAMASK_spectral_Utilities, only: & + use DAMASK_spectral_utilities, only: & tSolutionState implicit none @@ -26,13 +24,14 @@ module DAMASK_spectral_SolverAL DAMASK_spectral_SolverAL_label = 'al' !-------------------------------------------------------------------------------------------------- -! derived types +! derived types ToDo: use here the type definition for a full loadcase including mask type tSolutionParams real(pReal), dimension(3,3) :: P_BC, rotation_BC real(pReal) :: timeinc end type tSolutionParams type(tSolutionParams), private :: params + real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal !-------------------------------------------------------------------------------------------------- ! PETSc data @@ -42,27 +41,34 @@ module DAMASK_spectral_SolverAL !-------------------------------------------------------------------------------------------------- ! common pointwise data - real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, F_lambda_lastInc, F_lambdaDot, Fdot + real(pReal), private, dimension(:,:,:,:,:), allocatable :: & + F_lastInc, & !< field of previous compatible deformation gradients + F_lambda_lastInc, & !< field of previous incompatible deformation gradient + Fdot, & !< field of assumed rate of compatible deformation gradient + F_lambdaDot !< field of assumed rate of incopatible deformation gradient real(pReal), private, dimension(:,:,:,:), allocatable :: coordinates - real(pReal), private, dimension(:,:,:), allocatable :: temperature + real(pReal), private :: temperature !< temperature, no spatial quantity at the moment !-------------------------------------------------------------------------------------------------- ! stress, stiffness and compliance average etc. real(pReal), private, dimension(3,3) :: & - F_aimDot, & - F_aim = math_I3, & - F_aim_lastInc = math_I3, & - P_av - character(len=1024), private :: incInfo + F_aimDot, & !< assumed rate of average deformation gradient + F_aim = math_I3, & !< current prescribed deformation gradient + F_aim_lastInc = math_I3, & !< previous average deformation gradient + P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress + character(len=1024), private :: incInfo !< time and increment information real(pReal), private, dimension(3,3,3,3) :: & - C = 0.0_pReal, C_lastInc = 0.0_pReal, & - S = 0.0_pReal, & - C_scale = 0.0_pReal, & + C = 0.0_pReal, & !< current average stiffness + C_lastInc = 0.0_pReal, & !< previous average stiffness + S = 0.0_pReal, & !< current compliance (filled up with zeros) + C_scale = 0.0_pReal, & S_scale = 0.0_pReal - real(pReal), private :: err_stress, err_f, err_p - logical, private :: ForwardData - real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal + real(pReal), private :: & + err_stress, & !< deviation from stress BC + err_f, & !< difference between compatible and incompatible deformation gradient + err_p !< difference of stress resulting from compatible and incompatible F + logical, private :: ForwardData integer(pInt) :: reportIter = 0_pInt contains @@ -70,215 +76,209 @@ contains !> @brief allocates all neccessary fields and fills them with data, potentially from restart info !-------------------------------------------------------------------------------------------------- subroutine AL_init() - 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_read_JobBinaryFile, & - IO_write_JobBinaryFile - - use FEsolving, only: & - restartInc + 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_read_JobBinaryFile, & + IO_write_JobBinaryFile + use FEsolving, only: & + restartInc + use DAMASK_interface, only: & + getSolverJobName + use DAMASK_spectral_Utilities, only: & + Utilities_init, & + Utilities_constitutiveResponse, & + Utilities_updateGamma, & + debugRestart + use numerics, only: & + petsc_options + use mesh, only: & + res, & + geomdim, & + mesh_NcpElems + use math, only: & + math_invSym3333 - use DAMASK_interface, only: & - getSolverJobName - - use DAMASK_spectral_Utilities, only: & - Utilities_init, & - Utilities_constitutiveResponse, & - Utilities_updateGamma, & - debugRestart - - use numerics, only: & - petsc_options - - use mesh, only: & - res, & - geomdim, & - mesh_NcpElems - - use math, only: & - math_invSym3333 - - implicit none + implicit none #include #include - integer(pInt) :: i,j,k - real(pReal), dimension(:,:,:,:,:), allocatable :: P - - PetscErrorCode :: ierr - PetscObject :: dummy - PetscScalar, pointer, dimension(:,:,:,:) :: xx_psc, F, F_lambda - - call Utilities_init() - - write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverAL init -+>>>' - write(6,'(a)') ' $Id: DAMASK_spectral_SolverAL.f90 1654 2012-08-03 09:25:48Z MPIE\m.diehl $' -#include "compilation_info.f90" - write(6,'(a)') '' - - allocate (F_lastInc (3,3, res(1), res(2),res(3)), source = 0.0_pReal) - allocate (Fdot (3,3, res(1), res(2),res(3)), source = 0.0_pReal) - ! allocate (Fdot,source = F_lastInc) somethin like that should be possible - allocate (F_lambda_lastInc(3,3, res(1), res(2),res(3)), source = 0.0_pReal) - allocate (F_lambdaDot(3,3, res(1), res(2),res(3)), source = 0.0_pReal) - allocate (P (3,3, res(1), res(2),res(3)), source = 0.0_pReal) - allocate (coordinates( res(1), res(2),res(3),3), source = 0.0_pReal) - allocate (temperature( res(1), res(2),res(3)), source = 0.0_pReal) - - !-------------------------------------------------------------------------------------------------- - ! PETSc Init - call SNESCreate(PETSC_COMM_WORLD,snes,ierr) - call DMDACreate3d(PETSC_COMM_WORLD, & - DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, & - DMDA_STENCIL_BOX,res(1),res(2),res(3),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, & - 18,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr) - call DMCreateGlobalVector(da,solution_vec,ierr) - call DMDASetLocalFunction(da,AL_formResidual,ierr) - call SNESSetDM(snes,da,ierr) - call SNESSetConvergenceTest(snes,AL_converged,dummy,PETSC_NULL_FUNCTION,ierr) - call SNESSetFromOptions(snes,ierr) - - !-------------------------------------------------------------------------------------------------- - ! init fields - call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr) - F => xx_psc(0:8,:,:,:) - F_lambda => xx_psc(9:17,:,:,:) - if (restartInc == 1_pInt) then ! no deformation (no restart) - F_lastInc = spread(spread(spread(math_I3,3,res(1)),4,res(2)),5,res(3)) ! initialize to identity - F_lambda_lastInc = F_lastInc - F = reshape(F_lastInc,[9,res(1),res(2),res(3)]) - F_lambda = F - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - coordinates(i,j,k,1:3) = geomdim/real(res,pReal)*real([i,j,k],pReal) & - - geomdim/real(2_pInt*res,pReal) - enddo; enddo; enddo - elseif (restartInc > 1_pInt) then ! using old values from file - if (debugRestart) write(6,'(a,i6,a)') 'Reading values of increment ',& - restartInc - 1_pInt,' from file' - call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',& - trim(getSolverJobName()),size(F_lastInc)) - read (777,rec=1) F - close (777) - call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc',& - trim(getSolverJobName()),size(F_lastInc)) - read (777,rec=1) F_lastInc - close (777) - call IO_read_jobBinaryFile(777,'convergedSpectralDefgradLambda',& - trim(getSolverJobName()),size(F_lambda_lastInc)) - read (777,rec=1) F_lambda - close (777) - call IO_read_jobBinaryFile(777,'convergedSpectralDefgradLambda_lastInc',& - trim(getSolverJobName()),size(F_lambda_lastInc)) - read (777,rec=1) F_lambda_lastInc - close (777) - call IO_read_jobBinaryFile(777,'F_aim',trim(getSolverJobName()),size(F_aim)) - read (777,rec=1) F_aim - close (777) - call IO_read_jobBinaryFile(777,'F_aim_lastInc',trim(getSolverJobName()),size(F_aim_lastInc)) - read (777,rec=1) F_aim_lastInc - close (777) - - coordinates = 0.0 ! change it later!!! - endif - !print*, shape(F) - call Utilities_constitutiveResponse(coordinates,F,F,temperature,0.0_pReal,P,C,P_av,.false.,math_I3) - call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr) - - !-------------------------------------------------------------------------------------------------- - ! reference stiffness - if (restartInc == 1_pInt) then - call IO_write_jobBinaryFile(777,'C_ref',size(C)) - write (777,rec=1) C - close(777) - elseif (restartInc > 1_pInt) then - call IO_read_jobBinaryFile(777,'C_ref',trim(getSolverJobName()),size(C)) - read (777,rec=1) C - close (777) - endif - - call Utilities_updateGamma(C,.True.) - C_scale = C - S_scale = math_invSym3333(C) + integer(pInt) :: i,j,k + real(pReal), dimension(3,3, res(1), res(2),res(3)) :: P - end subroutine AL_init - + PetscErrorCode :: ierr + PetscObject :: dummy + PetscScalar, pointer, dimension(:,:,:,:) :: xx_psc, F, F_lambda + + call Utilities_init() + write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverAL init -+>>>' + write(6,'(a)') ' $Id: DAMASK_spectral_SolverAL.f90 1654 2012-08-03 09:25:48Z MPIE\m.diehl $' +#include "compilation_info.f90" + write(6,'(a)') '' + + allocate (F_lastInc (3,3, res(1), res(2),res(3)), source = 0.0_pReal) + allocate (Fdot (3,3, res(1), res(2),res(3)), source = 0.0_pReal) + ! allocate (Fdot,source = F_lastInc) somethin like that should be possible + allocate (F_lambda_lastInc(3,3, res(1), res(2),res(3)), source = 0.0_pReal) + allocate (F_lambdaDot(3,3, res(1), res(2),res(3)), source = 0.0_pReal) + allocate (coordinates( res(1), res(2),res(3),3), source = 0.0_pReal) + +!-------------------------------------------------------------------------------------------------- +! PETSc Init + call SNESCreate(PETSC_COMM_WORLD,snes,ierr) + CHKERRQ(ierr) + call DMDACreate3d(PETSC_COMM_WORLD, & + DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, & + DMDA_STENCIL_BOX,res(1),res(2),res(3),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, & + 18,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr) + CHKERRQ(ierr) + call DMCreateGlobalVector(da,solution_vec,ierr) + CHKERRQ(ierr) + call DMDASetLocalFunction(da,AL_formResidual,ierr) + CHKERRQ(ierr) + call SNESSetDM(snes,da,ierr) + CHKERRQ(ierr) + call SNESSetConvergenceTest(snes,AL_converged,dummy,PETSC_NULL_FUNCTION,ierr) + CHKERRQ(ierr) + call SNESSetFromOptions(snes,ierr) + CHKERRQ(ierr) + +!-------------------------------------------------------------------------------------------------- +! init fields + call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr) ! places pointer xx_psc on PETSc data + F => xx_psc(0:8,:,:,:) + F_lambda => xx_psc(9:17,:,:,:) + if (restartInc == 1_pInt) then ! no deformation (no restart) + F_lastInc = spread(spread(spread(math_I3,3,res(1)),4,res(2)),5,res(3)) ! initialize to identity + F_lambda_lastInc = F_lastInc + F = reshape(F_lastInc,[9,res(1),res(2),res(3)]) + F_lambda = F + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + coordinates(i,j,k,1:3) = geomdim/real(res,pReal)*real([i,j,k],pReal) & + - geomdim/real(2_pInt*res,pReal) + enddo; enddo; enddo + elseif (restartInc > 1_pInt) then ! using old values from file + if (debugRestart) write(6,'(a,i6,a)') 'Reading values of increment ',& + restartInc - 1_pInt,' from file' + call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',& + trim(getSolverJobName()),size(F_lastInc)) + read (777,rec=1) F + close (777) + call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc',& + trim(getSolverJobName()),size(F_lastInc)) + read (777,rec=1) F_lastInc + close (777) + call IO_read_jobBinaryFile(777,'convergedSpectralDefgradLambda',& + trim(getSolverJobName()),size(F_lambda_lastInc)) + read (777,rec=1) F_lambda + close (777) + call IO_read_jobBinaryFile(777,'convergedSpectralDefgradLambda_lastInc',& + trim(getSolverJobName()),size(F_lambda_lastInc)) + read (777,rec=1) F_lambda_lastInc + close (777) + call IO_read_jobBinaryFile(777,'F_aim',trim(getSolverJobName()),size(F_aim)) + read (777,rec=1) F_aim + close (777) + call IO_read_jobBinaryFile(777,'F_aim_lastInc',trim(getSolverJobName()),size(F_aim_lastInc)) + read (777,rec=1) F_aim_lastInc + close (777) + coordinates = 0.0 ! change it later!!! + endif + call Utilities_constitutiveResponse(coordinates,F,F,temperature,0.0_pReal,P,C,P_av,.false.,math_I3) + call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr) + CHKERRQ(ierr) + +!-------------------------------------------------------------------------------------------------- +! reference stiffness + if (restartInc == 1_pInt) then + call IO_write_jobBinaryFile(777,'C_ref',size(C)) + write (777,rec=1) C + close(777) + elseif (restartInc > 1_pInt) then + call IO_read_jobBinaryFile(777,'C_ref',trim(getSolverJobName()),size(C)) + read (777,rec=1) C + close (777) + endif + + call Utilities_updateGamma(C,.True.) + C_scale = C + S_scale = math_invSym3333(C) + +end subroutine AL_init + + !-------------------------------------------------------------------------------------------------- !> @brief solution for the AL scheme with internal iterations !-------------------------------------------------------------------------------------------------- - type(tSolutionState) function & +type(tSolutionState) function & AL_solution(incInfoIn,guess,timeinc,timeinc_old,P_BC,F_BC,temperature_bc,rotation_BC) - - use numerics, only: & - update_gamma - use math, only: & - math_mul33x33 ,& - math_rotate_backward33 - use mesh, only: & - res,& - geomdim,& - deformed_fft - use IO, only: & - IO_write_JobBinaryFile - use DAMASK_spectral_Utilities, only: & - tBoundaryCondition, & - Utilities_forwardField, & - Utilities_calculateRate, & - Utilities_maskedCompliance, & - Utilities_updateGamma, & - cutBack - - use FEsolving, only: & - restartWrite, & - terminallyIll - - implicit none + + use numerics, only: & + update_gamma + use math, only: & + math_mul33x33 ,& + math_rotate_backward33 + use mesh, only: & + res,& + geomdim,& + deformed_fft + use IO, only: & + IO_write_JobBinaryFile + use DAMASK_spectral_Utilities, only: & + tBoundaryCondition, & + Utilities_forwardField, & + Utilities_calculateRate, & + Utilities_maskedCompliance, & + Utilities_updateGamma, & + cutBack + use FEsolving, only: & + restartWrite, & + terminallyIll + + implicit none #include #include !-------------------------------------------------------------------------------------------------- ! input data for solution - real(pReal), intent(in) :: timeinc, timeinc_old, temperature_bc - logical, intent(in) :: guess - type(tBoundaryCondition), intent(in) :: P_BC,F_BC - character(len=*), intent(in) :: incInfoIn - real(pReal), dimension(3,3), intent(in) :: rotation_BC - real(pReal), dimension(3,3) ,save :: F_aimDot - real(pReal), dimension(3,3) :: F_aim_lab + real(pReal), intent(in) :: timeinc, timeinc_old, temperature_bc + logical, intent(in) :: guess + type(tBoundaryCondition), intent(in) :: P_BC,F_BC + character(len=*), intent(in) :: incInfoIn + real(pReal), dimension(3,3), intent(in) :: rotation_BC + real(pReal), dimension(3,3) ,save :: F_aimDot + real(pReal), dimension(3,3) :: F_aim_lab + !-------------------------------------------------------------------------------------------------- ! loop variables, convergence etc. - real(pReal), dimension(3,3) :: temp33_Real + real(pReal), dimension(3,3) :: temp33_Real !-------------------------------------------------------------------------------------------------- -! - PetscScalar, dimension(:,:,:,:), pointer :: xx_psc, F, F_lambda - PetscErrorCode :: ierr - SNESConvergedReason ::reason +! PETSc Data + PetscScalar, dimension(:,:,:,:), pointer :: xx_psc, F, F_lambda + PetscErrorCode :: ierr + SNESConvergedReason ::reason !-------------------------------------------------------------------------------------------------- ! restart information for spectral solver - incInfo = incInfoIn - if (restartWrite) then - write(6,'(a)') 'writing converged results for restart' - call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(F_lastInc)) - write (777,rec=1) F_LastInc - close (777) - call IO_write_jobBinaryFile(777,'C',size(C)) - write (777,rec=1) C - close(777) - endif - AL_solution%converged =.false. - call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr) - F => xx_psc(0:8,:,:,:) - F_lambda => xx_psc(9:17,:,:,:) + incInfo = incInfoIn + if (restartWrite) then + write(6,'(a)') 'writing converged results for restart' + call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(F_lastInc)) + write (777,rec=1) F_LastInc + close (777) + call IO_write_jobBinaryFile(777,'C',size(C)) + write (777,rec=1) C + close(777) + endif + AL_solution%converged =.false. + call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr) + F => xx_psc(0:8,:,:,:) + F_lambda => xx_psc(9:17,:,:,:) - -if ( cutBack) then - F_aim = F_aim_lastInc - F_lambda= reshape(F_lambda_lastInc,[9,res(1),res(2),res(3)]) - F = reshape(F_lastInc,[9,res(1),res(2),res(3)]) - C = C_lastInc -else - -!-------------------------------------------------------------------------------------------------- + if ( cutBack) then + F_aim = F_aim_lastInc + F_lambda= reshape(F_lambda_lastInc,[9,res(1),res(2),res(3)]) + F = reshape(F_lastInc,[9,res(1),res(2),res(3)]) + C = C_lastInc + else C_lastInc = C !-------------------------------------------------------------------------------------------------- ! calculate rate for aim @@ -287,7 +287,7 @@ else elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed f_aimDot = F_BC%maskFloat * F_BC%values endif - if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old + if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old F_aim_lastInc = F_aim !-------------------------------------------------------------------------------------------------- @@ -300,94 +300,97 @@ else F_lastInc = reshape(F,[3,3,res(1),res(2),res(3)]) F_lambda_lastInc = reshape(F_lambda,[3,3,res(1),res(2),res(3)]) - call deformed_fft(res,geomdim,math_rotate_backward33(F_aim_lastInc,rotation_BC), & + call deformed_fft(res,geomdim,math_rotate_backward33(F_aim_lastInc,rotation_BC), & 1.0_pReal,F_lastInc,coordinates) endif F_aim = F_aim + f_aimDot * timeinc !-------------------------------------------------------------------------------------------------- ! update local deformation gradient and coordinates - ! deltaF_aim = math_rotate_backward33(deltaF_aim,rotation_BC) +! deltaF_aim = math_rotate_backward33(deltaF_aim,rotation_BC) + F = reshape(Utilities_forwardField(timeinc,F_aim,F_lastInc,Fdot),[9,res(1),res(2),res(3)]) + F_lambda = reshape(Utilities_forwardField(timeinc,F_aim,F_lambda_lastInc,F_lambdadot),[9,res(1),res(2),res(3)]) - F = reshape(Utilities_forwardField(timeinc,F_aim,F_lastInc,Fdot),[9,res(1),res(2),res(3)]) - F_lambda = reshape(Utilities_forwardField(timeinc,F_aim,F_lambda_lastInc,F_lambdadot),[9,res(1),res(2),res(3)]) + call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr) + CHKERRQ(ierr) - call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr) - !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C) - if (update_gamma) call Utilities_updateGamma(C,restartWrite) - - ForwardData = .True. - mask_stress = P_BC%maskFloat - params%P_BC = P_BC%values - params%rotation_BC = rotation_BC - params%timeinc = timeinc + S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C) + if (update_gamma) call Utilities_updateGamma(C,restartWrite) + + ForwardData = .True. + mask_stress = P_BC%maskFloat + params%P_BC = P_BC%values + params%rotation_BC = rotation_BC + params%timeinc = timeinc - call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr) - call SNESGetConvergedReason(snes,reason,ierr) - - AL_solution%termIll = terminallyIll - terminallyIll = .false. - + call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr) + CHKERRQ(ierr) + call SNESGetConvergedReason(snes,reason,ierr) + CHKERRQ(ierr) + + AL_solution%termIll = terminallyIll + terminallyIll = .false. + if (reason > 0 ) then + AL_solution%converged = .true. + AL_solution%iterationsNeeded = reportIter + endif + +end function AL_solution - if (reason > 0 ) then - AL_solution%converged = .true. - AL_solution%iterationsNeeded = reportIter - endif - end function AL_solution !-------------------------------------------------------------------------------------------------- !> @brief forms the AL residual vector !-------------------------------------------------------------------------------------------------- - subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) - - use numerics, only: & - itmax, & - itmin - use math, only: & - math_rotate_backward33, & - math_transpose33, & - math_mul3333xx33 - use mesh, only: & - res, & - wgt - use DAMASK_spectral_Utilities, only: & - field_real, & - Utilities_FFTforward, & - Utilities_fourierConvolution, & - Utilities_FFTbackward, & - Utilities_constitutiveResponse - use IO, only: IO_intOut - - implicit none +subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) - integer(pInt) :: i,j,k - integer(pInt), save :: callNo = 3_pInt - real(pReal), dimension(3,3) :: temp33_Real - logical :: report - - DMDALocalInfo :: in(DMDA_LOCAL_INFO_SIZE) - PetscScalar, target :: x_scal(3,3,2,XG_RANGE,YG_RANGE,ZG_RANGE) - PetscScalar, target :: f_scal(3,3,2,X_RANGE,Y_RANGE,Z_RANGE) - PetscScalar, pointer :: F(:,:,:,:,:), F_lambda(:,:,:,:,:) - PetscScalar, pointer :: residual_F(:,:,:,:,:), residual_F_lambda(:,:,:,:,:) - PetscInt :: iter, nfuncs - PetscObject :: dummy - PetscErrorCode :: ierr + use numerics, only: & + itmax, & + itmin + use math, only: & + math_rotate_backward33, & + math_transpose33, & + math_mul3333xx33 + use mesh, only: & + res, & + wgt + use DAMASK_spectral_Utilities, only: & + field_real, & + Utilities_FFTforward, & + Utilities_fourierConvolution, & + Utilities_FFTbackward, & + Utilities_constitutiveResponse + use IO, only: IO_intOut - F => x_scal(:,:,1,:,:,:) - F_lambda => x_scal(:,:,2,:,:,:) - residual_F => f_scal(:,:,1,:,:,:) - residual_F_lambda => f_scal(:,:,2,:,:,:) - - call SNESGetNumberFunctionEvals(snes,nfuncs,ierr) - call SNESGetIterationNumber(snes,iter,ierr) - - !-------------------------------------------------------------------------------------------------- - ! report begin of new iteration + implicit none + integer(pInt) :: i,j,k + integer(pInt), save :: callNo = 3_pInt + real(pReal), dimension(3,3) :: temp33_Real + logical :: report + + DMDALocalInfo :: in(DMDA_LOCAL_INFO_SIZE) + PetscScalar, target :: x_scal(3,3,2,XG_RANGE,YG_RANGE,ZG_RANGE) + PetscScalar, target :: f_scal(3,3,2,X_RANGE,Y_RANGE,Z_RANGE) + PetscScalar, pointer :: F(:,:,:,:,:), F_lambda(:,:,:,:,:) + PetscScalar, pointer :: residual_F(:,:,:,:,:), residual_F_lambda(:,:,:,:,:) + PetscInt :: iter, nfuncs + PetscObject :: dummy + PetscErrorCode :: ierr + + F => x_scal(:,:,1,:,:,:) + F_lambda => x_scal(:,:,2,:,:,:) + residual_F => f_scal(:,:,1,:,:,:) + residual_F_lambda => f_scal(:,:,2,:,:,:) + + call SNESGetNumberFunctionEvals(snes,nfuncs,ierr) + CHKERRQ(ierr) + call SNESGetIterationNumber(snes,iter,ierr) + CHKERRQ(ierr) + +!-------------------------------------------------------------------------------------------------- +! report begin of new iteration if (iter == 0 .and. callNo>2) then callNo = 0_pInt reportIter = 0_pInt @@ -400,105 +403,112 @@ else reportIter = reportIter + 1_pInt endif callNo = callNo +1_pInt - !-------------------------------------------------------------------------------------------------- - ! evaluate constitutive response - call Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,params%timeinc, & - residual_F,C,P_av,ForwardData,params%rotation_BC) - ForwardData = .False. - +!-------------------------------------------------------------------------------------------------- +! evaluate constitutive response + call Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,params%timeinc, & + residual_F,C,P_av,ForwardData,params%rotation_BC) + ForwardData = .False. + !-------------------------------------------------------------------------------------------------- ! 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 + +!-------------------------------------------------------------------------------------------------- +! + field_real = 0.0_pReal + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + temp33_Real = math_mul3333xx33(S_scale,residual_F(1:3,1:3,i,j,k)) + math_I3 + residual_F(1:3,1:3,i,j,k) = temp33_Real + field_real(i,j,k,1:3,1:3) = -math_mul3333xx33(C_scale,F_lambda(1:3,1:3,i,j,k)-F(1:3,1:3,i,j,k)) + enddo; enddo; enddo + +!-------------------------------------------------------------------------------------------------- +! doing convolution in Fourier space + call Utilities_FFTforward() + call Utilities_fourierConvolution(math_rotate_backward33(F_aim,params%rotation_BC)) + call Utilities_FFTbackward() + +!-------------------------------------------------------------------------------------------------- +! constructing residual + residual_F_lambda = F - reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),& + [3,3,res(1),res(2),res(3)],order=[3,4,5,1,2]) + residual_F = residual_F - F_lambda + residual_F_lambda + +!-------------------------------------------------------------------------------------------------- +! calculating errors + err_f = wgt*sqrt(sum(residual_F_lambda**2.0_pReal)) + err_p = wgt*sqrt(sum((residual_F - residual_F_lambda)**2.0_pReal)) + +end subroutine AL_formResidual - - !-------------------------------------------------------------------------------------------------- - ! doing Fourier transform - field_real = 0.0_pReal - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - temp33_Real = math_mul3333xx33(S_scale,residual_F(1:3,1:3,i,j,k)) + math_I3 - residual_F(1:3,1:3,i,j,k) = temp33_Real - field_real(i,j,k,1:3,1:3) = -math_mul3333xx33(C_scale,F_lambda(1:3,1:3,i,j,k)-F(1:3,1:3,i,j,k)) - enddo; enddo; enddo - - !-------------------------------------------------------------------------------------------------- - ! doing Fourier transform - call Utilities_FFTforward() - call Utilities_fourierConvolution(math_rotate_backward33(F_aim,params%rotation_BC)) - call Utilities_FFTbackward() - - !-------------------------------------------------------------------------------------------------- - ! constructing residual - residual_F_lambda = F - reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),& - [3,3,res(1),res(2),res(3)],order=[3,4,5,1,2]) - residual_F = residual_F - F_lambda + residual_F_lambda - - !-------------------------------------------------------------------------------------------------- - ! calculating errors - err_f = wgt*sqrt(sum(residual_F_lambda**2.0_pReal)) - err_p = wgt*sqrt(sum((residual_F - residual_F_lambda)**2.0_pReal)) - - end subroutine AL_formResidual !-------------------------------------------------------------------------------------------------- !> @brief convergence check !-------------------------------------------------------------------------------------------------- - subroutine AL_converged(snes_local,it,xnorm,snorm,fnorm,reason,dummy,ierr) - - use numerics, only: & - itmax, & - itmin, & - err_f_tol, & - err_p_tol, & - err_stress_tolrel, & - err_stress_tolabs - - implicit none +subroutine AL_converged(snes_local,it,xnorm,snorm,fnorm,reason,dummy,ierr) - SNES :: snes_local - PetscInt :: it - PetscReal :: xnorm, snorm, fnorm - SNESConvergedReason :: reason - PetscObject :: dummy - PetscErrorCode ::ierr - logical :: Converged - - Converged = (it > itmin .and. & - all([ err_f/sqrt(sum((F_aim-math_I3)*(F_aim-math_I3)))/err_f_tol, & - err_p/sqrt(sum((F_aim-math_I3)*(F_aim-math_I3)))/err_p_tol, & - err_stress/min(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs)] < 1.0_pReal)) - - if (Converged) then - reason = 1 - elseif (it > itmax) then - reason = -1 - else - reason = 0 - endif - - write(6,'(a,es14.7)') 'error stress BC = ', err_stress/min(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs) - write(6,'(a,es14.7)') 'error F = ', err_f/sqrt(sum((F_aim-math_I3)*(F_aim-math_I3)))/err_f_tol - write(6,'(a,es14.7)') 'error P = ', err_p/sqrt(sum((F_aim-math_I3)*(F_aim-math_I3)))/err_p_tol - write(6,'(/,a)') '==========================================================================' - end subroutine AL_converged + use numerics, only: & + itmax, & + itmin, & + err_f_tol, & + err_p_tol, & + err_stress_tolrel, & + err_stress_tolabs + + implicit none + + SNES :: snes_local + PetscInt :: it + PetscReal :: xnorm, snorm, fnorm + SNESConvergedReason :: reason + PetscObject :: dummy + PetscErrorCode ::ierr + logical :: Converged + + Converged = (it > itmin .and. & + all([ err_f/sqrt(sum((F_aim-math_I3)*(F_aim-math_I3)))/err_f_tol, & + err_p/sqrt(sum((F_aim-math_I3)*(F_aim-math_I3)))/err_p_tol, & + err_stress/min(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs)] < 1.0_pReal)) + + if (Converged) then + reason = 1 + elseif (it > itmax) then + reason = -1 + else + reason = 0 + endif + + write(6,'(a,es14.7)') 'error stress BC = ', & + err_stress/min(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs) + write(6,'(a,es14.7)') 'error F = ',& + err_f/sqrt(sum((F_aim-math_I3)*(F_aim-math_I3)))/err_f_tol + write(6,'(a,es14.7)') 'error P = ', & + err_p/sqrt(sum((F_aim-math_I3)*(F_aim-math_I3)))/err_p_tol + write(6,'(/,a)') '==========================================================================' +end subroutine AL_converged !-------------------------------------------------------------------------------------------------- !> @brief destroy routine !-------------------------------------------------------------------------------------------------- - subroutine AL_destroy() - - use DAMASK_spectral_Utilities, only: & - Utilities_destroy - implicit none - PetscErrorCode :: ierr - - call VecDestroy(solution_vec,ierr) - call SNESDestroy(snes,ierr) - call DMDestroy(da,ierr) - call PetscFinalize(ierr) - call Utilities_destroy() +subroutine AL_destroy() - end subroutine AL_destroy + use DAMASK_spectral_Utilities, only: & + Utilities_destroy + implicit none + PetscErrorCode :: ierr - end module DAMASK_spectral_SolverAL + call VecDestroy(solution_vec,ierr) + CHKERRQ(ierr) + call SNESDestroy(snes,ierr) + CHKERRQ(ierr) + call DMDestroy(da,ierr) + CHKERRQ(ierr) + call PetscFinalize(ierr) + CHKERRQ(ierr) + call Utilities_destroy() + +end subroutine AL_destroy + +end module DAMASK_spectral_SolverAL diff --git a/code/DAMASK_spectral_solverBasic.f90 b/code/DAMASK_spectral_solverBasic.f90 index 67300dec6..b5d9b3409 100644 --- a/code/DAMASK_spectral_solverBasic.f90 +++ b/code/DAMASK_spectral_solverBasic.f90 @@ -12,10 +12,8 @@ module DAMASK_spectral_SolverBasic use prec, only: & pInt, & pReal - use math, only: & math_I3 - use DAMASK_spectral_Utilities, only: & tSolutionState @@ -24,19 +22,23 @@ module DAMASK_spectral_SolverBasic DAMASK_spectral_SolverBasic_label = 'basic' !-------------------------------------------------------------------------------------------------- -! pointwise data - real(pReal), private, dimension(:,:,:,:,:), allocatable :: F, F_lastInc, Fdot, P +! pointwise global data + real(pReal), private, dimension(:,:,:,:,:), allocatable :: & + F, & !< deformation gradient field + F_lastInc, & !< deformation gradient field last increment + Fdot !< assumed rate for F n to F n+1 real(pReal), private, dimension(:,:,:,:), allocatable :: coordinates - real(pReal), private, dimension(:,:,:), allocatable :: temperature + real(pReal), private :: temperature !< not pointwise !-------------------------------------------------------------------------------------------------- ! stress, stiffness and compliance average etc. real(pReal), private, dimension(3,3) :: & - F_aim = math_I3, & - F_aim_lastInc = math_I3, & - F_aimDot = 0.0_pReal + F_aim = math_I3, & !< deformation gradient aim + F_aim_lastInc = math_I3, & !< deformation gradient aim last increment + F_aimDot = 0.0_pReal !< assumed rate real(pReal), private,dimension(3,3,3,3) :: & - C = 0.0_pReal, C_lastInc = 0.0_pReal + C = 0.0_pReal, & !< average stiffness + C_lastInc = 0.0_pReal !< average stiffness last increment contains @@ -50,27 +52,27 @@ subroutine basic_init() IO_read_JobBinaryFile, & IO_write_JobBinaryFile, & IO_intOut - use FEsolving, only: & restartInc - use DAMASK_interface, only: & getSolverJobName - use DAMASK_spectral_Utilities, only: & Utilities_init, & Utilities_constitutiveResponse, & Utilities_updateGamma, & debugRestart - use mesh, only: & res, & geomdim implicit none - integer(pInt) :: i,j,k - real(pReal), dimension(3,3) :: temp33_Real - real(pReal), dimension(3,3,3,3) :: temp3333_Real + real(pReal), dimension(3,3,res(1),res(2),res(3)) :: P + integer(pInt) :: & + i, j, k + real(pReal), dimension(3,3) :: & + temp33_Real + real(pReal), dimension(3,3,3,3) :: & + temp3333_Real call Utilities_Init() @@ -79,15 +81,15 @@ subroutine basic_init() #include "compilation_info.f90" write(6,'(a)') '' +!-------------------------------------------------------------------------------------------------- +! allocate global fields allocate (F ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal) allocate (F_lastInc ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal) allocate (Fdot ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal) - allocate (P ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal) allocate (coordinates( res(1), res(2),res(3),3),source = 0.0_pReal) - allocate (temperature( res(1), res(2),res(3)), source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- -! init fields +! init fields and average quantities if (restartInc == 1_pInt) then ! no deformation (no restart) F = spread(spread(spread(math_I3,3,res(1)),4,res(2)),5,res(3)) ! initialize to identity F_lastInc = F @@ -127,13 +129,9 @@ subroutine basic_init() close (777) coordinates = 0.0 ! change it later!!! endif - call Utilities_constitutiveResponse(coordinates,F,F,temperature,0.0_pReal,& - P,C,temp33_Real,.false.,math_I3) + call Utilities_constitutiveResponse(coordinates,F,F,temperature,0.0_pReal,P,C,temp33_Real,.false.,math_I3) ! constitutive response with no deformation in no time to get reference stiffness !no rotation bc call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,rotation_BC),1.0_pReal,F_lastInc,coordinates) - -!-------------------------------------------------------------------------------------------------- -! reference stiffness if (restartInc == 1_pInt) then ! use initial stiffness as reference stiffness temp3333_Real = C endif @@ -148,7 +146,6 @@ end subroutine basic_init !-------------------------------------------------------------------------------------------------- type(tSolutionState) function & basic_solution(incInfo,guess,timeinc,timeinc_old,P_BC,F_BC,temperature_bc,rotation_BC) - use numerics, only: & itmax, & itmin, & @@ -162,11 +159,10 @@ type(tSolutionState) function & res,& geomdim, & deformed_fft, & - wgt + wgt use IO, only: & IO_write_JobBinaryFile, & IO_intOut - use DAMASK_spectral_Utilities, only: & tBoundaryCondition, & field_real, & @@ -179,7 +175,6 @@ type(tSolutionState) function & Utilities_updateGamma, & Utilities_constitutiveResponse, & Utilities_calculateRate - use FEsolving, only: & restartWrite, & restartRead, & @@ -190,26 +185,40 @@ type(tSolutionState) function & implicit none !-------------------------------------------------------------------------------------------------- ! input data for solution - real(pReal), intent(in) :: timeinc, timeinc_old, temperature_bc - logical, intent(in) :: guess - type(tBoundaryCondition), intent(in) :: P_BC,F_BC - character(len=*), intent(in) :: incInfo - real(pReal), dimension(3,3), intent(in) :: rotation_BC + real(pReal), intent(in) :: & + timeinc, & !< increment in time for current solution + timeinc_old, & !< increment in time of last increment + temperature_bc !< temperature (I wonder, why it is intent(in) + logical, intent(in) :: & + guess !< if .false., assume homogeneous addon + type(tBoundaryCondition), intent(in) :: & + P_BC,& !< stress boundary conditions + F_BC !< deformation boundary conditions + character(len=*), intent(in) :: & + incInfo !< string with information of current increment for output to screen + real(pReal), dimension(3,3), intent(in) :: & + rotation_BC !< rotation load to lab - real(pReal), dimension(3,3,3,3) :: S - real(pReal), dimension(3,3) :: F_aim_lab, & - F_aim_lab_lastIter, & - P_av + real(pReal), dimension(3,3,3,3) :: & + S !< current average compliance + real(pReal), dimension(3,3) :: & + F_aim_lab, & + F_aim_lab_lastIter, & !< aim of last iteration + P_av + real(pReal), dimension(3,3,res(1),res(2),res(3)) :: P !-------------------------------------------------------------------------------------------------- ! loop variables, convergence etc. real(pReal) :: err_div, err_stress integer(pInt) :: iter, row, column, i, j, k logical :: ForwardData - real(pReal) :: defgradDet, defgradDetMax, defgradDetMin + real(pReal) :: & + defgradDet, & + defgradDetMax, & + defgradDetMin real(pReal), dimension(3,3) :: temp33_Real !-------------------------------------------------------------------------------------------------- -! restart information for spectral solver +! write restart information for spectral solver if (restartWrite) then write(6,'(a)') 'writing converged results for restart' call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(F)) ! writing deformation gradient field to file @@ -235,14 +244,14 @@ type(tSolutionState) function & close(777) endif - if ( cutBack) then +!-------------------------------------------------------------------------------------------------- +! forward data + if (cutBack) then F_aim = F_aim_lastInc F = F_lastInc C = C_lastInc else - C_lastInc = C -!-------------------------------------------------------------------------------------------------- -! calculate rate for aim + C_lastInc = C if (F_BC%myType=='l') then ! calculate f_aimDot from given L and current F f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim) elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed @@ -251,9 +260,7 @@ type(tSolutionState) function & if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old F_aim_lastInc = F_aim -!-------------------------------------------------------------------------------------------------- -! update coordinates and rate and forward last inc - call deformed_fft(res,geomdim,math_rotate_backward33(F_aim_lastInc,rotation_BC), & + call deformed_fft(res,geomdim,math_rotate_backward33(F_aim_lastInc,rotation_BC), & ! update coordinates and rate and forward last inc 1.0_pReal,F_lastInc,coordinates) Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & timeinc,timeinc_old,guess,F_lastInc,F) @@ -265,13 +272,13 @@ type(tSolutionState) function & !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C) - if (update_gamma) call Utilities_updateGamma(C,restartWrite) +!-------------------------------------------------------------------------------------------------- +! iteration till converged if (.not. restartRead) ForwardData = .True. iter = 0_pInt convergenceLoop: do while(iter < itmax) - iter = iter + 1_pInt !-------------------------------------------------------------------------------------------------- ! report begin of new iteration @@ -279,16 +286,16 @@ type(tSolutionState) function & ' @ Iter. ', itmin, '≤',iter, '≤', itmax write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim =', & math_transpose33(F_aim) - F_aim_lab_lastIter = math_rotate_backward33(F_aim,rotation_BC) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response + F_aim_lab_lastIter = math_rotate_backward33(F_aim,rotation_BC) basic_solution%termIll = .false. call Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,timeinc,& P,C,P_av,ForwardData,rotation_BC) basic_solution%termIll = terminallyIll - terminallyIll = .False. - ForwardData = .False. + terminallyIll = .false. + ForwardData = .false. !-------------------------------------------------------------------------------------------------- ! stress BC handling @@ -326,7 +333,6 @@ logical function basic_Converged(err_div,pAvgDiv,err_stress,pAvgStress) err_div_tol, & err_stress_tolrel, & err_stress_tolabs - use math, only: & math_mul33x33, & math_eigenvalues33, & @@ -344,8 +350,6 @@ logical function basic_Converged(err_div,pAvgDiv,err_stress,pAvgStress) real(pReal) :: & err_stress_tol, & pAvgDivL2 - - pAvgDivL2 = sqrt(maxval(math_eigenvalues33(math_mul33x33(pAvgDiv,math_transpose33(pAvgDiv))))) ! L_2 norm of average stress (http://mathworld.wolfram.com/SpectralNorm.html) err_stress_tol = min(maxval(abs(pAvgStress))*err_stress_tolrel,err_stress_tolabs) @@ -363,8 +367,8 @@ end function basic_Converged subroutine basic_destroy() - use DAMASK_spectral_Utilities, only: & - Utilities_destroy +use DAMASK_spectral_Utilities, only: & + Utilities_destroy implicit none call Utilities_destroy() @@ -372,44 +376,4 @@ subroutine basic_destroy() end subroutine basic_destroy - - end module DAMASK_spectral_SolverBasic - - -!-------------------------------------------------------------------------------------------------- -! calculate some additional output - ! if(debugGeneral) then - ! maxCorrectionSkew = 0.0_pReal - ! maxCorrectionSym = 0.0_pReal - ! temp33_Real = 0.0_pReal - ! do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - ! maxCorrectionSym = max(maxCorrectionSym,& - ! maxval(math_symmetric33(field_real(i,j,k,1:3,1:3)))) - ! maxCorrectionSkew = max(maxCorrectionSkew,& - ! maxval(math_skew33(field_real(i,j,k,1:3,1:3)))) - ! temp33_Real = temp33_Real + field_real(i,j,k,1:3,1:3) - ! enddo; enddo; enddo - ! write(6,'(a,1x,es11.4)') 'max symmetric correction of deformation =',& - ! maxCorrectionSym*wgt - ! write(6,'(a,1x,es11.4)') 'max skew correction of deformation =',& - ! maxCorrectionSkew*wgt - ! write(6,'(a,1x,es11.4)') 'max sym/skew of avg correction = ',& - ! maxval(math_symmetric33(temp33_real))/& - ! maxval(math_skew33(temp33_real)) - ! endif - -!-------------------------------------------------------------------------------------------------- -! calculate bounds of det(F) and report - ! if(debugGeneral) then - ! defgradDetMax = -huge(1.0_pReal) - ! defgradDetMin = +huge(1.0_pReal) - ! do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - ! defgradDet = math_det33(F(i,j,k,1:3,1:3)) - ! defgradDetMax = max(defgradDetMax,defgradDet) - ! defgradDetMin = min(defgradDetMin,defgradDet) - ! enddo; enddo; enddo - - ! write(6,'(a,1x,es11.4)') 'max determinant of deformation =', defgradDetMax - ! write(6,'(a,1x,es11.4)') 'min determinant of deformation =', defgradDetMin - ! endif diff --git a/code/DAMASK_spectral_solverBasicPETSc.f90 b/code/DAMASK_spectral_solverBasicPETSc.f90 index 9ff517de8..d54173f7b 100644 --- a/code/DAMASK_spectral_solverBasicPETSc.f90 +++ b/code/DAMASK_spectral_solverBasicPETSc.f90 @@ -37,36 +37,37 @@ module DAMASK_spectral_SolverBasicPETSc !-------------------------------------------------------------------------------------------------- ! PETSc data - DM, private :: da - SNES, private :: snes - Vec, private :: solution_vec - + DM, private :: da + SNES, private :: snes + Vec, private :: solution_vec + !-------------------------------------------------------------------------------------------------- ! common pointwise data - real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, Fdot - real(pReal), private, dimension(:,:,:,:), allocatable :: coordinates - real(pReal), private, dimension(:,:,:), allocatable :: temperature - + real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, Fdot + real(pReal), private, dimension(:,:,:,:), allocatable :: coordinates + real(pReal) :: temperature + !-------------------------------------------------------------------------------------------------- ! stress, stiffness and compliance average etc. - real(pReal), private, dimension(3,3) :: & - F_aim = math_I3, & - F_aim_lastInc = math_I3, & - P_av - character(len=1024), private :: incInfo - real(pReal), private, dimension(3,3,3,3) :: & - C = 0.0_pReal, C_lastInc= 0.0_pReal, & - S = 0.0_pReal - - real(pReal), private :: err_stress, err_div - logical, private :: ForwardData - real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal - + real(pReal), private, dimension(3,3) :: & + F_aim = math_I3, & + F_aim_lastInc = math_I3, & + P_av + character(len=1024), private :: incInfo + real(pReal), private, dimension(3,3,3,3) :: & + C = 0.0_pReal, C_lastInc= 0.0_pReal, & + S = 0.0_pReal + + real(pReal), private :: err_stress, err_div + logical, private :: ForwardData + real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal + + + public :: basicPETSc_init, & + basicPETSc_solution ,& + basicPETSc_destroy + contains - public :: basicPETSc_init, & - basicPETSc_solution ,& - basicPETSc_destroy - contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all neccessary fields and fills them with data, potentially from restart info !-------------------------------------------------------------------------------------------------- @@ -111,17 +112,19 @@ subroutine basicPETSc_init() PetscObject :: dummy call Utilities_init() - write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasicPETSc init -+>>>' write(6,'(a)') ' $Id: DAMASK_spectral_SolverBasicPETSC.f90 1654 2012-08-03 09:25:48Z MPIE\m.diehl $' #include "compilation_info.f90" write(6,'(a)') '' +!-------------------------------------------------------------------------------------------------- +! allocate global fields allocate (F_lastInc (3,3, res(1), res(2),res(3)), source = 0.0_pReal) allocate (Fdot (3,3, res(1), res(2),res(3)), source = 0.0_pReal) allocate (coordinates( res(1), res(2),res(3),3), source = 0.0_pReal) - allocate (temperature( res(1), res(2),res(3)), source = 0.0_pReal) +!-------------------------------------------------------------------------------------------------- +! initialize solver specific parts of PETSc call SNESCreate(PETSC_COMM_WORLD,snes,ierr) CHKERRQ(ierr) call DMDACreate3d(PETSC_COMM_WORLD, & @@ -140,133 +143,130 @@ subroutine basicPETSc_init() call SNESSetFromOptions(snes,ierr) CHKERRQ(ierr) - !-------------------------------------------------------------------------------------------------- - ! init fields - call DMDAVecGetArrayF90(da,solution_vec,F,ierr) - CHKERRQ(ierr) - if (restartInc == 1_pInt) then ! no deformation (no restart) - F_lastInc = spread(spread(spread(math_I3,3,res(1)),4,res(2)),5,res(3)) ! initialize to identity - F = reshape(F_lastInc,[9,res(1),res(2),res(3)]) - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - coordinates(i,j,k,1:3) = geomdim/real(res,pReal)*real([i,j,k],pReal) & - - geomdim/real(2_pInt*res,pReal) - enddo; enddo; enddo - elseif (restartInc > 1_pInt) then ! using old values from file - if (debugRestart) write(6,'(a,i6,a)') 'Reading values of increment ',& - restartInc - 1_pInt,' from file' - flush(6) - call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',& - trim(getSolverJobName()),size(F)) - read (777,rec=1) F - close (777) - call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc',& - trim(getSolverJobName()),size(F_lastInc)) - read (777,rec=1) F_lastInc - close (777) - call IO_read_jobBinaryFile(777,'F_aim',trim(getSolverJobName()),size(F_aim)) - read (777,rec=1) F_aim - close (777) - call IO_read_jobBinaryFile(777,'F_aim_lastInc',trim(getSolverJobName()),size(F_aim_lastInc)) - read (777,rec=1) F_aim_lastInc - close (777) - - coordinates = 0.0 ! change it later!!! - endif +!-------------------------------------------------------------------------------------------------- +! init fields + call DMDAVecGetArrayF90(da,solution_vec,F,ierr) ! get the data out of PETSc to work with + CHKERRQ(ierr) + if (restartInc == 1_pInt) then ! no deformation (no restart) + F_lastInc = spread(spread(spread(math_I3,3,res(1)),4,res(2)),5,res(3)) ! initialize to identity + F = reshape(F_lastInc,[9,res(1),res(2),res(3)]) + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + coordinates(i,j,k,1:3) = geomdim/real(res,pReal)*real([i,j,k],pReal) & + - geomdim/real(2_pInt*res,pReal) + enddo; enddo; enddo + elseif (restartInc > 1_pInt) then ! using old values from file + if (debugRestart) write(6,'(a,i6,a)') 'Reading values of increment ',& + restartInc - 1_pInt,' from file' + flush(6) + call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',& + trim(getSolverJobName()),size(F)) + read (777,rec=1) F + close (777) + call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc',& + trim(getSolverJobName()),size(F_lastInc)) + read (777,rec=1) F_lastInc + close (777) + call IO_read_jobBinaryFile(777,'F_aim',trim(getSolverJobName()),size(F_aim)) + read (777,rec=1) F_aim + close (777) + call IO_read_jobBinaryFile(777,'F_aim_lastInc',trim(getSolverJobName()),size(F_aim_lastInc)) + read (777,rec=1) F_aim_lastInc + close (777) - call Utilities_constitutiveResponse(coordinates,& - reshape(F(0:8,0:res(1)-1_pInt,0:res(2)-1_pInt,0:res(3)-1_pInt),[3,3,res(1),res(2),res(3)]),& - reshape(F(0:8,0:res(1)-1_pInt,0:res(2)-1_pInt,0:res(3)-1_pInt),[3,3,res(1),res(2),res(3)]),& - temperature,0.0_pReal,P,C,P_av,.false.,math_I3) + coordinates = 0.0 ! change it later!!! + endif - call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr) + call Utilities_constitutiveResponse(coordinates,& + reshape(F(0:8,0:res(1)-1_pInt,0:res(2)-1_pInt,0:res(3)-1_pInt),[3,3,res(1),res(2),res(3)]),& + reshape(F(0:8,0:res(1)-1_pInt,0:res(2)-1_pInt,0:res(3)-1_pInt),[3,3,res(1),res(2),res(3)]),& + temperature,0.0_pReal,P,C,P_av,.false.,math_I3) + call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr) ! write data back into PETSc + CHKERRQ(ierr) +!-------------------------------------------------------------------------------------------------- +! reference stiffness and Gamma operator + if (restartInc == 1_pInt) then + call IO_write_jobBinaryFile(777,'C_ref',size(C)) + write (777,rec=1) C + close(777) + elseif (restartInc > 1_pInt) then + call IO_read_jobBinaryFile(777,'C_ref',trim(getSolverJobName()),size(C)) + read (777,rec=1) C + close (777) + endif + call Utilities_updateGamma(C,.True.) - !-------------------------------------------------------------------------------------------------- - ! reference stiffness - if (restartInc == 1_pInt) then - call IO_write_jobBinaryFile(777,'C_ref',size(C)) - write (777,rec=1) C - close(777) - elseif (restartInc > 1_pInt) then - call IO_read_jobBinaryFile(777,'C_ref',trim(getSolverJobName()),size(C)) - read (777,rec=1) C - close (777) - endif - - call Utilities_updateGamma(C,.True.) - - end subroutine basicPETSc_init +end subroutine basicPETSc_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the Basic PETSC scheme with internal iterations !-------------------------------------------------------------------------------------------------- - type(tSolutionState) function & - basicPETSc_solution(incInfoIn,guess,timeinc,timeinc_old,P_BC,F_BC,temperature_bc,rotation_BC) - use numerics, only: & - update_gamma - use math, only: & - math_mul33x33 ,& - math_rotate_backward33 - use mesh, only: & - res,& - geomdim,& - deformed_fft - use IO, only: & - IO_write_JobBinaryFile - use DAMASK_spectral_Utilities, only: & - tBoundaryCondition, & - Utilities_calculateRate, & - Utilities_forwardField, & - Utilities_maskedCompliance, & - Utilities_updateGamma, & - cutBack - use FEsolving, only: & - restartWrite, & - terminallyIll - implicit none +type(tSolutionState) function & + basicPETSc_solution(incInfoIn,guess,timeinc,timeinc_old,P_BC,F_BC,temperature_bc,rotation_BC) + use numerics, only: & + update_gamma + use math, only: & + math_mul33x33 ,& + math_rotate_backward33 + use mesh, only: & + res,& + geomdim,& + deformed_fft + use IO, only: & + IO_write_JobBinaryFile + use DAMASK_spectral_Utilities, only: & + tBoundaryCondition, & + Utilities_calculateRate, & + Utilities_forwardField, & + Utilities_maskedCompliance, & + Utilities_updateGamma, & + cutBack + use FEsolving, only: & + restartWrite, & + terminallyIll + implicit none #include #include !-------------------------------------------------------------------------------------------------- ! input data for solution - real(pReal), intent(in) :: timeinc, timeinc_old, temperature_bc - logical, intent(in):: guess - type(tBoundaryCondition), intent(in) :: P_BC,F_BC - real(pReal), dimension(3,3), intent(in) :: rotation_BC + real(pReal), intent(in) :: timeinc, timeinc_old, temperature_bc + logical, intent(in):: guess + type(tBoundaryCondition), intent(in) :: P_BC,F_BC + real(pReal), dimension(3,3), intent(in) :: rotation_BC character(len=*), intent(in) :: incInfoIn - real(pReal), dimension(3,3),save :: F_aimDot=0.0_pReal - real(pReal), dimension(3,3) :: & - F_aim_lab + real(pReal), dimension(3,3),save :: F_aimDot=0.0_pReal + real(pReal), dimension(3,3) :: & + F_aim_lab !-------------------------------------------------------------------------------------------------- ! loop variables, convergence etc. - real(pReal), dimension(3,3) :: temp33_Real + real(pReal), dimension(3,3) :: temp33_Real !-------------------------------------------------------------------------------------------------- ! - PetscScalar, pointer :: F(:,:,:,:) - PetscErrorCode :: ierr - SNESConvergedReason :: reason + PetscScalar, pointer :: F(:,:,:,:) + PetscErrorCode :: ierr + SNESConvergedReason :: reason !-------------------------------------------------------------------------------------------------- ! restart information for spectral solver - incInfo = incInfoIn - if (restartWrite) then - write(6,'(a)') 'writing converged results for restart' - call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(F_lastInc)) - write (777,rec=1) F_LastInc - close (777) - call IO_write_jobBinaryFile(777,'C',size(C)) - write (777,rec=1) C - close(777) - endif - call DMDAVecGetArrayF90(da,solution_vec,F,ierr) + incInfo = incInfoIn + if (restartWrite) then + write(6,'(a)') 'writing converged results for restart' + call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(F_lastInc)) + write (777,rec=1) F_LastInc + close (777) + call IO_write_jobBinaryFile(777,'C',size(C)) + write (777,rec=1) C + close(777) + endif + call DMDAVecGetArrayF90(da,solution_vec,F,ierr) -if ( cutBack) then - F_aim = F_aim_lastInc - F = reshape(F_lastInc,[9,res(1),res(2),res(3)]) - C = C_lastInc -else - - C_lastInc = C + if ( cutBack) then + F_aim = F_aim_lastInc + F = reshape(F_lastInc,[9,res(1),res(2),res(3)]) + C = C_lastInc + else + C_lastInc = C !-------------------------------------------------------------------------------------------------- ! calculate rate for aim @@ -277,7 +277,7 @@ else endif if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old F_aim_lastInc = F_aim - + !-------------------------------------------------------------------------------------------------- ! update coordinates and rate and forward last inc call deformed_fft(res,geomdim,math_rotate_backward33(F_aim_lastInc,rotation_BC), & @@ -291,167 +291,177 @@ else F = reshape(Utilities_forwardField(timeinc,F_aim,F_lastInc,Fdot),[9,res(1),res(2),res(3)]) call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr) + CHKERRQ(ierr) call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,rotation_BC),1.0_pReal,F_lastInc,coordinates) !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C) - if (update_gamma) call Utilities_updateGamma(C,restartWrite) - - ForwardData = .True. - mask_stress = P_BC%maskFloat - params%P_BC = P_BC%values - params%rotation_BC = rotation_BC - params%timeinc = timeinc + S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C) + if (update_gamma) call Utilities_updateGamma(C,restartWrite) + + ForwardData = .True. + mask_stress = P_BC%maskFloat + params%P_BC = P_BC%values + params%rotation_BC = rotation_BC + params%timeinc = timeinc - call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr) - call SNESGetConvergedReason(snes,reason,ierr) - basicPETSc_solution%termIll = terminallyIll - terminallyIll = .false. - BasicPETSC_solution%converged =.false. - if (reason > 0 ) then - BasicPETSC_solution%converged = .true. - endif + call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr) + CHKERRQ(ierr) + call SNESGetConvergedReason(snes,reason,ierr) + CHKERRQ(ierr) + basicPETSc_solution%termIll = terminallyIll + terminallyIll = .false. + BasicPETSC_solution%converged =.false. + if (reason > 0 ) then + BasicPETSC_solution%converged = .true. + endif - end function BasicPETSc_solution +end function BasicPETSc_solution !-------------------------------------------------------------------------------------------------- !> @brief forms the AL residual vector !-------------------------------------------------------------------------------------------------- - subroutine BasicPETSC_formResidual(myIn,x_scal,f_scal,dummy,ierr) - - use numerics, only: & - itmax, & - itmin - use math, only: & - math_rotate_backward33, & - math_transpose33, & - math_mul3333xx33 - use mesh, only: & - res, & - wgt - use DAMASK_spectral_Utilities, only: & - field_real, & - Utilities_FFTforward, & - Utilities_FFTbackward, & - Utilities_fourierConvolution, & - Utilities_constitutiveResponse, & - Utilities_divergenceRMS - use IO, only : IO_intOut - implicit none +subroutine BasicPETSC_formResidual(myIn,x_scal,f_scal,dummy,ierr) - real(pReal), dimension(3,3) :: F_aim_lab_lastIter, F_aim_lab - - DMDALocalInfo, dimension(*) :: myIn - PetscScalar, dimension(3,3,res(1),res(2),res(3)) :: x_scal - PetscScalar, dimension(3,3,res(1),res(2),res(3)):: f_scal - PetscInt :: iter, nfuncs - PetscObject :: dummy - PetscErrorCode :: ierr - - call SNESGetNumberFunctionEvals(snes,nfuncs,ierr) - call SNESGetIterationNumber(snes,iter,ierr) + use numerics, only: & + itmax, & + itmin + use math, only: & + math_rotate_backward33, & + math_transpose33, & + math_mul3333xx33 + use mesh, only: & + res, & + wgt + use DAMASK_spectral_Utilities, only: & + field_real, & + Utilities_FFTforward, & + Utilities_FFTbackward, & + Utilities_fourierConvolution, & + Utilities_constitutiveResponse, & + Utilities_divergenceRMS + use IO, only : IO_intOut + + implicit none + real(pReal), dimension(3,3) :: F_aim_lab_lastIter, F_aim_lab + + DMDALocalInfo, dimension(*) :: myIn + PetscScalar, dimension(3,3,res(1),res(2),res(3)) :: x_scal + PetscScalar, dimension(3,3,res(1),res(2),res(3)):: f_scal + PetscInt :: iter, nfuncs + PetscObject :: dummy + PetscErrorCode :: ierr + + call SNESGetNumberFunctionEvals(snes,nfuncs,ierr) + CHKERRQ(ierr) + call SNESGetIterationNumber(snes,iter,ierr) + CHKERRQ(ierr) +!-------------------------------------------------------------------------------------------------- +! report begin of new iteration + write(6,'(/,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & + ' @ Iter. ', itmin, '<',iter, '≤', itmax + write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim =',& + math_transpose33(F_aim) + F_aim_lab_lastIter = math_rotate_backward33(F_aim,params%rotation_BC) + +!-------------------------------------------------------------------------------------------------- +! evaluate constitutive response + call Utilities_constitutiveResponse(coordinates,F_lastInc,x_scal,temperature,params%timeinc, & + f_scal,C,P_av,ForwardData,params%rotation_BC) + ForwardData = .false. - !-------------------------------------------------------------------------------------------------- - ! report begin of new iteration - write(6,'(/,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & - ' @ Iter. ', itmin, '<',iter, '≤', itmax - write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim =',& - math_transpose33(F_aim) - F_aim_lab_lastIter = math_rotate_backward33(F_aim,params%rotation_BC) - - !-------------------------------------------------------------------------------------------------- - ! evaluate constitutive response - call Utilities_constitutiveResponse(coordinates,F_lastInc,x_scal,temperature,params%timeinc, & - f_scal,C,P_av,ForwardData,params%rotation_BC) - ForwardData = .false. - !-------------------------------------------------------------------------------------------------- ! 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_lab = math_rotate_backward33(F_aim,params%rotation_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 + F_aim_lab = math_rotate_backward33(F_aim,params%rotation_BC) + !-------------------------------------------------------------------------------------------------- ! updated deformation gradient using fix point algorithm of basic scheme - field_real = 0.0_pReal - field_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = reshape(f_scal,[res(1),res(2),res(3),3,3],& - order=[4,5,1,2,3]) ! field real has a different order - call Utilities_FFTforward() - err_div = Utilities_divergenceRMS() - call Utilities_fourierConvolution(F_aim_lab_lastIter - F_aim_lab) - call Utilities_FFTbackward() - - !-------------------------------------------------------------------------------------------------- - ! constructing residual - f_scal = reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),shape(x_scal),order=[3,4,5,1,2]) - write(6,'(/,a)') '==========================================================================' - end subroutine BasicPETSc_formResidual + field_real = 0.0_pReal + field_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = reshape(f_scal,[res(1),res(2),res(3),3,3],& + order=[4,5,1,2,3]) ! field real has a different order + call Utilities_FFTforward() + err_div = Utilities_divergenceRMS() + call Utilities_fourierConvolution(F_aim_lab_lastIter - F_aim_lab) + call Utilities_FFTbackward() + +!-------------------------------------------------------------------------------------------------- +! constructing residual + f_scal = reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),shape(x_scal),order=[3,4,5,1,2]) + write(6,'(/,a)') '==========================================================================' +end subroutine BasicPETSc_formResidual + !-------------------------------------------------------------------------------------------------- !> @brief convergence check !-------------------------------------------------------------------------------------------------- - subroutine BasicPETSc_converged(snes_local,it,xnorm,snorm,fnorm,reason,dummy,ierr) - - use numerics, only: & - itmax, & - itmin, & - err_div_tol, & - err_stress_tolrel, & - err_stress_tolabs - - use math, only: & - math_mul33x33, & - math_eigenvalues33, & - math_transpose33 - - implicit none +subroutine BasicPETSc_converged(snes_local,it,xnorm,snorm,fnorm,reason,dummy,ierr) - SNES :: snes_local - PetscInt :: it - PetscReal :: xnorm, snorm, fnorm - SNESConvergedReason :: reason - PetscObject :: dummy - PetscErrorCode :: ierr - logical :: Converged - real(pReal) :: pAvgDivL2 - - pAvgDivL2 = sqrt(maxval(math_eigenvalues33(math_mul33x33(P_av,math_transpose33(P_av))))) - Converged = (it > itmin .and. & - all([ err_div/pAvgDivL2/err_div_tol, & - err_stress/min(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs)] < 1.0_pReal)) - - if (Converged) then - reason = 1 - elseif (it > itmax) then - reason = -1 - else - reason = 0 - endif - - write(6,'(a,f6.2,a,es11.4,a)') 'error divergence = ', err_div/pAvgDivL2/err_div_tol,& - ' (',err_div/pAvgDivL2,' N/m³)' - write(6,'(a,f6.2,a,es11.4,a)') 'error stress = ', err_stress/min(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs), & - ' (',err_stress,' Pa)' + use numerics, only: & + itmax, & + itmin, & + err_div_tol, & + err_stress_tolrel, & + err_stress_tolabs + + use math, only: & + math_mul33x33, & + math_eigenvalues33, & + math_transpose33 + + implicit none + + SNES :: snes_local + PetscInt :: it + PetscReal :: xnorm, snorm, fnorm + SNESConvergedReason :: reason + PetscObject :: dummy + PetscErrorCode :: ierr + logical :: Converged + real(pReal) :: pAvgDivL2 + + pAvgDivL2 = sqrt(maxval(math_eigenvalues33(math_mul33x33(P_av,math_transpose33(P_av))))) + Converged = (it > itmin .and. & + all([ err_div/pAvgDivL2/err_div_tol, & + err_stress/min(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs)] < 1.0_pReal)) + + if (Converged) then + reason = 1 + elseif (it > itmax) then + reason = -1 + else + reason = 0 + endif + + write(6,'(a,f6.2,a,es11.4,a)') 'error divergence = ', err_div/pAvgDivL2/err_div_tol,& + ' (',err_div/pAvgDivL2,' N/m³)' + write(6,'(a,f6.2,a,es11.4,a)') 'error stress = ', err_stress/min(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs), & + ' (',err_stress,' Pa)' +end subroutine BasicPETSc_converged - end subroutine BasicPETSc_converged !-------------------------------------------------------------------------------------------------- !> @brief destroy routine !-------------------------------------------------------------------------------------------------- - subroutine BasicPETSc_destroy() +subroutine BasicPETSc_destroy() - use DAMASK_spectral_Utilities, only: & - Utilities_destroy - implicit none - PetscErrorCode :: ierr + use DAMASK_spectral_Utilities, only: & + Utilities_destroy + + implicit none + PetscErrorCode :: ierr - call VecDestroy(solution_vec,ierr) - call SNESDestroy(snes,ierr) - call DMDestroy(da,ierr) - call PetscFinalize(ierr) - call Utilities_destroy() + call VecDestroy(solution_vec,ierr) + CHKERRQ(ierr) + call SNESDestroy(snes,ierr) + CHKERRQ(ierr) + call DMDestroy(da,ierr) + CHKERRQ(ierr) + call PetscFinalize(ierr) + CHKERRQ(ierr) + call Utilities_destroy() + CHKERRQ(ierr) +end subroutine BasicPETSc_destroy - end subroutine BasicPETSc_destroy - - end module DAMASK_spectral_SolverBasicPETSc +end module DAMASK_spectral_SolverBasicPETSc diff --git a/code/DAMASK_spectral_utilities.f90 b/code/DAMASK_spectral_utilities.f90 index c9fcbc880..00f736331 100644 --- a/code/DAMASK_spectral_utilities.f90 +++ b/code/DAMASK_spectral_utilities.f90 @@ -17,26 +17,29 @@ module DAMASK_spectral_utilities !-------------------------------------------------------------------------------------------------- ! variables storing information for spectral method and FFTW - real(pReal), public, dimension(:,:,:,:,:), pointer :: field_real !< real representation (some stress of deformation) of field_fourier - type(C_PTR), private :: plan_forward, plan_backward !< plans for FFTW + real(pReal), public, dimension(:,:,:,:,:), pointer :: field_real !< real representation (some stress or deformation) of field_fourier + complex(pReal),private, dimension(:,:,:,:,:), pointer :: field_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 - complex(pReal),private, dimension(:,:,:,:,:), pointer :: field_fourier !< field on which the Fourier transform operates real(pReal), private, dimension(3,3,3,3) :: C_ref !< reference stiffness !-------------------------------------------------------------------------------------------------- ! debug fftw - type(C_PTR), private :: plan_scalarField_forth, plan_scalarField_back !< plans for FFTW in case of debugging the Fourier transform - complex(pReal),private, dimension(:,:,:), pointer :: scalarField_real !< scalar field real representation for debug of FFTW - complex(pReal),private, dimension(:,:,:), pointer :: scalarField_fourier !< scalar field complex representation for debug of FFTW + complex(pReal),private, dimension(:,:,:), pointer :: scalarField_real, & !< scalar field real representation for debug of FFTW + scalarField_fourier !< scalar field complex representation for debug of FFTW !-------------------------------------------------------------------------------------------------- ! debug divergence - type(C_PTR), private :: plan_divergence !< plan for FFTW in case of debugging divergence calculation 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(:,:,:,:), allocatable :: divergence_post !< data of divergence calculation using function from core modules (serves as a reference) +!-------------------------------------------------------------------------------------------------- +! 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 + !-------------------------------------------------------------------------------------------------- ! variables controlling debugging logical, public :: & @@ -80,8 +83,8 @@ module DAMASK_spectral_utilities contains !-------------------------------------------------------------------------------------------------- -!> @brief allocates all neccessary fields, sets debug flags, create plans for fftw -!> @details Sets the debug levels for general, divergence, restart and fftw from the biwise coding +!> @brief allocates all neccessary fields, sets debug flags, create plans for FFTW +!> @details Sets the debug levels for general, divergence, restart and FFTW from the biwise coding !> provided by the debug module to logicals. !> Allocates all fields used by FFTW and create the corresponding plans depending on the debug !> level chosen. @@ -132,8 +135,8 @@ subroutine utilities_init() !-------------------------------------------------------------------------------------------------- ! allocation - allocate (xi (3,res1_red,res(2),res(3)), source = 0.0_pReal) ! start out isothermally - tensorField = fftw_alloc_complex(int(res1_red*res(2)*res(3)*9_pInt,C_SIZE_T)) ! allocate continous data using a C function, C_SIZE_T is of type integer(8) + 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 @@ -148,7 +151,7 @@ subroutine utilities_init() call fftw_set_timelimit(fftw_timelimit) ! set timelimit for plan creation !-------------------------------------------------------------------------------------------------- -! creating plans +! 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 @@ -162,7 +165,7 @@ subroutine utilities_init() 1, res(3)*res(2)*(res(1)+2_pInt), fftw_planner_flag) ! striding, product of physical length in the 3 dimensions, planner precision !-------------------------------------------------------------------------------------------------- -! depending on (debug) options, allocate more memory and create additional plans +! 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]) @@ -267,7 +270,7 @@ end subroutine utilities_updateGamma !-------------------------------------------------------------------------------------------------- !> @brief forward FFT of data in field_real to field_fourier with highest freqs. removed -!> Does an unweighted FFT transform from real to complex. +!> @detailed 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 !-------------------------------------------------------------------------------------------------- @@ -326,7 +329,7 @@ end subroutine utilities_FFTforward !-------------------------------------------------------------------------------------------------- !> @brief backward FFT of data in field_fourier to field_real -!> Does an inverse FFT transform from complex to real +!> @detailed 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 @@ -343,7 +346,7 @@ subroutine utilities_FFTbackward(row,column) integer(pInt) :: i, j, k, m, n !-------------------------------------------------------------------------------------------------- -! unpack FFT data for conj complex symmetric part +! unpack FFT data for conj complex symmetric part. This data is not transformed when using c2r if (debugFFTW) 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 @@ -377,6 +380,29 @@ subroutine utilities_FFTbackward(row,column) field_real = field_real * wgt ! normalize the result by number of elements +!-------------------------------------------------------------------------------------------------- +! calculate some additional output + ! if(debugGeneral) then + ! maxCorrectionSkew = 0.0_pReal + ! maxCorrectionSym = 0.0_pReal + ! temp33_Real = 0.0_pReal + ! do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + ! maxCorrectionSym = max(maxCorrectionSym,& + ! maxval(math_symmetric33(field_real(i,j,k,1:3,1:3)))) + ! maxCorrectionSkew = max(maxCorrectionSkew,& + ! maxval(math_skew33(field_real(i,j,k,1:3,1:3)))) + ! temp33_Real = temp33_Real + field_real(i,j,k,1:3,1:3) + ! enddo; enddo; enddo + ! write(6,'(a,1x,es11.4)') 'max symmetric correction of deformation =',& + ! maxCorrectionSym*wgt + ! write(6,'(a,1x,es11.4)') 'max skew correction of deformation =',& + ! maxCorrectionSkew*wgt + ! write(6,'(a,1x,es11.4)') 'max sym/skew of avg correction = ',& + ! maxval(math_symmetric33(temp33_real))/& + ! maxval(math_skew33(temp33_real)) + ! endif + + end subroutine utilities_FFTbackward @@ -397,7 +423,7 @@ subroutine utilities_fourierConvolution(fieldAim) real(pReal), intent(in), dimension(3,3) :: fieldAim !< desired average value of the field after convolution real(pReal), dimension(3,3) :: xiDyad, temp33_Real real(pReal) :: filter !< weighting of current component - complex(pReal), dimension(3,3) :: temp33_complex + complex(pReal), dimension(3,3) :: temp33_complex integer(pInt) :: & i, j, k, & l, m, n, o @@ -621,7 +647,7 @@ subroutine utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,ti CPFEM_general implicit none - real(pReal), intent(inout), dimension(res(1),res(2),res(3)) :: temperature !< temperature field + real(pReal), intent(inout) :: temperature !< temperature (no field) real(pReal), intent(in), dimension(res(1),res(2),res(3),3) :: coordinates !< coordinates field real(pReal), intent(in), dimension(3,3,res(1),res(2),res(3)) :: & F_lastInc, & !< target deformation gradient @@ -655,7 +681,20 @@ subroutine utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,ti calcMode = 2_pInt collectMode = 5_pInt endif - +!-------------------------------------------------------------------------------------------------- +! calculate bounds of det(F) and report + ! if(debugGeneral) then + ! defgradDetMax = -huge(1.0_pReal) + ! defgradDetMin = +huge(1.0_pReal) + ! do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + ! defgradDet = math_det33(F(i,j,k,1:3,1:3)) + ! defgradDetMax = max(defgradDetMax,defgradDet) + ! defgradDetMin = min(defgradDetMin,defgradDet) + ! enddo; enddo; enddo + + ! write(6,'(a,1x,es11.4)') 'max determinant of deformation =', defgradDetMax + ! write(6,'(a,1x,es11.4)') 'min determinant of deformation =', defgradDetMin + ! endif if (DebugGeneral) write(6,*) 'collect mode: ', collectMode,' calc mode: ', calcMode flush(6) @@ -664,7 +703,7 @@ subroutine utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,ti ielem = ielem + 1_pInt call CPFEM_general(collectMode,& ! collect cycle coordinates(i,j,k,1:3), F_lastInc(1:3,1:3,i,j,k),F(1:3,1:3,i,j,k), & - temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde,P(1:3,1:3,i,j,k),dPdF) + temperature,timeinc,ielem,1_pInt,sigma,dsde,P(1:3,1:3,i,j,k),dPdF) collectMode = 3_pInt enddo; enddo; enddo @@ -676,7 +715,7 @@ subroutine utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,ti ielem = ielem + 1_pInt call CPFEM_general(calcMode,& ! first element in first iteration retains CPFEM_mode 1, coordinates(i,j,k,1:3),F_lastInc(1:3,1:3,i,j,k), F(1:3,1:3,i,j,k), & ! others get 2 (saves winding forward effort) - temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde,P(1:3,1:3,i,j,k),dPdF) + temperature,timeinc,ielem,1_pInt,sigma,dsde,P(1:3,1:3,i,j,k),dPdF) calcMode = 2_pInt C = C + dPdF enddo; enddo; enddo @@ -711,7 +750,7 @@ pure function utilities_calculateRate(delta_aim,timeinc,timeinc_old,guess,field_ field !< data of current step real(pReal), dimension(3,3,res(1),res(2),res(3)) :: utilities_calculateRate - if (guess) then + if(guess) then utilities_calculateRate = (field-field_lastInc) / timeinc_old else utilities_calculateRate = spread(spread(spread(delta_aim,3,res(1)),4,res(2)),5,res(3))/timeinc