2013-07-24 18:36:16 +05:30
! Copyright 2011-13 Max-Planck-Institut für Eisenforschung GmbH
!
! This file is part of DAMASK,
! the Düsseldorf Advanced MAterial Simulation Kit.
!
! DAMASK is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! DAMASK is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with DAMASK. If not, see <http://www.gnu.org/licenses/>.
!
!--------------------------------------------------------------------------------------------------
! $Id$
!--------------------------------------------------------------------------------------------------
!> @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 Polarisation scheme solver
!--------------------------------------------------------------------------------------------------
module DAMASK_spectral_solverPolarisation
use prec , only : &
pInt , &
pReal
use math , only : &
math_I3
use DAMASK_spectral_utilities , only : &
tSolutionState
implicit none
private
#include <finclude/petscsys.h>
#include <finclude/petscdmda.h>
#include <finclude/petscsnes.h>
character ( len = * ) , parameter , public :: &
DAMASK_spectral_solverPolarisation_label = 'polarisation'
!--------------------------------------------------------------------------------------------------
! derived types
type tSolutionParams !< @todo use here the type definition for a full loadcase including mask
real ( pReal ) , dimension ( 3 , 3 ) :: P_BC , rotation_BC
real ( pReal ) :: timeinc
2013-07-26 21:55:37 +05:30
real ( pReal ) :: timeincOld
2013-07-24 18:36:16 +05:30
real ( pReal ) :: temperature
2013-07-26 21:55:37 +05:30
real ( pReal ) :: density
2013-07-24 18:36:16 +05:30
end type tSolutionParams
type ( tSolutionParams ) , private :: params
real ( pReal ) , private , dimension ( 3 , 3 ) :: mask_stress = 0.0_pReal
!--------------------------------------------------------------------------------------------------
! PETSc data
DM , private :: da
SNES , private :: snes
Vec , private :: solution_vec
!--------------------------------------------------------------------------------------------------
! common pointwise data
real ( pReal ) , private , dimension ( : , : , : , : , : ) , allocatable :: &
F_lastInc , & !< field of previous compatible deformation gradients
2013-07-26 21:55:37 +05:30
F_lastInc2 , & !< field of 2nd previous compatible deformation gradients
2013-07-24 18:36:16 +05:30
F_tau_lastInc , & !< field of previous incompatible deformation gradient
Fdot , & !< field of assumed rate of compatible deformation gradient
F_tauDot !< field of assumed rate of incopatible deformation gradient
2013-07-26 21:55:37 +05:30
complex ( pReal ) , private , dimension ( : , : , : , : , : ) , allocatable :: inertiaField_fourier
2013-07-24 18:36:16 +05:30
!--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc.
real ( pReal ) , private , dimension ( 3 , 3 ) :: &
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
2013-08-09 21:55:13 +05:30
F_av = 0.0_pReal , & !< average incompatible def grad field
2013-07-30 21:02:55 +05:30
P_av = 0.0_pReal , & !< average 1st Piola--Kirchhoff stress
P_avLastEval = 0.0_pReal !< average 1st Piola--Kirchhoff stress last call of CPFEM_general
2013-08-07 22:50:05 +05:30
character ( len = 1024 ) , private :: incInfo !< time and increment information
2013-07-24 18:36:16 +05:30
real ( pReal ) , private , dimension ( 3 , 3 , 3 , 3 ) :: &
C_volAvg = 0.0_pReal , & !< current volume average stiffness
C_volAvgLastInc = 0.0_pReal , & !< previous volume average stiffness
C_minMaxAvg = 0.0_pReal , & !< current (min+max)/2 stiffness
S = 0.0_pReal , & !< current compliance (filled up with zeros)
C_scale = 0.0_pReal , &
S_scale = 0.0_pReal
real ( pReal ) , private :: &
2013-08-09 21:55:13 +05:30
err_BC , & !< deviation from stress BC
2013-08-07 22:50:05 +05:30
err_curl , & !< RMS of curl of F
err_div !< RMS of div of P
2013-07-24 18:36:16 +05:30
logical , private :: ForwardData
integer ( pInt ) , private :: &
totalIter = 0_pInt !< total iteration in current increment
public :: &
Polarisation_init , &
Polarisation_solution , &
2013-12-18 15:05:05 +05:30
Polarisation_forward , &
2013-07-24 18:36:16 +05:30
Polarisation_destroy
external :: &
VecDestroy , &
DMDestroy , &
DMDACreate3D , &
DMCreateGlobalVector , &
2013-12-17 21:07:14 +05:30
DMDASNESSetFunctionLocal , &
2013-07-24 18:36:16 +05:30
PETScFinalize , &
SNESDestroy , &
SNESGetNumberFunctionEvals , &
SNESGetIterationNumber , &
SNESSolve , &
SNESSetDM , &
SNESGetConvergedReason , &
SNESSetConvergenceTest , &
SNESSetFromOptions , &
SNESCreate , &
MPI_Abort
contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields and fills them with data, potentially from restart info
!> @todo use sourced allocation, e.g. allocate(Fdot,source = F_lastInc)
!--------------------------------------------------------------------------------------------------
subroutine Polarisation_init ( temperature )
use , intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment)
use IO , only : &
IO_intOut , &
2013-09-18 19:37:55 +05:30
IO_read_realFile , &
2013-07-24 18:36:16 +05:30
IO_timeStamp
use debug , only : &
debug_level , &
debug_spectral , &
debug_spectralRestart
use FEsolving , only : &
restartInc
use DAMASK_interface , only : &
getSolverJobName
use DAMASK_spectral_Utilities , only : &
Utilities_init , &
Utilities_constitutiveResponse , &
Utilities_updateGamma , &
grid , &
2013-07-26 21:55:37 +05:30
grid1Red , &
2013-07-24 18:36:16 +05:30
geomSize , &
wgt
use mesh , only : &
mesh_ipCoordinates , &
mesh_deformedCoordsFFT
use math , only : &
math_invSym3333
implicit none
real ( pReal ) , intent ( inout ) :: &
temperature
#include <finclude/petscdmda.h90>
#include <finclude/petscsnes.h90>
2013-12-17 21:07:14 +05:30
#include <finclude/petscvec.h>
2013-07-24 18:36:16 +05:30
real ( pReal ) , dimension ( : , : , : , : , : ) , allocatable :: P
real ( pReal ) , dimension ( 3 , 3 ) :: &
temp33_Real = 0.0_pReal
PetscErrorCode :: ierr
PetscObject :: dummy
PetscScalar , pointer , dimension ( : , : , : , : ) :: xx_psc , F , F_tau
call Utilities_init ( )
write ( 6 , '(/,a)' ) ' <<<+- DAMASK_spectral_solverPolarisation init -+>>>'
write ( 6 , '(a)' ) ' $Id$'
2014-03-26 14:39:28 +05:30
write ( 6 , '(a15,a)' ) ' Current time: ' , IO_timeStamp ( )
2013-07-24 18:36:16 +05:30
#include "compilation_info.f90"
allocate ( P ( 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ) , source = 0.0_pReal )
!--------------------------------------------------------------------------------------------------
! allocate global fields
allocate ( F_lastInc ( 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ) , source = 0.0_pReal )
2013-07-26 21:55:37 +05:30
allocate ( F_lastInc2 ( 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ) , source = 0.0_pReal )
2013-07-24 18:36:16 +05:30
allocate ( Fdot ( 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ) , source = 0.0_pReal )
allocate ( F_tau_lastInc ( 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ) , source = 0.0_pReal )
allocate ( F_tauDot ( 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ) , source = 0.0_pReal )
2013-07-26 21:55:37 +05:30
allocate ( inertiaField_fourier ( grid1Red , grid ( 2 ) , grid ( 3 ) , 3 , 3 ) , source = cmplx ( 0.0_pReal , 0.0_pReal , pReal ) )
2013-07-24 18:36:16 +05:30
!--------------------------------------------------------------------------------------------------
! 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 , grid ( 1 ) , grid ( 2 ) , grid ( 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 )
2013-12-17 21:07:14 +05:30
call DMDASNESSetFunctionLocal ( da , INSERT_VALUES , Polarisation_formResidual , dummy , ierr )
CHKERRQ ( ierr )
2013-07-24 18:36:16 +05:30
call SNESSetDM ( snes , da , ierr ) ; CHKERRQ ( ierr )
call SNESSetConvergenceTest ( snes , Polarisation_converged , dummy , PETSC_NULL_FUNCTION , ierr )
CHKERRQ ( ierr )
call SNESSetFromOptions ( snes , ierr ) ; CHKERRQ ( ierr )
!--------------------------------------------------------------------------------------------------
! init fields
call DMDAVecGetArrayF90 ( da , solution_vec , xx_psc , ierr ) ; CHKERRQ ( ierr ) ! places pointer xx_psc on PETSc data
F = > xx_psc ( 0 : 8 , : , : , : )
F_tau = > xx_psc ( 9 : 17 , : , : , : )
if ( restartInc == 1_pInt ) then ! no deformation (no restart)
F_lastInc = spread ( spread ( spread ( math_I3 , 3 , grid ( 1 ) ) , 4 , grid ( 2 ) ) , 5 , grid ( 3 ) ) ! initialize to identity
2013-07-26 21:55:37 +05:30
F_lastInc2 = F_lastInc
2013-07-24 18:36:16 +05:30
F = reshape ( F_lastInc , [ 9 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ] )
2013-07-30 21:02:55 +05:30
F_tau = 2.0_pReal * F
2014-03-25 21:14:16 +05:30
F_tau_lastInc = 2.0_pReal * F_lastInc
2013-07-24 18:36:16 +05:30
elseif ( restartInc > 1_pInt ) then
if ( iand ( debug_level ( debug_spectral ) , debug_spectralRestart ) / = 0 ) &
write ( 6 , '(/,a,' / / IO_intOut ( restartInc - 1_pInt ) / / ',a)' ) &
'reading values of increment' , restartInc - 1_pInt , 'from file'
flush ( 6 )
2013-09-20 19:52:37 +05:30
call IO_read_realFile ( 777 , 'F' , trim ( getSolverJobName ( ) ) , size ( F ) )
2013-07-24 18:36:16 +05:30
read ( 777 , rec = 1 ) F
close ( 777 )
2013-09-20 19:52:37 +05:30
call IO_read_realFile ( 777 , 'F_lastInc' , trim ( getSolverJobName ( ) ) , size ( F_lastInc ) )
2013-07-24 18:36:16 +05:30
read ( 777 , rec = 1 ) F_lastInc
close ( 777 )
2013-09-20 19:52:37 +05:30
call IO_read_realFile ( 777 , 'F_lastInc2' , trim ( getSolverJobName ( ) ) , size ( F_lastInc2 ) )
2013-07-26 21:55:37 +05:30
read ( 777 , rec = 1 ) F_lastInc2
close ( 777 )
2013-09-20 19:52:37 +05:30
call IO_read_realFile ( 777 , 'F_tau' , trim ( getSolverJobName ( ) ) , size ( F_tau ) )
2013-07-24 18:36:16 +05:30
read ( 777 , rec = 1 ) F_tau
close ( 777 )
2013-09-18 19:37:55 +05:30
call IO_read_realFile ( 777 , 'F_tau_lastInc' , &
2013-07-24 18:36:16 +05:30
trim ( getSolverJobName ( ) ) , size ( F_tau_lastInc ) )
read ( 777 , rec = 1 ) F_tau_lastInc
close ( 777 )
2014-03-24 15:29:30 +05:30
call IO_read_realFile ( 777 , 'F_aim' , trim ( getSolverJobName ( ) ) , size ( F_aim ) )
read ( 777 , rec = 1 ) F_aim
close ( 777 )
call IO_read_realFile ( 777 , 'F_aim_lastInc' , trim ( getSolverJobName ( ) ) , size ( F_aim_lastInc ) )
read ( 777 , rec = 1 ) F_aim_lastInc
close ( 777 )
2013-09-18 19:37:55 +05:30
call IO_read_realFile ( 777 , 'F_aimDot' , trim ( getSolverJobName ( ) ) , size ( f_aimDot ) )
2013-07-24 18:36:16 +05:30
read ( 777 , rec = 1 ) f_aimDot
close ( 777 )
2014-03-25 21:14:16 +05:30
endif
mesh_ipCoordinates = reshape ( mesh_deformedCoordsFFT ( geomSize , reshape ( &
F , [ 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ] ) ) , [ 3 , 1 , product ( grid ) ] )
call Utilities_constitutiveResponse ( F_lastInc , reshape ( F , [ 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ] ) , &
temperature , 0.0_pReal , P , C_volAvg , C_minMaxAvg , temp33_Real , . false . , math_I3 )
nullify ( F )
nullify ( F_tau )
call DMDAVecRestoreArrayF90 ( da , solution_vec , xx_psc , ierr ) ; CHKERRQ ( ierr ) ! write data back to PETSc
if ( restartInc > 1_pInt ) then ! using old values from files
if ( iand ( debug_level ( debug_spectral ) , debug_spectralRestart ) / = 0 ) &
write ( 6 , '(/,a,' / / IO_intOut ( restartInc - 1_pInt ) / / ',a)' ) &
'reading more values of increment' , restartInc - 1_pInt , 'from file'
flush ( 6 )
2013-09-18 19:37:55 +05:30
call IO_read_realFile ( 777 , 'C_volAvg' , trim ( getSolverJobName ( ) ) , size ( C_volAvg ) )
2013-07-24 18:36:16 +05:30
read ( 777 , rec = 1 ) C_volAvg
close ( 777 )
2013-09-18 19:37:55 +05:30
call IO_read_realFile ( 777 , 'C_volAvgLastInc' , trim ( getSolverJobName ( ) ) , size ( C_volAvgLastInc ) )
2013-07-24 18:36:16 +05:30
read ( 777 , rec = 1 ) C_volAvgLastInc
close ( 777 )
2014-03-25 21:14:16 +05:30
call IO_read_realFile ( 777 , 'C_ref' , trim ( getSolverJobName ( ) ) , size ( C_minMaxAvg ) )
2013-07-24 18:36:16 +05:30
read ( 777 , rec = 1 ) C_minMaxAvg
close ( 777 )
endif
2014-03-25 21:14:16 +05:30
call Utilities_updateGamma ( C_minMaxAvg , . True . )
C_scale = C_minMaxAvg
S_scale = math_invSym3333 ( C_minMaxAvg )
2013-07-24 18:36:16 +05:30
end subroutine Polarisation_init
!--------------------------------------------------------------------------------------------------
!> @brief solution for the Polarisation scheme with internal iterations
!--------------------------------------------------------------------------------------------------
type ( tSolutionState ) function &
2013-07-26 21:55:37 +05:30
Polarisation_solution ( incInfoIn , guess , timeinc , timeinc_old , loadCaseTime , P_BC , F_BC , temperature_bc , rotation_BC , density )
2013-07-24 18:36:16 +05:30
use numerics , only : &
2013-07-30 21:02:55 +05:30
update_gamma
2013-07-24 18:36:16 +05:30
use math , only : &
2013-12-18 15:05:05 +05:30
math_invSym3333
2013-07-24 18:36:16 +05:30
use DAMASK_spectral_Utilities , only : &
tBoundaryCondition , &
Utilities_maskedCompliance , &
2013-12-18 15:05:05 +05:30
Utilities_updateGamma
2013-07-24 18:36:16 +05:30
use FEsolving , only : &
restartWrite , &
terminallyIll
implicit none
!--------------------------------------------------------------------------------------------------
! input data for solution
2013-07-30 21:02:55 +05:30
real ( pReal ) , intent ( in ) :: &
2013-07-24 18:36:16 +05:30
timeinc , & !< increment in time for current solution
timeinc_old , & !< increment in time of last increment
loadCaseTime , & !< remaining time of current load case
2013-07-30 21:02:55 +05:30
temperature_bc , &
density
2013-07-24 18:36:16 +05:30
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
2013-07-30 21:02:55 +05:30
2013-07-24 18:36:16 +05:30
!--------------------------------------------------------------------------------------------------
! PETSc Data
PetscErrorCode :: ierr
SNESConvergedReason :: reason
2013-08-02 19:25:44 +05:30
incInfo = incInfoIn
2013-07-24 18:36:16 +05:30
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
S = Utilities_maskedCompliance ( rotation_BC , P_BC % maskLogical , C_volAvg )
if ( update_gamma ) then
call Utilities_updateGamma ( C_minMaxAvg , restartWrite )
C_scale = C_minMaxAvg
S_scale = math_invSym3333 ( C_minMaxAvg )
endif
2013-12-20 16:19:14 +05:30
Polarisation_solution % converged = . false .
!--------------------------------------------------------------------------------------------------
! set module wide availabe data
2013-07-24 18:36:16 +05:30
mask_stress = P_BC % maskFloat
params % P_BC = P_BC % values
params % rotation_BC = rotation_BC
params % timeinc = timeinc
2013-07-26 21:55:37 +05:30
params % timeincOld = timeinc_old
2013-07-24 18:36:16 +05:30
params % temperature = temperature_bc
2013-07-26 21:55:37 +05:30
params % density = density
2013-07-24 18:36:16 +05:30
!--------------------------------------------------------------------------------------------------
! solve BVP
2013-07-30 21:02:55 +05:30
call SNESSolve ( snes , PETSC_NULL_OBJECT , solution_vec , ierr )
CHKERRQ ( ierr )
2013-07-24 18:36:16 +05:30
!--------------------------------------------------------------------------------------------------
! check convergence
2013-07-30 21:02:55 +05:30
call SNESGetConvergedReason ( snes , reason , ierr )
CHKERRQ ( ierr )
Polarisation_solution % termIll = terminallyIll
terminallyIll = . false .
Polarisation_solution % converged = . true .
if ( reason < 1 ) Polarisation_solution % converged = . false .
Polarisation_solution % iterationsNeeded = totalIter
2013-07-24 18:36:16 +05:30
end function Polarisation_solution
!--------------------------------------------------------------------------------------------------
!> @brief forms the Polarisation residual vector
!--------------------------------------------------------------------------------------------------
subroutine Polarisation_formResidual ( in , x_scal , f_scal , dummy , ierr )
use numerics , only : &
itmax , &
itmin , &
polarAlpha , &
2013-08-02 19:25:44 +05:30
polarBeta
2013-07-24 18:36:16 +05:30
use IO , only : &
IO_intOut
use math , only : &
math_rotate_backward33 , &
math_transpose33 , &
math_mul3333xx33 , &
2013-07-30 21:02:55 +05:30
math_invSym3333 , &
2013-08-07 22:50:05 +05:30
math_mul33x33 , &
2013-07-30 21:02:55 +05:30
PI
2013-07-24 18:36:16 +05:30
use DAMASK_spectral_Utilities , only : &
grid , &
2013-07-26 21:55:37 +05:30
geomSize , &
2013-07-24 18:36:16 +05:30
wgt , &
field_real , &
2013-07-26 21:55:37 +05:30
field_fourier , &
2013-07-24 18:36:16 +05:30
Utilities_FFTforward , &
Utilities_fourierConvolution , &
2013-07-26 21:55:37 +05:30
Utilities_inverseLaplace , &
2013-07-24 18:36:16 +05:30
Utilities_FFTbackward , &
2013-08-07 22:50:05 +05:30
Utilities_constitutiveResponse , &
Utilities_divergenceRMS , &
Utilities_curlRMS
2013-07-24 18:36:16 +05:30
use debug , only : &
debug_level , &
debug_spectral , &
debug_spectralRotation
use homogenization , only : &
materialpoint_dPdF
implicit none
!--------------------------------------------------------------------------------------------------
! strange syntax in the next line because otherwise macros expand beyond 132 character limit
DMDALocalInfo , dimension ( &
DMDA_LOCAL_INFO_SIZE ) :: &
in
PetscScalar , target , dimension ( 3 , 3 , 2 , &
XG_RANGE , YG_RANGE , ZG_RANGE ) :: &
x_scal
PetscScalar , target , dimension ( 3 , 3 , 2 , &
X_RANGE , Y_RANGE , Z_RANGE ) :: &
f_scal
PetscScalar , pointer , dimension ( : , : , : , : , : ) :: &
F , &
F_tau , &
residual_F , &
residual_F_tau
PetscInt :: &
PETScIter , &
nfuncs
PetscObject :: dummy
PetscErrorCode :: ierr
integer ( pInt ) :: &
i , j , k , e
2013-07-30 21:02:55 +05:30
2013-07-24 18:36:16 +05:30
F = > x_scal ( 1 : 3 , 1 : 3 , 1 , &
XG_RANGE , YG_RANGE , ZG_RANGE )
F_tau = > x_scal ( 1 : 3 , 1 : 3 , 2 , &
XG_RANGE , YG_RANGE , ZG_RANGE )
residual_F = > f_scal ( 1 : 3 , 1 : 3 , 1 , &
X_RANGE , Y_RANGE , Z_RANGE )
residual_F_tau = > f_scal ( 1 : 3 , 1 : 3 , 2 , &
X_RANGE , Y_RANGE , Z_RANGE )
call SNESGetNumberFunctionEvals ( snes , nfuncs , ierr ) ; CHKERRQ ( ierr )
call SNESGetIterationNumber ( snes , PETScIter , ierr ) ; CHKERRQ ( ierr )
2013-08-09 21:55:13 +05:30
F_av = sum ( sum ( sum ( F , dim = 5 ) , dim = 4 ) , dim = 3 ) * wgt
2013-07-30 21:02:55 +05:30
if ( nfuncs == 0 . and . PETScIter == 0 ) totalIter = - 1_pInt ! new increment
if ( totalIter < = PETScIter ) then ! new iteration
2013-07-24 18:36:16 +05:30
!--------------------------------------------------------------------------------------------------
! report begin of new iteration
totalIter = totalIter + 1_pInt
write ( 6 , '(1x,a,3(a,' / / IO_intOut ( itmax ) / / '))' ) trim ( incInfo ) , &
' @ Iteration ' , itmin , '≤' , totalIter , '≤' , itmax
if ( iand ( debug_level ( debug_spectral ) , debug_spectralRotation ) / = 0 ) &
write ( 6 , '(/,a,/,3(3(f12.7,1x)/))' , advance = 'no' ) ' deformation gradient aim (lab) =' , &
math_transpose33 ( math_rotate_backward33 ( F_aim , params % rotation_BC ) )
write ( 6 , '(/,a,/,3(3(f12.7,1x)/))' , advance = 'no' ) ' deformation gradient aim =' , &
math_transpose33 ( F_aim )
flush ( 6 )
endif
2013-07-26 21:55:37 +05:30
!--------------------------------------------------------------------------------------------------
! evaluate inertia
2013-07-30 21:02:55 +05:30
dynamic : if ( params % density > 0.0_pReal ) then
residual_F = ( ( F - F_lastInc ) / params % timeinc - ( F_lastInc - F_lastInc2 ) / params % timeincOld ) / &
( ( params % timeinc + params % timeincOld ) / 2.0_pReal )
2013-07-26 21:55:37 +05:30
residual_F = params % density * product ( geomSize / grid ) * residual_F
field_real = 0.0_pReal
field_real ( 1 : grid ( 1 ) , 1 : grid ( 2 ) , 1 : grid ( 3 ) , 1 : 3 , 1 : 3 ) = reshape ( residual_F , [ grid ( 1 ) , grid ( 2 ) , grid ( 3 ) , 3 , 3 ] , &
2013-07-30 21:02:55 +05:30
order = [ 4 , 5 , 1 , 2 , 3 ] ) ! field real has a different order
2013-07-26 21:55:37 +05:30
call Utilities_FFTforward ( )
call Utilities_inverseLaplace ( )
inertiaField_fourier = field_fourier
2013-07-30 21:02:55 +05:30
else dynamic
2013-07-26 21:55:37 +05:30
inertiaField_fourier = cmplx ( 0.0_pReal , 0.0_pReal , pReal )
2013-07-30 21:02:55 +05:30
endif dynamic
2013-07-24 18:36:16 +05:30
!--------------------------------------------------------------------------------------------------
!
field_real = 0.0_pReal
do k = 1_pInt , grid ( 3 ) ; do j = 1_pInt , grid ( 2 ) ; do i = 1_pInt , grid ( 1 )
2013-08-07 22:50:05 +05:30
field_real ( i , j , k , 1 : 3 , 1 : 3 ) = &
polarBeta * math_mul3333xx33 ( C_scale , F ( 1 : 3 , 1 : 3 , i , j , k ) - math_I3 ) - &
polarAlpha * math_mul33x33 ( F ( 1 : 3 , 1 : 3 , i , j , k ) , &
math_mul3333xx33 ( C_scale , F_tau ( 1 : 3 , 1 : 3 , i , j , k ) - F ( 1 : 3 , 1 : 3 , i , j , k ) - math_I3 ) )
2013-07-24 18:36:16 +05:30
enddo ; enddo ; enddo
!--------------------------------------------------------------------------------------------------
! doing convolution in Fourier space
call Utilities_FFTforward ( )
2013-07-26 21:55:37 +05:30
field_fourier = field_fourier + polarAlpha * inertiaField_fourier
2013-07-24 18:36:16 +05:30
call Utilities_fourierConvolution ( math_rotate_backward33 ( polarBeta * F_aim , params % rotation_BC ) )
call Utilities_FFTbackward ( )
2013-08-07 22:50:05 +05:30
2013-07-24 18:36:16 +05:30
!--------------------------------------------------------------------------------------------------
! constructing residual
residual_F_tau = polarBeta * F - reshape ( field_real ( 1 : grid ( 1 ) , 1 : grid ( 2 ) , 1 : grid ( 3 ) , 1 : 3 , 1 : 3 ) , &
[ 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ] , order = [ 3 , 4 , 5 , 1 , 2 ] )
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
2013-07-30 21:02:55 +05:30
P_avLastEval = P_av
2013-07-24 18:36:16 +05:30
call Utilities_constitutiveResponse ( F_lastInc , F - residual_F_tau / polarBeta , params % temperature , params % timeinc , &
residual_F , C_volAvg , C_minMaxAvg , P_av , ForwardData , params % rotation_BC )
ForwardData = . False .
2013-08-07 22:50:05 +05:30
2013-07-30 21:02:55 +05:30
!--------------------------------------------------------------------------------------------------
2013-08-07 22:50:05 +05:30
! calculate divergence
field_real = 0.0_pReal
field_real = reshape ( residual_F , [ grid ( 1 ) , grid ( 2 ) , grid ( 3 ) , 3 , 3 ] , order = [ 4 , 5 , 1 , 2 , 3 ] )
call Utilities_FFTforward ( )
err_div = Utilities_divergenceRMS ( )
call Utilities_FFTbackward ( )
2013-07-30 21:02:55 +05:30
2013-07-24 18:36:16 +05:30
!--------------------------------------------------------------------------------------------------
! constructing residual
e = 0_pInt
do k = 1_pInt , grid ( 3 ) ; do j = 1_pInt , grid ( 2 ) ; do i = 1_pInt , grid ( 1 )
e = e + 1_pInt
2013-08-08 23:40:58 +05:30
residual_F ( 1 : 3 , 1 : 3 , i , j , k ) = &
2014-03-25 21:14:16 +05:30
math_mul3333xx33 ( math_invSym3333 ( materialpoint_dPdF ( 1 : 3 , 1 : 3 , 1 : 3 , 1 : 3 , 1 , e ) + C_scale ) , &
2013-08-09 09:46:16 +05:30
residual_F ( 1 : 3 , 1 : 3 , i , j , k ) - math_mul33x33 ( F ( 1 : 3 , 1 : 3 , i , j , k ) , &
math_mul3333xx33 ( C_scale , F_tau ( 1 : 3 , 1 : 3 , i , j , k ) - F ( 1 : 3 , 1 : 3 , i , j , k ) - math_I3 ) ) ) &
+ residual_F_tau ( 1 : 3 , 1 : 3 , i , j , k )
2013-07-24 18:36:16 +05:30
enddo ; enddo ; enddo
!--------------------------------------------------------------------------------------------------
2013-08-07 22:50:05 +05:30
! calculating curl
field_real = 0.0_pReal
field_real = reshape ( F , [ grid ( 1 ) , grid ( 2 ) , grid ( 3 ) , 3 , 3 ] , order = [ 4 , 5 , 1 , 2 , 3 ] )
call Utilities_FFTforward ( )
err_curl = Utilities_curlRMS ( )
call Utilities_FFTbackward ( )
2013-07-30 21:02:55 +05:30
2013-07-24 18:36:16 +05:30
end subroutine Polarisation_formResidual
!--------------------------------------------------------------------------------------------------
!> @brief convergence check
!--------------------------------------------------------------------------------------------------
subroutine Polarisation_converged ( snes_local , PETScIter , xnorm , snorm , fnorm , reason , dummy , ierr )
use numerics , only : &
2013-08-08 14:43:29 +05:30
itmax , &
itmin , &
err_div_tolRel , &
err_div_tolAbs , &
err_curl_tolRel , &
err_curl_tolAbs , &
err_stress_tolabs , &
err_stress_tolrel
use math , only : &
math_mul3333xx33
2013-07-24 18:36:16 +05:30
use FEsolving , only : &
terminallyIll
implicit none
SNES :: snes_local
PetscInt :: PETScIter
PetscReal :: &
xnorm , &
snorm , &
fnorm
SNESConvergedReason :: reason
PetscObject :: dummy
PetscErrorCode :: ierr
2013-07-30 21:02:55 +05:30
real ( pReal ) :: &
2013-08-07 22:50:05 +05:30
curlTol , &
divTol , &
2013-08-09 21:55:13 +05:30
BC_tol
2013-08-07 22:50:05 +05:30
!--------------------------------------------------------------------------------------------------
! stress BC handling
F_aim = F_aim - math_mul3333xx33 ( S , ( ( P_av - params % P_BC ) ) ) ! S = 0.0 for no bc
2013-09-18 19:37:55 +05:30
err_BC = maxval ( abs ( ( - mask_stress + 1.0_pReal ) * math_mul3333xx33 ( C_scale , F_aim - F_av ) + &
mask_stress * ( P_av - params % P_BC ) ) ) ! mask = 0.0 for no bc
2013-08-08 14:43:29 +05:30
!--------------------------------------------------------------------------------------------------
! error calculation
curlTol = max ( maxval ( abs ( F_aim - math_I3 ) ) * err_curl_tolRel , err_curl_tolAbs )
divTol = max ( maxval ( abs ( P_av ) ) * err_div_tolRel , err_div_tolAbs )
2013-08-09 21:55:13 +05:30
BC_tol = max ( maxval ( abs ( P_av ) ) * err_stress_tolrel , err_stress_tolabs )
2013-08-07 22:50:05 +05:30
2013-07-30 21:02:55 +05:30
converged : if ( ( totalIter > = itmin . and . &
2013-08-07 22:50:05 +05:30
all ( [ err_div / divTol , &
err_curl / curlTol , &
2013-08-09 21:55:13 +05:30
err_BC / BC_tol ] < 1.0_pReal ) ) &
2013-08-08 14:43:29 +05:30
. or . terminallyIll ) then
2013-07-30 21:02:55 +05:30
reason = 1
elseif ( totalIter > = itmax ) then converged
reason = - 1
else converged
reason = 0
endif converged
2013-08-07 22:50:05 +05:30
!--------------------------------------------------------------------------------------------------
! report
write ( 6 , '(1/,a)' ) ' ... reporting .............................................................'
2013-08-08 14:43:29 +05:30
write ( 6 , '(/,a,f12.2,a,es8.2,a,es9.2,a)' ) ' error curl = ' , &
err_curl / curlTol , ' (' , err_curl , ' -, tol =' , curlTol , ')'
2013-09-18 19:37:55 +05:30
write ( 6 , ' (a,f12.2,a,es8.2,a,es9.2,a)' ) ' error divergence = ' , &
err_div / divTol , ' (' , err_div , ' / m, tol =' , divTol , ')'
write ( 6 , ' (a,f12.2,a,es8.2,a,es9.2,a)' ) ' error BC = ' , &
err_BC / BC_tol , ' (' , err_BC , ' Pa, tol =' , BC_tol , ')'
2013-08-07 22:50:05 +05:30
write ( 6 , '(/,a)' ) ' ==========================================================================='
flush ( 6 )
2013-07-24 18:36:16 +05:30
end subroutine Polarisation_converged
2013-12-18 15:05:05 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief forwarding routine
!--------------------------------------------------------------------------------------------------
subroutine Polarisation_forward ( guess , timeinc , timeinc_old , loadCaseTime , F_BC , P_BC , rotation_BC )
use math , only : &
math_mul33x33 , &
math_mul3333xx33 , &
math_transpose33 , &
math_rotate_backward33
use DAMASK_spectral_Utilities , only : &
grid , &
geomSize , &
Utilities_calculateRate , &
Utilities_forwardField , &
tBoundaryCondition , &
cutBack
use mesh , only : &
mesh_ipCoordinates , &
mesh_deformedCoordsFFT
2013-12-20 16:19:14 +05:30
use IO , only : &
IO_write_JobRealFile
use FEsolving , only : &
restartWrite
2013-12-18 15:05:05 +05:30
implicit none
#include <finclude/petscdmda.h90>
real ( pReal ) , intent ( in ) :: &
timeinc_old , &
timeinc , &
loadCaseTime !< remaining time of current load case
type ( tBoundaryCondition ) , intent ( in ) :: &
P_BC , &
F_BC
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: rotation_BC
logical , intent ( in ) :: &
guess
PetscErrorCode :: ierr
PetscScalar , dimension ( : , : , : , : ) , pointer :: xx_psc , F , F_tau
integer ( pInt ) :: i , j , k
real ( pReal ) , dimension ( 3 , 3 ) :: F_lambda33
!--------------------------------------------------------------------------------------------------
! update coordinates and rate and forward last inc
call DMDAVecGetArrayF90 ( da , solution_vec , xx_psc , ierr )
F = > xx_psc ( 0 : 8 , : , : , : )
F_tau = > xx_psc ( 9 : 17 , : , : , : )
2013-12-20 16:19:14 +05:30
if ( restartWrite ) then
write ( 6 , '(/,a)' ) ' writing converged results for restart'
flush ( 6 )
call IO_write_jobRealFile ( 777 , 'F' , size ( F ) ) ! writing deformation gradient field to file
write ( 777 , rec = 1 ) F
close ( 777 )
call IO_write_jobRealFile ( 777 , 'F_lastInc' , size ( F_lastInc ) ) ! writing F_lastInc field to file
write ( 777 , rec = 1 ) F_lastInc
close ( 777 )
call IO_write_jobRealFile ( 777 , 'F_lastInc2' , size ( F_lastInc2 ) ) ! writing F_lastInc field to file
write ( 777 , rec = 1 ) F_lastInc2
close ( 777 )
call IO_write_jobRealFile ( 777 , 'F_tau' , size ( F_tau ) ) ! writing deformation gradient field to file
write ( 777 , rec = 1 ) F_tau
close ( 777 )
call IO_write_jobRealFile ( 777 , 'F_tau_lastInc' , size ( F_tau_lastInc ) ) ! writing F_lastInc field to file
write ( 777 , rec = 1 ) F_tau_lastInc
close ( 777 )
2014-03-24 15:29:30 +05:30
call IO_write_jobRealFile ( 777 , 'F_aim' , size ( F_aim ) )
write ( 777 , rec = 1 ) F_aim
close ( 777 )
call IO_write_jobRealFile ( 777 , 'F_aim_lastInc' , size ( F_aim_lastInc ) )
write ( 777 , rec = 1 ) F_aim_lastInc
close ( 777 )
2013-12-20 16:19:14 +05:30
call IO_write_jobRealFile ( 777 , 'F_aimDot' , size ( F_aimDot ) )
write ( 777 , rec = 1 ) F_aimDot
close ( 777 )
call IO_write_jobRealFile ( 777 , 'C_volAvg' , size ( C_volAvg ) )
write ( 777 , rec = 1 ) C_volAvg
close ( 777 )
call IO_write_jobRealFile ( 777 , 'C_volAvgLastInc' , size ( C_volAvgLastInc ) )
write ( 777 , rec = 1 ) C_volAvgLastInc
close ( 777 )
endif
mesh_ipCoordinates = reshape ( mesh_deformedCoordsFFT ( geomSize , reshape ( &
F , [ 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ] ) ) , [ 3 , 1 , product ( grid ) ] )
2013-12-18 15:05:05 +05:30
if ( cutBack ) then
F_aim = F_aim_lastInc
F_tau = reshape ( F_tau_lastInc , [ 9 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ] )
F = reshape ( F_lastInc , [ 9 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ] )
C_volAvg = C_volAvgLastInc
else
2014-01-27 16:35:36 +05:30
ForwardData = . True .
2013-12-18 15:05:05 +05:30
C_volAvgLastInc = C_volAvg
!--------------------------------------------------------------------------------------------------
! calculate rate for aim
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
f_aimDot = F_BC % maskFloat * F_BC % values
elseif ( F_BC % myType == 'f' ) then ! aim at end of load case is prescribed
f_aimDot = F_BC % maskFloat * ( F_BC % values - F_aim ) / loadCaseTime
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
mesh_ipCoordinates = reshape ( mesh_deformedCoordsFFT ( geomSize , reshape ( &
F , [ 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ] ) ) , [ 3 , 1 , product ( grid ) ] )
Fdot = Utilities_calculateRate ( math_rotate_backward33 ( f_aimDot , rotation_BC ) , &
timeinc_old , guess , F_lastInc , reshape ( F , [ 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ] ) )
F_tauDot = Utilities_calculateRate ( math_rotate_backward33 ( 2.0_pReal * f_aimDot , rotation_BC ) , &
timeinc_old , guess , F_tau_lastInc , reshape ( F_tau , [ 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ] ) )
F_lastInc2 = F_lastInc
F_lastInc = reshape ( F , [ 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ] )
F_tau_lastInc = reshape ( F_tau , [ 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ] )
endif
F_aim = F_aim + f_aimDot * timeinc
!--------------------------------------------------------------------------------------------------
! update local deformation gradient
F = reshape ( Utilities_forwardField ( timeinc , F_lastInc , Fdot , & ! ensure that it matches rotated F_aim
math_rotate_backward33 ( F_aim , rotation_BC ) ) , [ 9 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ] )
F_tau = reshape ( Utilities_forwardField ( timeinc , F_tau_lastInc , F_taudot ) , [ 9 , grid ( 1 ) , grid ( 2 ) , grid ( 3 ) ] ) ! does not have any average value as boundary condition
if ( . not . guess ) then ! large strain forwarding
do k = 1_pInt , grid ( 3 ) ; do j = 1_pInt , grid ( 2 ) ; do i = 1_pInt , grid ( 1 )
2014-03-26 14:39:28 +05:30
F_lambda33 = reshape ( F_tau ( 1 : 9 , i , j , k ) - F ( 1 : 9 , i , j , k ) , [ 3 , 3 ] )
2013-12-18 15:05:05 +05:30
F_lambda33 = math_mul3333xx33 ( S_scale , math_mul33x33 ( F_lambda33 , &
math_mul3333xx33 ( C_scale , &
math_mul33x33 ( math_transpose33 ( F_lambda33 ) , &
F_lambda33 ) - math_I3 ) ) * 0.5_pReal ) &
+ math_I3
2014-03-26 14:39:28 +05:30
F_tau ( 1 : 9 , i , j , k ) = reshape ( F_lambda33 , [ 9 ] ) + F ( 1 : 9 , i , j , k )
2013-12-18 15:05:05 +05:30
enddo ; enddo ; enddo
endif
call DMDAVecRestoreArrayF90 ( da , solution_vec , xx_psc , ierr ) ; CHKERRQ ( ierr )
end subroutine Polarisation_forward
2013-07-24 18:36:16 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief destroy routine
!--------------------------------------------------------------------------------------------------
subroutine Polarisation_destroy ( )
use DAMASK_spectral_Utilities , only : &
Utilities_destroy
implicit none
PetscErrorCode :: ierr
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 Polarisation_destroy
end module DAMASK_spectral_SolverPolarisation