From 70c4e11742f2b430d036370d37ad0f07741de6d5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 12 Nov 2012 14:14:39 +0000 Subject: [PATCH] added comments and structured the code. Temperature is not longer stored for each point but only to simplify restart behavior as long as it is not fully supported by the constitutive models --- code/DAMASK_spectral_solverAL.f90 | 740 +++++++++++----------- code/DAMASK_spectral_solverBasic.f90 | 156 ++--- code/DAMASK_spectral_solverBasicPETSc.f90 | 546 ++++++++-------- code/DAMASK_spectral_utilities.f90 | 83 ++- 4 files changed, 774 insertions(+), 751 deletions(-) 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