223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

added pheno+ module 

modify crystallite microstructure call
to pass orientations
This commit is contained in:
Chen Zhang 2015-10-14 18:36:19 +00:00
parent 121c471455
commit 484a34b7f1
4 changed files with 217 additions and 200 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

227
code/DAMASK_spectral_utilities.f90
View File

@ -13,7 +13,7 @@ module DAMASK_spectral_utilities
pInt
use math, only: &
math_I3
use numerics, only: &
use numerics, only: &
spectral_filter
implicit none
@ -22,11 +22,11 @@ module DAMASK_spectral_utilities
#include <petsc/finclude/petscsys.h>
#endif
include 'fftw3-mpi.f03'
logical, public :: cutBack =.false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill
integer(pInt), public, parameter :: maxPhaseFields = 2_pInt
integer(pInt), public :: nActiveFields = 0_pInt
!--------------------------------------------------------------------------------------------------
! field labels information
enum, bind(c)
@ -36,11 +36,11 @@ module DAMASK_spectral_utilities
FIELD_DAMAGE_ID, &
FIELD_VACANCYDIFFUSION_ID
end enum
!--------------------------------------------------------------------------------------------------
! grid related information information
real(pReal), public :: wgt !< weighting factor 1/Nelems
!--------------------------------------------------------------------------------------------------
! variables storing information for spectral method and FFTW
integer(pInt), public :: grid1Red !< grid(1)/2
@ -55,7 +55,7 @@ module DAMASK_spectral_utilities
real(pReal), private, dimension(:,:,:,:), allocatable :: xi2nd !< wave vector field for second derivatives
real(pReal), private, dimension(3,3,3,3) :: C_ref !< mechanic reference stiffness
real(pReal), protected, public, dimension(3) :: scaledGeomSize !< scaled geometry size for calculation of divergence (Basic, Basic PETSc)
!--------------------------------------------------------------------------------------------------
! plans for FFTW
type(C_PTR), private :: &
@ -76,10 +76,10 @@ module DAMASK_spectral_utilities
!--------------------------------------------------------------------------------------------------
! derived types
type, public :: tSolutionState !< return type of solution from spectral solver variants
logical :: converged = .true.
logical :: regrid = .false.
logical :: stagConverged = .true.
logical :: termIll = .false.
logical :: converged = .true.
logical :: regrid = .false.
logical :: stagConverged = .true.
logical :: termIll = .false.
integer(pInt) :: iterationsNeeded = 0_pInt
end type tSolutionState
@ -99,35 +99,35 @@ module DAMASK_spectral_utilities
outputfrequency = 1_pInt, & !< frequency of result writes
restartfrequency = 0_pInt, & !< frequency of restart writes
logscale = 0_pInt !< linear/logarithmic time inc flag
logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase
logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase
integer(kind(FIELD_UNDEFINED_ID)), allocatable :: ID(:)
end type tLoadCase
type, public :: 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
real(pReal) :: timeincOld
real(pReal) :: density
end type tSolutionParams
type(tSolutionParams), private :: params
type, public :: phaseFieldDataBin !< set of parameters defining a phase field
real(pReal) :: diffusion = 0.0_pReal, & !< thermal conductivity
mobility = 0.0_pReal, & !< thermal mobility
phaseField0 = 0.0_pReal !< homogeneous damage field starting condition
phaseField0 = 0.0_pReal !< homogeneous damage field starting condition
logical :: active = .false.
character(len=64) :: label = ''
end type phaseFieldDataBin
enum, bind(c)
enumerator :: FILTER_NONE_ID, &
FILTER_GRADIENT_ID, &
enum, bind(c)
enumerator :: FILTER_NONE_ID, &
FILTER_GRADIENT_ID, &
FILTER_COSINE_ID
end enum
integer(kind(FILTER_NONE_ID)) :: &
spectral_filter_ID
public :: &
utilities_init, &
utilities_updateGamma, &
@ -156,11 +156,11 @@ module DAMASK_spectral_utilities
private :: &
utilities_getFilter
contains
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
!> @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.
@ -193,7 +193,7 @@ subroutine utilities_init()
use debug, only: &
PETSCDEBUG
#endif
use math
use math
use mesh, only: &
grid, &
grid3, &
@ -207,7 +207,7 @@ subroutine utilities_init()
PETScOptionsInsertString, &
MPI_Abort
PetscErrorCode :: ierr
#endif
#endif
integer(pInt) :: i, j, k
integer(pInt), dimension(3) :: k_s
type(C_PTR) :: &
@ -220,7 +220,7 @@ subroutine utilities_init()
scalarSize = 1_C_INTPTR_T, &
vecSize = 3_C_INTPTR_T, &
tensorSize = 9_C_INTPTR_T
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- DAMASK_spectral_utilities init -+>>>'
write(6,'(a)') ' $Id$'
@ -251,7 +251,7 @@ subroutine utilities_init()
write(6,'(a,3(i12 ))') ' grid a b c: ', grid
write(6,'(a,3(es12.5))') ' size x y z: ', geomSize
endif
!--------------------------------------------------------------------------------------------------
! scale dimension to calculate either uncorrected, dimension-independent, or dimension- and
! resolution-independent divergence
@ -277,7 +277,7 @@ subroutine utilities_init()
MPI_COMM_WORLD, local_K, local_K_offset)
allocate (xi1st(3,grid1Red,grid(2),grid3),source = 0.0_pReal) ! frequencies, only half the size for first dimension
allocate (xi2nd(3,grid1Red,grid(2),grid3),source = 0.0_pReal) ! frequencies, only half the size for first dimension
tensorField = fftw_alloc_complex(tensorSize*alloc_local)
call c_f_pointer(tensorField, tensorField_real, [3_C_INTPTR_T,3_C_INTPTR_T, &
2_C_INTPTR_T*(gridFFTW(1)/2_C_INTPTR_T + 1_C_INTPTR_T),gridFFTW(2),local_K]) ! place a pointer for a real tensor representation
@ -295,7 +295,7 @@ subroutine utilities_init()
[2_C_INTPTR_T*(gridFFTW(1)/2_C_INTPTR_T + 1),gridFFTW(2),local_K]) ! place a pointer for a real scalar representation
call c_f_pointer(scalarField, scalarField_fourier, &
[ gridFFTW(1)/2_C_INTPTR_T + 1 ,gridFFTW(2),local_K]) ! place a pointer for a fourier scarlar representation
!--------------------------------------------------------------------------------------------------
! tensor MPI fftw plans
planTensorForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order
@ -308,7 +308,7 @@ subroutine utilities_init()
tensorField_fourier,tensorField_real, & ! input data, output data
MPI_COMM_WORLD, fftw_planner_flag) ! all processors, planer precision
if (.not. C_ASSOCIATED(planTensorBack)) call IO_error(810, ext_msg='planTensorBack')
!--------------------------------------------------------------------------------------------------
! vector MPI fftw plans
planVectorForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order
@ -321,10 +321,10 @@ subroutine utilities_init()
vectorField_fourier,vectorField_real, & ! input data, output data
MPI_COMM_WORLD, fftw_planner_flag) ! all processors, planer precision
if (.not. C_ASSOCIATED(planVectorBack)) call IO_error(810, ext_msg='planVectorBack')
!--------------------------------------------------------------------------------------------------
! scalar MPI fftw plans
planScalarForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order
! scalar MPI fftw plans
planScalarForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order
scalarSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, & ! no. of transforms, default iblock and oblock
scalarField_real, scalarField_fourier, & ! input data, output data
MPI_COMM_WORLD, fftw_planner_flag) ! use all processors, planer precision
@ -342,7 +342,7 @@ subroutine utilities_init()
if (debugGeneral .and. worldrank == 0_pInt) write(6,'(/,a)') ' FFTW initialized'
flush(6)
!--------------------------------------------------------------------------------------------------
! calculation of discrete angular frequencies, ordered as in FFTW (wrap around)
do k = grid3Offset+1_pInt, grid3Offset+grid3
@ -360,7 +360,7 @@ subroutine utilities_init()
xi1st(1:3,i,j,k-grid3Offset) = xi2nd(1:3,i,j,k-grid3Offset)
endwhere
enddo; enddo; enddo
if(memory_efficient) then ! allocate just single fourth order tensor
allocate (gamma_hat(3,3,3,3,1,1,1), source = 0.0_pReal)
else ! precalculation of gamma_hat field
@ -408,7 +408,7 @@ subroutine utilities_updateGamma(C,saveReference)
integer(pInt) :: &
i, j, k, &
l, m, n, o
C_ref = C
if (saveReference) then
if (worldrank == 0_pInt) then
@ -419,8 +419,8 @@ subroutine utilities_updateGamma(C,saveReference)
close(777)
endif
endif
if(.not. memory_efficient) then
if(.not. memory_efficient) then
do k = grid3Offset+1_pInt, grid3Offset+grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid1Red
if (any([i,j,k] /= 1_pInt)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) &
@ -430,9 +430,9 @@ subroutine utilities_updateGamma(C,saveReference)
temp33_Real = math_inv33(temp33_Real)
forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, o=1_pInt:3_pInt)&
gamma_hat(l,m,n,o,i,j,k-grid3Offset) = temp33_Real(l,n)*xiDyad(m,o)
endif
endif
enddo; enddo; enddo
endif
endif
end subroutine utilities_updateGamma
@ -451,13 +451,13 @@ subroutine utilities_FFTtensorForward()
!--------------------------------------------------------------------------------------------------
! doing the tensor FFT
call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier)
!--------------------------------------------------------------------------------------------------
! applying filter
forall(k = 1_pInt:grid3, j = 1_pInt:grid(2), i = 1_pInt:grid1Red) &
tensorField_fourier(1:3,1:3,i,j,k) = utilities_getFilter(xi2nd(1:3,i,j,k))* &
tensorField_fourier(1:3,1:3,i,j,k)
end subroutine utilities_FFTtensorForward
@ -481,10 +481,10 @@ subroutine utilities_FFTscalarForward()
use mesh, only: &
grid3, &
grid
implicit none
integer(pInt) :: i, j, k
!--------------------------------------------------------------------------------------------------
! doing the scalar FFT
call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_real,scalarField_fourier)
@ -494,7 +494,7 @@ subroutine utilities_FFTscalarForward()
forall(k = 1_pInt:grid3, j = 1_pInt:grid(2), i = 1_pInt:grid1Red) &
scalarField_fourier(i,j,k) = utilities_getFilter(xi2nd(1:3,i,j,k))* &
scalarField_fourier(i,j,k)
end subroutine utilities_FFTscalarForward
@ -519,7 +519,7 @@ subroutine utilities_FFTvectorForward()
use mesh, only: &
grid3, &
grid
implicit none
integer(pInt) :: i, j, k
@ -532,7 +532,7 @@ subroutine utilities_FFTvectorForward()
forall(k = 1_pInt:grid3, j = 1_pInt:grid(2), i = 1_pInt:grid1Red) &
vectorField_fourier(1:3,i,j,k) = utilities_getFilter(xi2nd(1:3,i,j,k))* &
vectorField_fourier(1:3,i,j,k)
end subroutine utilities_FFTvectorForward
@ -558,31 +558,31 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
use math, only: &
math_inv33
use numerics, only: &
worldrank
worldrank
use mesh, only: &
grid3, &
grid, &
grid3Offset
implicit none
implicit none
real(pReal), intent(in), dimension(3,3) :: fieldAim !< desired average value of the field after convolution
real(pReal), dimension(3,3) :: xiDyad, temp33_Real
complex(pReal), dimension(3,3) :: temp33_complex
integer(pInt) :: &
i, j, k, &
l, m, n, o
if (worldrank == 0_pInt) then
if (worldrank == 0_pInt) then
write(6,'(/,a)') ' ... doing gamma convolution ...............................................'
flush(6)
endif
!--------------------------------------------------------------------------------------------------
! do the actual spectral method calculation (mechanical equilibrium)
memoryEfficient: if(memory_efficient) then
do k = 1_pInt, grid3; do j = 1_pInt, grid(2) ;do i = 1_pInt, grid1Red
if(any([i,j,k+grid3Offset] /= 1_pInt)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
if(any([i,j,k+grid3Offset] /= 1_pInt)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) &
xiDyad(l,m) = xi1st(l, i,j,k)*xi1st(m, i,j,k)
forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) &
@ -593,8 +593,8 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) &
temp33_Complex(l,m) = sum(gamma_hat(l,m,1:3,1:3, 1,1,1) * &
tensorField_fourier(1:3,1:3,i,j,k))
tensorField_fourier(1:3,1:3,i,j,k) = temp33_Complex
endif
tensorField_fourier(1:3,1:3,i,j,k) = temp33_Complex
endif
enddo; enddo; enddo
else memoryEfficient
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid1Red
@ -604,12 +604,12 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
tensorField_fourier(1:3,1:3,i,j,k) = temp33_Complex
enddo; enddo; enddo
endif memoryEfficient
if (grid3Offset == 0_pInt) &
tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim/wgt,0.0_pReal,pReal) ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim/wgt,0.0_pReal,pReal) ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
end subroutine utilities_fourierGammaConvolution
!--------------------------------------------------------------------------------------------------
!> @brief doing convolution DamageGreenOp_hat * field_real
@ -624,13 +624,13 @@ subroutine utilities_fourierGreenConvolution(D_ref, mobility_ref, deltaT)
grid3, &
geomSize
implicit none
implicit none
real(pReal), dimension(3,3), intent(in) :: D_ref !< desired average value of the field after convolution
real(pReal), intent(in) :: mobility_ref, deltaT !< desired average value of the field after convolution
real(pReal), dimension(3) :: k_s
real(pReal) :: GreenOp_hat
integer(pInt) :: i, j, k
!--------------------------------------------------------------------------------------------------
! do the actual spectral method calculation
do k = 1_pInt, grid3; do j = 1_pInt, grid(2) ;do i = 1_pInt, grid1Red
@ -638,7 +638,7 @@ subroutine utilities_fourierGreenConvolution(D_ref, mobility_ref, deltaT)
GreenOp_hat = 1.0_pReal/ &
(mobility_ref + deltaT*sum((2.0_pReal*PI*k_s/geomSize)* &
math_mul33x3(D_ref,(2.0_pReal*PI*k_s/geomSize)))) !< GreenOp_hat = iK^{T} * D_ref * iK, K is frequency
scalarField_fourier(i,j,k) = scalarField_fourier(i,j,k)*GreenOp_hat
scalarField_fourier(i,j,k) = scalarField_fourier(i,j,k)*GreenOp_hat
enddo; enddo; enddo
end subroutine utilities_fourierGreenConvolution
@ -647,22 +647,22 @@ end subroutine utilities_fourierGreenConvolution
!--------------------------------------------------------------------------------------------------
!> @brief calculate root mean square of divergence of field_fourier
!--------------------------------------------------------------------------------------------------
real(pReal) function utilities_divergenceRMS()
real(pReal) function utilities_divergenceRMS()
use math, only: &
TWOPIIMG, &
math_mul33x3_complex
use numerics, only: &
worldrank
worldrank
use mesh, only: &
grid, &
grid3
implicit none
integer(pInt) :: i, j, k
integer(pInt) :: i, j, k
PetscErrorCode :: ierr
if (worldrank == 0_pInt) then
if (worldrank == 0_pInt) then
write(6,'(/,a)') ' ... calculating divergence ................................................'
flush(6)
endif
@ -674,8 +674,8 @@ real(pReal) function utilities_divergenceRMS()
do i = 2_pInt, grid1Red -1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice.
utilities_divergenceRMS = utilities_divergenceRMS &
+ 2.0_pReal*(sum (real(math_mul33x3_complex(tensorField_fourier(1:3,1:3,i,j,k),& ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2. do not take square root and square again
xi1st(1:3,i,j,k))*TWOPIIMG)**2.0_pReal)& ! --> sum squared L_2 norm of vector
+sum(aimag(math_mul33x3_complex(tensorField_fourier(1:3,1:3,i,j,k),&
xi1st(1:3,i,j,k))*TWOPIIMG)**2.0_pReal)& ! --> sum squared L_2 norm of vector
+sum(aimag(math_mul33x3_complex(tensorField_fourier(1:3,1:3,i,j,k),&
xi1st(1:3,i,j,k))*TWOPIIMG)**2.0_pReal))
enddo
utilities_divergenceRMS = utilities_divergenceRMS & ! these two layers (DC and Nyquist) do not have a conjugate complex counterpart (if grid(1) /= 1)
@ -694,25 +694,25 @@ real(pReal) function utilities_divergenceRMS()
end function utilities_divergenceRMS
!--------------------------------------------------------------------------------------------------
!> @brief calculate max of curl of field_fourier
!--------------------------------------------------------------------------------------------------
real(pReal) function utilities_curlRMS()
use math
use math
use numerics, only: &
worldrank
worldrank
use mesh, only: &
grid, &
grid3
implicit none
integer(pInt) :: i, j, k, l
integer(pInt) :: i, j, k, l
complex(pReal), dimension(3,3) :: curl_fourier
PetscErrorCode :: ierr
if (worldrank == 0_pInt) then
if (worldrank == 0_pInt) then
write(6,'(/,a)') ' ... calculating curl ......................................................'
flush(6)
endif
@ -720,9 +720,9 @@ real(pReal) function utilities_curlRMS()
!--------------------------------------------------------------------------------------------------
! calculating max curl criterion in Fourier space
utilities_curlRMS = 0.0_pReal
do k = 1_pInt, grid3; do j = 1_pInt, grid(2);
do i = 2_pInt, grid1Red - 1_pInt
do k = 1_pInt, grid3; do j = 1_pInt, grid(2);
do i = 2_pInt, grid1Red - 1_pInt
do l = 1_pInt, 3_pInt
curl_fourier(l,1) = (+tensorField_fourier(l,3,i,j,k)*xi1st(2,i,j,k)&
-tensorField_fourier(l,2,i,j,k)*xi1st(3,i,j,k))*TWOPIIMG
@ -733,8 +733,8 @@ real(pReal) function utilities_curlRMS()
enddo
utilities_curlRMS = utilities_curlRMS + &
2.0_pReal*sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal)! Has somewhere a conj. complex counterpart. Therefore count it twice.
enddo
do l = 1_pInt, 3_pInt
enddo
do l = 1_pInt, 3_pInt
curl_fourier = (+tensorField_fourier(l,3,1,j,k)*xi1st(2,1,j,k)&
-tensorField_fourier(l,2,1,j,k)*xi1st(3,1,j,k))*TWOPIIMG
curl_fourier = (+tensorField_fourier(l,1,1,j,k)*xi1st(3,1,j,k)&
@ -744,7 +744,7 @@ real(pReal) function utilities_curlRMS()
enddo
utilities_curlRMS = utilities_curlRMS + &
sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal)! this layer (DC) does not have a conjugate complex counterpart (if grid(1) /= 1)
do l = 1_pInt, 3_pInt
do l = 1_pInt, 3_pInt
curl_fourier = (+tensorField_fourier(l,3,grid1Red,j,k)*xi1st(2,grid1Red,j,k)&
-tensorField_fourier(l,2,grid1Red,j,k)*xi1st(3,grid1Red,j,k))*TWOPIIMG
curl_fourier = (+tensorField_fourier(l,1,grid1Red,j,k)*xi1st(3,grid1Red,j,k)&
@ -753,7 +753,7 @@ real(pReal) function utilities_curlRMS()
-tensorField_fourier(l,1,grid1Red,j,k)*xi1st(2,grid1Red,j,k))*TWOPIIMG
enddo
utilities_curlRMS = utilities_curlRMS + &
sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal)! this layer (Nyquist) does not have a conjugate complex counterpart (if grid(1) /= 1)
sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal)! this layer (Nyquist) does not have a conjugate complex counterpart (if grid(1) /= 1)
enddo; enddo
call MPI_Allreduce(MPI_IN_PLACE,utilities_curlRMS,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
@ -783,17 +783,17 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
real(pReal), intent(in) , dimension(3,3,3,3) :: C !< current average stiffness
real(pReal), intent(in) , dimension(3,3) :: rot_BC !< rotation of load frame
logical, intent(in), dimension(3,3) :: mask_stress !< mask of stress BC
integer(pInt) :: j, k, m, n
integer(pInt) :: j, k, m, n
logical, dimension(9) :: mask_stressVector
real(pReal), dimension(9,9) :: temp99_Real
integer(pInt) :: size_reduced = 0_pInt
real(pReal), dimension(9,9) :: temp99_Real
integer(pInt) :: size_reduced = 0_pInt
real(pReal), dimension(:,:), allocatable :: &
s_reduced, & !< reduced compliance matrix (depending on number of stress BC)
c_reduced, & !< reduced stiffness (depending on number of stress BC)
c_reduced, & !< reduced stiffness (depending on number of stress BC)
sTimesC !< temp variable to check inversion
logical :: errmatinv
character(len=1024):: formatString
mask_stressVector = reshape(transpose(mask_stress), [9])
size_reduced = int(count(mask_stressVector), pInt)
if(size_reduced > 0_pInt )then
@ -803,7 +803,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
temp99_Real = math_Plain3333to99(math_rotate_forward3333(C,rot_BC))
if(debugGeneral .and. worldrank == 0_pInt) then
if(debugGeneral .and. worldrank == 0_pInt) then
write(6,'(/,a)') ' ... updating masked compliance ............................................'
write(6,'(/,a,/,9(9(2x,f12.7,1x)/))',advance='no') ' Stiffness C (load) / GPa =',&
transpose(temp99_Real)/1.e9_pReal
@ -819,6 +819,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
j = j + 1_pInt
c_reduced(k,j) = temp99_Real(n,m)
endif; enddo; endif; enddo
call math_invert(size_reduced, c_reduced, s_reduced, errmatinv) ! invert reduced stiffness
if(errmatinv) call IO_error(error_ID=400_pInt,ext_msg='utilities_maskedCompliance')
temp99_Real = 0.0_pReal ! fill up compliance with zeros
@ -832,7 +833,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
j = j + 1_pInt
temp99_Real(n,m) = s_reduced(k,j)
endif; enddo; endif; enddo
!--------------------------------------------------------------------------------------------------
! check if inversion was successful
sTimesC = matmul(c_reduced,s_reduced)
@ -862,7 +863,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
flush(6)
utilities_maskedCompliance = math_Plain99to3333(temp99_Real)
end function utilities_maskedCompliance
end function utilities_maskedCompliance
!--------------------------------------------------------------------------------------------------
@ -917,7 +918,7 @@ end subroutine utilities_fourierVectorDivergence
!--------------------------------------------------------------------------------------------------
!> @brief calculates constitutive response
!--------------------------------------------------------------------------------------------------
subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc,&
subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
P,C_volAvg,C_minmaxAvg,P_av,forwardData,rotation_BC)
use debug, only: &
debug_reset, &
@ -943,7 +944,7 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc,&
materialpoint_F, &
materialpoint_P, &
materialpoint_dPdF
implicit none
real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: &
F_lastInc, & !< target deformation gradient
@ -951,33 +952,33 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc,&
real(pReal), intent(in) :: timeinc !< loading time
logical, intent(in) :: forwardData !< age results
real(pReal), intent(in), dimension(3,3) :: rotation_BC !< rotation of load frame
real(pReal),intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness
real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress
real(pReal),intent(out), dimension(3,3,grid(1),grid(2),grid3) :: P !< PK stress
integer(pInt) :: &
calcMode, & !< CPFEM mode for calculation
j,k
real(pReal), dimension(3,3,3,3) :: max_dPdF, min_dPdF
real(pReal) :: max_dPdF_norm, min_dPdF_norm, defgradDetMin, defgradDetMax, defgradDet
PetscErrorCode :: ierr
external :: &
MPI_Allreduce
if (worldrank == 0_pInt) then
if (worldrank == 0_pInt) then
write(6,'(/,a)') ' ... evaluating constitutive response ......................................'
flush(6)
endif
calcMode = CPFEM_CALCRESULTS
if (forwardData) then ! aging results
calcMode = ior(calcMode, CPFEM_AGERESULTS)
calcMode = ior(calcMode, CPFEM_AGERESULTS)
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3])
endif
if (cutBack) then ! restore saved variables
calcMode = iand(calcMode, not(CPFEM_AGERESULTS))
calcMode = iand(calcMode, not(CPFEM_AGERESULTS))
endif
call CPFEM_general(CPFEM_COLLECT,F_lastInc(1:3,1:3,1,1,1),F(1:3,1:3,1,1,1), &
@ -994,11 +995,11 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc,&
do j = 1_pInt, product(grid(1:2))*grid3
defgradDet = math_det33(materialpoint_F(1:3,1:3,1,j))
defgradDetMax = max(defgradDetMax,defgradDet)
defgradDetMin = min(defgradDetMin,defgradDet)
defgradDetMin = min(defgradDetMin,defgradDet)
end do
call MPI_reduce(MPI_IN_PLACE,defgradDetMax,1,MPI_DOUBLE,MPI_MAX,0,PETSC_COMM_WORLD,ierr)
call MPI_reduce(MPI_IN_PLACE,defgradDetMin,1,MPI_DOUBLE,MPI_MIN,0,PETSC_COMM_WORLD,ierr)
if (worldrank == 0_pInt) then
if (worldrank == 0_pInt) then
write(6,'(a,1x,es11.4)') ' max determinant of deformation =', defgradDetMax
write(6,'(a,1x,es11.4)') ' min determinant of deformation =', defgradDetMin
flush(6)
@ -1016,26 +1017,28 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc,&
if (max_dPdF_norm < sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k)**2.0_pReal)) then
max_dPdF = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k)
max_dPdF_norm = sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k)**2.0_pReal)
endif
endif
if (min_dPdF_norm > sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k)**2.0_pReal)) then
min_dPdF = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k)
min_dPdF_norm = sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k)**2.0_pReal)
endif
endif
end do
call MPI_Allreduce(MPI_IN_PLACE,max_dPdF,81,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr)
call MPI_Allreduce(MPI_IN_PLACE,min_dPdF,81,MPI_DOUBLE,MPI_MIN,PETSC_COMM_WORLD,ierr)
C_minmaxAvg = 0.5_pReal*(max_dPdF + min_dPdF)
C_minmaxAvg = 0.5_pReal*(max_dPdF + min_dPdF)
C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) * wgt
call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
call debug_info()
restartWrite = .false. ! reset restartWrite status
cutBack = .false. ! reset cutBack status
P = reshape(materialpoint_P, [3,3,grid(1),grid(2),grid3])
P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P
P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P
call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
if (debugRotation .and. worldrank == 0_pInt) &
write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress (lab) / MPa =',&
@ -1057,7 +1060,7 @@ pure function utilities_calculateRate(avRate,timeinc_old,guess,field_lastInc,fie
use mesh, only: &
grid3, &
grid
implicit none
real(pReal), intent(in), dimension(3,3) :: avRate !< homogeneous addon
real(pReal), intent(in) :: &
@ -1069,7 +1072,7 @@ pure function utilities_calculateRate(avRate,timeinc_old,guess,field_lastInc,fie
field !< data of current step
real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: &
utilities_calculateRate
if (guess) then
utilities_calculateRate = (field-field_lastInc) / timeinc_old
else
@ -1080,16 +1083,16 @@ end function utilities_calculateRate
!--------------------------------------------------------------------------------------------------
!> @brief forwards a field with a pointwise given rate, if aim is given,
!> @brief forwards a field with a pointwise given rate, if aim is given,
!> ensures that the average matches the aim
!--------------------------------------------------------------------------------------------------
function utilities_forwardField(timeinc,field_lastInc,rate,aim)
use mesh, only: &
grid3, &
grid
implicit none
real(pReal), intent(in) :: &
real(pReal), intent(in) :: &
timeinc !< timeinc of current step
real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: &
field_lastInc, & !< initial field
@ -1100,10 +1103,10 @@ function utilities_forwardField(timeinc,field_lastInc,rate,aim)
utilities_forwardField
real(pReal), dimension(3,3) :: fieldDiff !< <a + adot*t> - aim
PetscErrorCode :: ierr
external :: &
MPI_Allreduce
utilities_forwardField = field_lastInc + rate*timeinc
if (present(aim)) then !< correct to match average
fieldDiff = sum(sum(sum(utilities_forwardField,dim=5),dim=4),dim=3)*wgt
@ -1158,7 +1161,7 @@ subroutine utilities_updateIPcoords(F)
grid3, &
grid3Offset, &
geomSize, &
mesh_ipCoordinates
mesh_ipCoordinates
implicit none
real(pReal), dimension(3,3,grid(1),grid(2),grid3), intent(in) :: F
@ -1175,7 +1178,7 @@ subroutine utilities_updateIPcoords(F)
integrator = geomSize * 0.5_pReal / PI
step = geomSize/real(grid, pReal)
!--------------------------------------------------------------------------------------------------
! average F
if (grid3Offset == 0_pInt) Favg = real(tensorField_fourier(1:3,1:3,1,1,1),pReal)*wgt

18
code/crystallite.f90
View File

@ -438,20 +438,16 @@ subroutine crystallite_init
call crystallite_orientations()
crystallite_orientation0 = crystallite_orientation ! store initial orientations for calculation of grain rotations
!***some debugging statement here
!write(6,*) 'CZ: before crystallite initialization'
!$OMP PARALLEL DO PRIVATE(myNgrains)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1_pInt,myNgrains
!***dirty way to pass orientation to constitutive module
call constitutive_microstructure( &
crystallite_orientation, &
crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fp(1:3,1:3,g,i,e), &
g,i,e) ! update dependent state variables to be consistent with basic states
call constitutive_microstructure(crystallite_orientation, &
crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fp(1:3,1:3,g,i,e), &
g,i,e) ! update dependent state variables to be consistent with basic states
enddo
enddo
enddo
@ -654,8 +650,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
if (iand(debug_level(debug_crystallite),debug_levelSelective) /= 0_pInt &
.and. FEsolving_execElem(1) <= debug_e &
.and. debug_e <= FEsolving_execElem(2)) then
write(6,'(/,a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> values at el (elFE) ip g ', &
write(6,'(/,a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> boundary values at el (elFE) ip g ', &
debug_e,'(',mesh_element(1,debug_e), ')',debug_i, debug_g
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> F ', &
math_transpose33(crystallite_partionedF(1:3,1:3,debug_g,debug_i,debug_e))
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> F0 ', &
math_transpose33(crystallite_partionedF0(1:3,1:3,debug_g,debug_i,debug_e))
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Fp0', &
@ -666,8 +664,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
math_transpose33(crystallite_partionedLp0(1:3,1:3,debug_g,debug_i,debug_e))
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Li0', &
math_transpose33(crystallite_partionedLi0(1:3,1:3,debug_g,debug_i,debug_e))
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> F ', &
math_transpose33(crystallite_partionedF(1:3,1:3,debug_g,debug_i,debug_e))
endif
!--------------------------------------------------------------------------------------------------

159
code/homogenization.f90
View File

@ -4,7 +4,7 @@
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @author Denny Tjahjanto, Max-Planck-Institut für Eisenforschung GmbH
!> @brief homogenization manager, organizing deformation partitioning and stress homogenization
!> @brief homogenization manager, organizing deformation partitioning and stress homogenization
!--------------------------------------------------------------------------------------------------
module homogenization
use prec, only: &
@ -91,7 +91,7 @@ subroutine homogenization_init
use mesh, only: &
mesh_maxNips, &
mesh_NcpElems, &
mesh_element, &
mesh_element, &
FE_Nips, &
FE_geomtype
#ifdef FEM
@ -151,7 +151,7 @@ subroutine homogenization_init
if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) &
call homogenization_RGC_init(FILEUNIT)
!--------------------------------------------------------------------------------------------------
! parse thermal from config file
call IO_checkAndRewind(FILEUNIT)
@ -162,7 +162,7 @@ subroutine homogenization_init
if (any(thermal_type == THERMAL_conduction_ID)) &
call thermal_conduction_init(FILEUNIT)
!--------------------------------------------------------------------------------------------------
! parse damage from config file
call IO_checkAndRewind(FILEUNIT)
@ -227,7 +227,7 @@ subroutine homogenization_init
thisSize => homogenization_RGC_sizePostResult
case default
knownHomogenization = .false.
end select
end select
write(FILEUNIT,'(/,a,/)') '['//trim(homogenization_name(p))//']'
if (knownHomogenization) then
write(FILEUNIT,'(a)') '(type)'//char(9)//trim(outputName)
@ -236,8 +236,8 @@ subroutine homogenization_init
do e = 1,thisNoutput(i)
write(FILEUNIT,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i)
enddo
endif
endif
endif
endif
i = thermal_typeInstance(p) ! which instance of this thermal type
knownThermal = .true. ! assume valid
select case(thermal_type(p)) ! split per thermal type
@ -258,15 +258,15 @@ subroutine homogenization_init
thisSize => thermal_conduction_sizePostResult
case default
knownThermal = .false.
end select
end select
if (knownThermal) then
write(FILEUNIT,'(a)') '(thermal)'//char(9)//trim(outputName)
if (thermal_type(p) /= THERMAL_isothermal_ID) then
do e = 1,thisNoutput(i)
write(FILEUNIT,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i)
enddo
endif
endif
endif
endif
i = damage_typeInstance(p) ! which instance of this damage type
knownDamage = .true. ! assume valid
select case(damage_type(p)) ! split per damage type
@ -287,15 +287,15 @@ subroutine homogenization_init
thisSize => damage_nonlocal_sizePostResult
case default
knownDamage = .false.
end select
end select
if (knownDamage) then
write(FILEUNIT,'(a)') '(damage)'//char(9)//trim(outputName)
if (damage_type(p) /= DAMAGE_none_ID) then
do e = 1,thisNoutput(i)
write(FILEUNIT,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i)
enddo
endif
endif
endif
endif
i = vacancyflux_typeInstance(p) ! which instance of this vacancy flux type
knownVacancyflux = .true. ! assume valid
select case(vacancyflux_type(p)) ! split per vacancy flux type
@ -316,15 +316,15 @@ subroutine homogenization_init
thisSize => vacancyflux_cahnhilliard_sizePostResult
case default
knownVacancyflux = .false.
end select
end select
if (knownVacancyflux) then
write(FILEUNIT,'(a)') '(vacancyflux)'//char(9)//trim(outputName)
if (vacancyflux_type(p) /= VACANCYFLUX_isoconc_ID) then
do e = 1,thisNoutput(i)
write(FILEUNIT,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i)
enddo
endif
endif
endif
endif
i = porosity_typeInstance(p) ! which instance of this porosity type
knownPorosity = .true. ! assume valid
select case(porosity_type(p)) ! split per porosity type
@ -340,15 +340,15 @@ subroutine homogenization_init
thisSize => porosity_phasefield_sizePostResult
case default
knownPorosity = .false.
end select
end select
if (knownPorosity) then
write(FILEUNIT,'(a)') '(porosity)'//char(9)//trim(outputName)
if (porosity_type(p) /= POROSITY_none_ID) then
do e = 1,thisNoutput(i)
write(FILEUNIT,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i)
enddo
endif
endif
endif
endif
i = hydrogenflux_typeInstance(p) ! which instance of this hydrogen flux type
knownHydrogenflux = .true. ! assume valid
select case(hydrogenflux_type(p)) ! split per hydrogen flux type
@ -364,15 +364,15 @@ subroutine homogenization_init
thisSize => hydrogenflux_cahnhilliard_sizePostResult
case default
knownHydrogenflux = .false.
end select
end select
if (knownHydrogenflux) then
write(FILEUNIT,'(a)') '(hydrogenflux)'//char(9)//trim(outputName)
if (hydrogenflux_type(p) /= HYDROGENFLUX_isoconc_ID) then
do e = 1,thisNoutput(i)
write(FILEUNIT,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i)
enddo
endif
endif
endif
endif
endif
enddo
close(FILEUNIT)
@ -421,13 +421,13 @@ subroutine homogenization_init
vacancyflux_maxSizePostResults = max(vacancyflux_maxSizePostResults ,vacancyfluxState (p)%sizePostResults)
porosity_maxSizePostResults = max(porosity_maxSizePostResults ,porosityState (p)%sizePostResults)
hydrogenflux_maxSizePostResults = max(hydrogenflux_maxSizePostResults ,hydrogenfluxState(p)%sizePostResults)
enddo
enddo
#ifdef FEM
allocate(homogOutput (material_Nhomogenization ))
allocate(crystalliteOutput(material_Ncrystallite, homogenization_maxNgrains))
allocate(phaseOutput (material_Nphase, homogenization_maxNgrains))
do p = 1, material_Nhomogenization
do p = 1, material_Nhomogenization
homogOutput(p)%sizeResults = homogState (p)%sizePostResults + &
thermalState (p)%sizePostResults + &
damageState (p)%sizePostResults + &
@ -436,19 +436,19 @@ subroutine homogenization_init
hydrogenfluxState(p)%sizePostResults
homogOutput(p)%sizeIpCells = count(material_homog==p)
allocate(homogOutput(p)%output(homogOutput(p)%sizeResults,homogOutput(p)%sizeIpCells))
enddo
enddo
do p = 1, material_Ncrystallite; do e = 1, homogenization_maxNgrains
crystalliteOutput(p,e)%sizeResults = crystallite_sizePostResults(p)
crystalliteOutput(p,e)%sizeIpCells = count(microstructure_crystallite(mesh_element(4,:)) == p .and. &
homogenization_Ngrains (mesh_element(3,:)) >= e)*mesh_maxNips
allocate(crystalliteOutput(p,e)%output(crystalliteOutput(p,e)%sizeResults,crystalliteOutput(p,e)%sizeIpCells))
enddo; enddo
enddo; enddo
do p = 1, material_Nphase; do e = 1, homogenization_maxNgrains
phaseOutput(p,e)%sizeResults = plasticState (p)%sizePostResults + &
sum(sourceState (p)%p(:)%sizePostResults)
phaseOutput(p,e)%sizeIpCells = count(material_phase(e,:,:) == p)
allocate(phaseOutput(p,e)%output(phaseOutput(p,e)%sizeResults,phaseOutput(p,e)%sizeIpCells))
enddo; enddo
enddo; enddo
#else
materialpoint_sizeResults = 1 & ! grain count
+ 1 + homogenization_maxSizePostResults & ! homogSize & homogResult
@ -459,11 +459,11 @@ subroutine homogenization_init
+ hydrogenflux_maxSizePostResults &
+ homogenization_maxNgrains * (1 + crystallite_maxSizePostResults & ! crystallite size & crystallite results
+ 1 + constitutive_plasticity_maxSizePostResults & ! constitutive size & constitutive results
+ constitutive_source_maxSizePostResults)
+ constitutive_source_maxSizePostResults)
allocate(materialpoint_results(materialpoint_sizeResults,mesh_maxNips,mesh_NcpElems))
#endif
mainProcess: if (worldrank == 0) then
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- homogenization init -+>>>'
write(6,'(a)') ' $Id$'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
@ -494,10 +494,10 @@ subroutine homogenization_init
write(6,'(a32,1x,7(i8,1x))') 'maxSizePostResults: ', homogenization_maxSizePostResults
endif
flush(6)
if (debug_g < 1 .or. debug_g > homogenization_Ngrains(mesh_element(3,debug_e))) &
call IO_error(602_pInt,ext_msg='component (grain)')
end subroutine homogenization_init
@ -511,7 +511,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
stepIncreaseHomog, &
nHomog, &
nMPstate
use math, only: &
use math, only: &
math_transpose33
use FEsolving, only: &
FEsolving_execElem, &
@ -529,11 +529,11 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
porosityState, &
hydrogenfluxState, &
phase_Nsources, &
mappingHomogenization, &
mappingHomogenization, &
mappingConstitutive, &
homogenization_Ngrains
use crystallite, only: &
crystallite_F0, &
crystallite_Fp0, &
@ -572,7 +572,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
debug_MaterialpointStateLoopDistribution
use math, only: &
math_pDecomposition
implicit none
real(pReal), intent(in) :: dt !< time increment
logical, intent(in) :: updateJaco !< initiating Jacobian update
@ -609,7 +609,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
do mySource = 1_pInt, phase_Nsources(mappingConstitutive(2,g,i,e))
sourceState(mappingConstitutive(2,g,i,e))%p(mySource)%partionedState0(:,mappingConstitutive(1,g,i,e)) = &
sourceState(mappingConstitutive(2,g,i,e))%p(mySource)%state0( :,mappingConstitutive(1,g,i,e))
enddo
enddo
crystallite_partionedFp0(1:3,1:3,g,i,e) = crystallite_Fp0(1:3,1:3,g,i,e) ! ...plastic def grads
crystallite_partionedLp0(1:3,1:3,g,i,e) = crystallite_Lp0(1:3,1:3,g,i,e) ! ...plastic velocity grads
@ -653,7 +653,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
hydrogenfluxState(mappingHomogenization(2,i,e))%State0( :,mappingHomogenization(1,i,e)) ! ...internal hydrogen transport state
enddo
NiterationHomog = 0_pInt
cutBackLooping: do while (.not. terminallyIll .and. &
any(materialpoint_subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMinHomog))
@ -661,43 +661,50 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(mesh_element(3,e))
IpLooping1: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
converged: if ( materialpoint_converged(i,e) ) then
#ifndef _OPENMP
if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i) &
.and. ((e == debug_e .and. i == debug_i) &
.or. .not. iand(debug_level(debug_homogenization),debug_levelSelective) /= 0_pInt)) then
write(6,'(a,1x,f12.8,1x,a,1x,f12.8,1x,a,i8,1x,i2/)') '<< HOMOG >> winding forward from', &
materialpoint_subFrac(i,e), 'to current materialpoint_subFrac', &
materialpoint_subFrac(i,e)+materialpoint_subStep(i,e),'in materialpoint_stressAndItsTangent at el ip',e,i
endif
#endif
!--------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
! calculate new subStep and new subFrac
materialpoint_subFrac(i,e) = materialpoint_subFrac(i,e) + materialpoint_subStep(i,e)
!$OMP FLUSH(materialpoint_subFrac)
materialpoint_subStep(i,e) = min(1.0_pReal-materialpoint_subFrac(i,e), &
stepIncreaseHomog*materialpoint_subStep(i,e)) ! introduce flexibility for step increase/acceleration
!$OMP FLUSH(materialpoint_subStep)
steppingNeeded: if (materialpoint_subStep(i,e) > subStepMinHomog) then
! wind forward grain starting point of...
crystallite_partionedF0(1:3,1:3,1:myNgrains,i,e) = &
crystallite_partionedF(1:3,1:3,1:myNgrains,i,e) ! ...def grads
crystallite_partionedFp0(1:3,1:3,1:myNgrains,i,e) = &
crystallite_Fp(1:3,1:3,1:myNgrains,i,e) ! ...plastic def grads
crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e) = &
crystallite_Lp(1:3,1:3,1:myNgrains,i,e) ! ...plastic velocity grads
crystallite_partionedFi0(1:3,1:3,1:myNgrains,i,e) = &
crystallite_Fi(1:3,1:3,1:myNgrains,i,e) ! ...intermediate def grads
crystallite_partionedLi0(1:3,1:3,1:myNgrains,i,e) = &
crystallite_Li(1:3,1:3,1:myNgrains,i,e) ! ...intermediate velocity grads
crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,1:myNgrains,i,e) = &
crystallite_dPdF(1:3,1:3,1:3,1:3,1:myNgrains,i,e) ! ...stiffness
crystallite_partionedTstar0_v(1:6,1:myNgrains,i,e) = &
crystallite_Tstar_v(1:6,1:myNgrains,i,e) ! ...2nd PK stress
do g = 1,myNgrains
plasticState (mappingConstitutive(2,g,i,e))%partionedState0(:,mappingConstitutive(1,g,i,e)) = &
plasticState (mappingConstitutive(2,g,i,e))%state( :,mappingConstitutive(1,g,i,e))
@ -705,7 +712,8 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
sourceState(mappingConstitutive(2,g,i,e))%p(mySource)%partionedState0(:,mappingConstitutive(1,g,i,e)) = &
sourceState(mappingConstitutive(2,g,i,e))%p(mySource)%state( :,mappingConstitutive(1,g,i,e))
enddo
enddo
enddo
forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), &
homogState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) &
homogState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) = &
@ -757,17 +765,17 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
else ! cutback makes sense
materialpoint_subStep(i,e) = subStepSizeHomog * materialpoint_subStep(i,e) ! crystallite had severe trouble, so do a significant cutback
!$OMP FLUSH(materialpoint_subStep)
#ifndef _OPENMP
if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i) &
.or. .not. iand(debug_level(debug_homogenization), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,1x,f12.8,a,i8,1x,i2/)') &
'<< HOMOG >> cutback step in materialpoint_stressAndItsTangent with new materialpoint_subStep:',&
materialpoint_subStep(i,e),' at el ip',e,i
materialpoint_subStep(i,e),' at el ip',e,i
endif
#endif
!--------------------------------------------------------------------------------------------------
! restore...
crystallite_Fp(1:3,1:3,1:myNgrains,i,e) = &
@ -789,7 +797,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
sourceState(mappingConstitutive(2,g,i,e))%p(mySource)%state( :,mappingConstitutive(1,g,i,e)) = &
sourceState(mappingConstitutive(2,g,i,e))%p(mySource)%partionedState0(:,mappingConstitutive(1,g,i,e))
enddo
enddo
enddo
forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), &
homogState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) &
homogState(mappingHomogenization(2,i,e))%State( :,mappingHomogenization(1,i,e)) = &
@ -814,9 +822,9 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
hydrogenfluxState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) &
hydrogenfluxState(mappingHomogenization(2,i,e))%State( :,mappingHomogenization(1,i,e)) = &
hydrogenfluxState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e))! ...internal hydrogen transport state
endif
endif
endif converged
if (materialpoint_subStep(i,e) > subStepMinHomog) then
materialpoint_requested(i,e) = .true.
materialpoint_subF(1:3,1:3,i,e) = materialpoint_subF0(1:3,1:3,i,e) + &
@ -829,7 +837,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
!$OMP END PARALLEL DO
NiterationMPstate = 0_pInt
convergenceLooping: do while (.not. terminallyIll .and. &
any( materialpoint_requested(:,FEsolving_execELem(1):FEsolving_execElem(2)) &
.and. .not. materialpoint_doneAndHappy(1,:,FEsolving_execELem(1):FEsolving_execElem(2)) &
@ -839,7 +847,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
!--------------------------------------------------------------------------------------------------
! deformation partitioning
! based on materialpoint_subF0,.._subF,crystallite_partionedF0, and homogenization_state,
! based on materialpoint_subF0,.._subF,crystallite_partionedF0, and homogenization_state,
! results in crystallite_partionedF
!$OMP PARALLEL DO PRIVATE(myNgrains)
elementLooping2: do e = FEsolving_execElem(1),FEsolving_execElem(2)
@ -856,7 +864,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
enddo IpLooping2
enddo elementLooping2
!$OMP END PARALLEL DO
!--------------------------------------------------------------------------------------------------
! crystallite integration
! based on crystallite_partionedF0,.._partionedF
@ -897,7 +905,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
enddo cutBackLooping
if (.not. terminallyIll ) then
if (.not. terminallyIll ) then
call crystallite_orientations() ! calculate crystal orientations
!$OMP PARALLEL DO
elementLooping4: do e = FEsolving_execElem(1),FEsolving_execElem(2)
@ -911,7 +919,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
write(6,'(/,a,/)') '<< HOMOG >> Material Point terminally ill'
!$OMP END CRITICAL (write2out)
endif
end subroutine materialpoint_stressAndItsTangent
@ -927,10 +935,10 @@ subroutine materialpoint_postResults
use material, only: &
mappingHomogenization, &
#ifdef FEM
mappingConstitutive, &
mappingConstitutive, &
homogenization_maxNgrains, &
material_Ncrystallite, &
material_Nphase, &
material_Nphase, &
#else
homogState, &
thermalState, &
@ -971,10 +979,10 @@ subroutine materialpoint_postResults
myHomog, &
myPhase, &
crystalliteCtr(material_Ncrystallite, homogenization_maxNgrains), &
phaseCtr (material_Nphase, homogenization_maxNgrains)
phaseCtr (material_Nphase, homogenization_maxNgrains)
real(pReal), dimension(1+crystallite_maxSizePostResults + &
1+constitutive_plasticity_maxSizePostResults + &
constitutive_source_maxSizePostResults) :: &
constitutive_source_maxSizePostResults) :: &
crystalliteResults
@ -989,7 +997,7 @@ subroutine materialpoint_postResults
homogOutput(myHomog)%output(1: &
homogOutput(myHomog)%sizeResults, &
thePos) = homogenization_postResults(i,e)
grainLooping :do g = 1,myNgrains
myPhase = mappingConstitutive(2,g,i,e)
crystalliteResults(1:1+crystallite_sizePostResults(myCrystallite) + &
@ -999,16 +1007,16 @@ subroutine materialpoint_postResults
homogenization_Ngrains (mesh_element(3,e)) >= g) then
crystalliteCtr(myCrystallite,g) = crystalliteCtr(myCrystallite,g) + 1_pInt
crystalliteOutput(myCrystallite,g)% &
output(1:crystalliteOutput(myCrystallite,g)%sizeResults,crystalliteCtr(myCrystallite,g)) = &
output(1:crystalliteOutput(myCrystallite,g)%sizeResults,crystalliteCtr(myCrystallite,g)) = &
crystalliteResults(2:1+crystalliteOutput(myCrystallite,g)%sizeResults)
endif
if (material_phase(g,i,e) == myPhase) then
phaseCtr(myPhase,g) = phaseCtr(myPhase,g) + 1_pInt
phaseOutput(myPhase,g)% &
output(1:phaseOutput(myPhase,g)%sizeResults,phaseCtr(myPhase,g)) = &
output(1:phaseOutput(myPhase,g)%sizeResults,phaseCtr(myPhase,g)) = &
crystalliteResults(3 + crystalliteOutput(myCrystallite,g)%sizeResults: &
1 + crystalliteOutput(myCrystallite,g)%sizeResults + &
1 + plasticState (myphase)%sizePostResults + &
1 + crystalliteOutput(myCrystallite,g)%sizeResults + &
1 + plasticState (myphase)%sizePostResults + &
sum(sourceState(myphase)%p(:)%sizePostResults))
endif
enddo grainLooping
@ -1022,7 +1030,7 @@ subroutine materialpoint_postResults
myCrystallite = microstructure_crystallite(mesh_element(4,e))
IpLooping: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
thePos = 0_pInt
theSize = homogState (mappingHomogenization(2,i,e))%sizePostResults &
+ thermalState (mappingHomogenization(2,i,e))%sizePostResults &
+ damageState (mappingHomogenization(2,i,e))%sizePostResults &
@ -1053,8 +1061,8 @@ subroutine materialpoint_postResults
#endif
end subroutine materialpoint_postResults
!--------------------------------------------------------------------------------------------------
!> @brief partition material point def grad onto constituents
!--------------------------------------------------------------------------------------------------
@ -1103,7 +1111,7 @@ end subroutine homogenization_partitionDeformation
!--------------------------------------------------------------------------------------------------
!> @brief update the internal state of the homogenization scheme and tell whether "done" and
!> @brief update the internal state of the homogenization scheme and tell whether "done" and
!> "happy" with result
!--------------------------------------------------------------------------------------------------
function homogenization_updateState(ip,el)
@ -1138,7 +1146,7 @@ function homogenization_updateState(ip,el)
ip, & !< integration point
el !< element number
logical, dimension(2) :: homogenization_updateState
homogenization_updateState = .true.
chosenHomogenization: select case(homogenization_type(mesh_element(3,el)))
case (HOMOGENIZATION_RGC_ID) chosenHomogenization
@ -1219,12 +1227,13 @@ subroutine homogenization_averageStressAndItsTangent(ip,el)
integer(pInt), intent(in) :: &
ip, & !< integration point
el !< element number
chosenHomogenization: select case(homogenization_type(mesh_element(3,el)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization
materialpoint_P(1:3,1:3,ip,el) = sum(crystallite_P(1:3,1:3,1:1,ip,el),3)
materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el) &
= sum(crystallite_dPdF(1:3,1:3,1:3,1:3,1:1,ip,el),5)
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
call homogenization_isostrain_averageStressAndItsTangent(&
materialpoint_P(1:3,1:3,ip,el), &
@ -1244,7 +1253,7 @@ subroutine homogenization_averageStressAndItsTangent(ip,el)
end subroutine homogenization_averageStressAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief return array of homogenization results for post file inclusion. call only,
!> @brief return array of homogenization results for post file inclusion. call only,
!> if homogenization_sizePostResults(i,e) > 0 !!
!--------------------------------------------------------------------------------------------------
function homogenization_postResults(ip,el)
@ -1300,7 +1309,7 @@ function homogenization_postResults(ip,el)
porosity_phasefield_postResults
use hydrogenflux_cahnhilliard, only: &
hydrogenflux_cahnhilliard_postResults
implicit none
integer(pInt), intent(in) :: &
ip, & !< integration point
@ -1314,7 +1323,7 @@ function homogenization_postResults(ip,el)
homogenization_postResults
integer(pInt) :: &
startPos, endPos
homogenization_postResults = 0.0_pReal
startPos = 1_pInt

13
code/plastic_phenoplus.f90
View File

@ -5,7 +5,7 @@
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @author Chen Zhang, Michigan State University
!> @brief material subroutine for phenomenological crystal plasticity formulation using a powerlaw
!! fitting
!... fitting
!--------------------------------------------------------------------------------------------------
module plastic_phenoplus
use prec, only: &
@ -830,7 +830,7 @@ subroutine plastic_phenoplus_microstructure(orientation,ipc,ip,el)
ns = plastic_phenoplus_totalNslip(instance)
nt = plastic_phenoplus_totalNtwin(instance)
offset_acshear_slip = ns + nt + 2_pInt
kappa_max = ns + nt + 2_pInt + ns + nt !location of kappa in plasticState
index_kappa = ns + nt + 2_pInt + ns + nt !location of kappa in plasticState
!***gather my accumulative shear from palsticState
findMyShear: do j = 1_pInt,ns
@ -876,6 +876,7 @@ subroutine plastic_phenoplus_microstructure(orientation,ipc,ip,el)
2.0_pReal * &
(kappa_max - 1.0_pReal) * &
(1.0_pReal - mprimeavg)
enddo loopMySlip
end subroutine plastic_phenoplus_microstructure
@ -980,6 +981,14 @@ subroutine plastic_phenoplus_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
((abs(tau_slip_neg)/(plasticState(ph)%state(j, of)* &
plasticState(ph)%state(j+index_kappa, of))) & !?should we make it direction aware
**plastic_phenoplus_n_slip(instance))*sign(1.0_pReal,tau_slip_neg)
!***in case for future use
! gdot_slip_pos = 0.5_pReal*plastic_phenoplus_gdot0_slip(instance)* &
! ((abs(tau_slip_pos)/(plasticState(ph)%state(j, of))) & !in-place modification of gdot
! **plastic_phenoplus_n_slip(instance))*sign(1.0_pReal,tau_slip_pos)
! gdot_slip_neg = 0.5_pReal*plastic_phenoplus_gdot0_slip(instance)* &
! ((abs(tau_slip_neg)/(plasticState(ph)%state(j, of))) & !?should we make it direction aware
! **plastic_phenoplus_n_slip(instance))*sign(1.0_pReal,tau_slip_neg)
Lp = Lp + (1.0_pReal-plasticState(ph)%state(index_F,of))*& ! 1-F
(gdot_slip_pos+gdot_slip_neg)*lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph)