Coupled phase field implementation for the basic petsc solver. needs clearing up… work in progress
This commit is contained in:
parent
a05d99a54c
commit
724380c624
|
@ -523,7 +523,6 @@ program DAMASK_spectral_Driver
|
||||||
' increment ', totalIncsCounter, ' NOT converged'
|
' increment ', totalIncsCounter, ' NOT converged'
|
||||||
notConvergedCounter = notConvergedCounter + 1_pInt
|
notConvergedCounter = notConvergedCounter + 1_pInt
|
||||||
endif; flush(6)
|
endif; flush(6)
|
||||||
call utilities_temperatureUpdate()
|
|
||||||
if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency
|
if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency
|
||||||
write(6,'(1/,a)') ' ... writing results to file ......................................'
|
write(6,'(1/,a)') ' ... writing results to file ......................................'
|
||||||
write(resUnit) materialpoint_results ! write result to file
|
write(resUnit) materialpoint_results ! write result to file
|
||||||
|
|
|
@ -62,6 +62,10 @@ module DAMASK_spectral_SolverBasicPETSc
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! common pointwise data
|
! common pointwise data
|
||||||
real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, Fdot, F_lastInc2
|
real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, Fdot, F_lastInc2
|
||||||
|
real(pReal), private, dimension(:,:,:), allocatable :: &
|
||||||
|
heatSource_lastInc, &
|
||||||
|
phaseField_lastInc, &
|
||||||
|
phaseFieldDot
|
||||||
complex(pReal), private, dimension(:,:,:,:,:), allocatable :: inertiaField_fourier
|
complex(pReal), private, dimension(:,:,:,:,:), allocatable :: inertiaField_fourier
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
|
@ -78,7 +82,7 @@ module DAMASK_spectral_SolverBasicPETSc
|
||||||
C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness
|
C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness
|
||||||
C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 stiffness
|
C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 stiffness
|
||||||
S = 0.0_pReal !< current compliance (filled up with zeros)
|
S = 0.0_pReal !< current compliance (filled up with zeros)
|
||||||
real(pReal), private :: err_stress, err_div, err_divPrev, err_divDummy
|
real(pReal), private :: err_stress, err_div, err_divPrev, err_divDummy, err_phaseField, phaseField_Avg
|
||||||
logical, private :: ForwardData
|
logical, private :: ForwardData
|
||||||
integer(pInt), private :: &
|
integer(pInt), private :: &
|
||||||
totalIter = 0_pInt !< total iteration in current increment
|
totalIter = 0_pInt !< total iteration in current increment
|
||||||
|
@ -144,8 +148,10 @@ subroutine basicPETSc_init(temperature)
|
||||||
temperature
|
temperature
|
||||||
#include <finclude/petscdmda.h90>
|
#include <finclude/petscdmda.h90>
|
||||||
#include <finclude/petscsnes.h90>
|
#include <finclude/petscsnes.h90>
|
||||||
|
#include <finclude/petscvec.h>
|
||||||
real(pReal), dimension(:,:,:,:,:), allocatable :: P
|
real(pReal), dimension(:,:,:,:,:), allocatable :: P
|
||||||
PetscScalar, dimension(:,:,:,:), pointer :: F
|
PetscScalar, dimension(:,:,:,:), pointer :: xx_psc, F
|
||||||
|
PetscScalar, dimension(:,:,:), pointer :: phaseField
|
||||||
PetscErrorCode :: ierr
|
PetscErrorCode :: ierr
|
||||||
PetscObject :: dummy
|
PetscObject :: dummy
|
||||||
real(pReal), dimension(3,3) :: &
|
real(pReal), dimension(3,3) :: &
|
||||||
|
@ -167,6 +173,9 @@ subroutine basicPETSc_init(temperature)
|
||||||
allocate (F_lastInc2(3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
|
allocate (F_lastInc2(3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
|
||||||
allocate (Fdot (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
|
allocate (Fdot (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
|
||||||
allocate (inertiaField_fourier (grid1Red,grid(2),grid(3),3,3),source = cmplx(0.0_pReal,0.0_pReal,pReal))
|
allocate (inertiaField_fourier (grid1Red,grid(2),grid(3),3,3),source = cmplx(0.0_pReal,0.0_pReal,pReal))
|
||||||
|
allocate (heatSource_lastInc (grid(1),grid(2),grid(3)),source = 0.0_pReal)
|
||||||
|
allocate (phaseField_lastInc(grid(1),grid(2),grid(3)),source = 0.0_pReal)
|
||||||
|
allocate (phaseFieldDot (grid(1),grid(2),grid(3)),source = 0.0_pReal)
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! initialize solver specific parts of PETSc
|
! initialize solver specific parts of PETSc
|
||||||
|
@ -174,9 +183,10 @@ subroutine basicPETSc_init(temperature)
|
||||||
call DMDACreate3d(PETSC_COMM_WORLD, &
|
call DMDACreate3d(PETSC_COMM_WORLD, &
|
||||||
DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, &
|
DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, &
|
||||||
DMDA_STENCIL_BOX,grid(1),grid(2),grid(3),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, &
|
DMDA_STENCIL_BOX,grid(1),grid(2),grid(3),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, &
|
||||||
9,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)
|
10,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)
|
||||||
CHKERRQ(ierr)
|
CHKERRQ(ierr)
|
||||||
call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr)
|
call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr)
|
||||||
|
!call DMDASNESSetFunctionLocal(da,INSERT_VALUES,BasicPETSC_formResidual,dummy,ierr); CHKERRQ(ierr) ! needed for newer versions of petsc
|
||||||
call DMDASetLocalFunction(da,BasicPETSC_formResidual,ierr); CHKERRQ(ierr)
|
call DMDASetLocalFunction(da,BasicPETSC_formResidual,ierr); CHKERRQ(ierr)
|
||||||
call SNESSetDM(snes,da,ierr); CHKERRQ(ierr)
|
call SNESSetDM(snes,da,ierr); CHKERRQ(ierr)
|
||||||
call SNESSetConvergenceTest(snes,BasicPETSC_converged,dummy,PETSC_NULL_FUNCTION,ierr)
|
call SNESSetConvergenceTest(snes,BasicPETSC_converged,dummy,PETSC_NULL_FUNCTION,ierr)
|
||||||
|
@ -188,12 +198,16 @@ subroutine basicPETSc_init(temperature)
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! init fields
|
! init fields
|
||||||
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with
|
call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with
|
||||||
|
F => xx_psc(0:8,:,:,:)
|
||||||
|
phaseField => xx_psc(9,:,:,:)
|
||||||
|
|
||||||
if (restartInc == 1_pInt) then ! no deformation (no restart)
|
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
|
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid(3)) ! initialize to identity
|
||||||
F = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)])
|
F = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)])
|
||||||
F_lastInc2 = F_lastInc
|
F_lastInc2 = F_lastInc
|
||||||
|
phaseField = temperature
|
||||||
|
phaseField_lastInc = temperature
|
||||||
elseif (restartInc > 1_pInt) then ! using old values from file
|
elseif (restartInc > 1_pInt) then ! using old values from file
|
||||||
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) &
|
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) &
|
||||||
write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
|
write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
|
||||||
|
@ -227,10 +241,10 @@ subroutine basicPETSc_init(temperature)
|
||||||
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(&
|
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(&
|
||||||
F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
|
F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
|
||||||
call Utilities_constitutiveResponse(&
|
call Utilities_constitutiveResponse(&
|
||||||
reshape(F(0:8,0:grid(1)-1_pInt,0:grid(2)-1_pInt,0:grid(3)-1_pInt),[3,3,grid(1),grid(2),grid(3)]),&
|
reshape(F,[3,3,grid(1),grid(2),grid(3)]),&
|
||||||
reshape(F(0:8,0:grid(1)-1_pInt,0:grid(2)-1_pInt,0:grid(3)-1_pInt),[3,3,grid(1),grid(2),grid(3)]),&
|
reshape(F,[3,3,grid(1),grid(2),grid(3)]),&
|
||||||
temperature,0.0_pReal,P,C_volAvg,C_minmaxAvg,temp33_Real,.false.,math_I3)
|
temperature,0.0_pReal,P,C_volAvg,C_minmaxAvg,temp33_Real,.false.,math_I3)
|
||||||
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! write data back into PETSc
|
call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! write data back into PETSc
|
||||||
if (restartInc == 1_pInt) then ! use initial stiffness as reference stiffness
|
if (restartInc == 1_pInt) then ! use initial stiffness as reference stiffness
|
||||||
temp3333_Real = C_minMaxAvg
|
temp3333_Real = C_minMaxAvg
|
||||||
endif
|
endif
|
||||||
|
@ -267,6 +281,8 @@ type(tSolutionState) function basicPETSc_solution( &
|
||||||
use FEsolving, only: &
|
use FEsolving, only: &
|
||||||
restartWrite, &
|
restartWrite, &
|
||||||
terminallyIll
|
terminallyIll
|
||||||
|
use homogenization, only: &
|
||||||
|
materialpoint_heat
|
||||||
|
|
||||||
implicit none
|
implicit none
|
||||||
#include <finclude/petscdmda.h90>
|
#include <finclude/petscdmda.h90>
|
||||||
|
@ -291,12 +307,14 @@ type(tSolutionState) function basicPETSc_solution( &
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! PETSc Data
|
! PETSc Data
|
||||||
PetscScalar, pointer :: F(:,:,:,:)
|
PetscScalar, pointer :: xx_psc(:,:,:,:), F(:,:,:,:), phaseField(:,:,:)
|
||||||
PetscErrorCode :: ierr
|
PetscErrorCode :: ierr
|
||||||
SNESConvergedReason :: reason
|
SNESConvergedReason :: reason
|
||||||
incInfo = incInfoIn
|
incInfo = incInfoIn
|
||||||
|
|
||||||
call DMDAVecGetArrayF90(da,solution_vec,F,ierr)
|
call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with
|
||||||
|
F => xx_psc(0:8,:,:,:)
|
||||||
|
phaseField => xx_psc(9,:,:,:)
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! restart information for spectral solver
|
! restart information for spectral solver
|
||||||
if (restartWrite) then
|
if (restartWrite) then
|
||||||
|
@ -327,6 +345,7 @@ type(tSolutionState) function basicPETSc_solution( &
|
||||||
F_aim = F_aim_lastInc
|
F_aim = F_aim_lastInc
|
||||||
F = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)])
|
F = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)])
|
||||||
C_volAvg = C_volAvgLastInc
|
C_volAvg = C_volAvgLastInc
|
||||||
|
phaseField = phaseField_lastInc
|
||||||
else
|
else
|
||||||
C_volAvgLastInc = C_volAvg
|
C_volAvgLastInc = C_volAvg
|
||||||
|
|
||||||
|
@ -349,14 +368,18 @@ type(tSolutionState) function basicPETSc_solution( &
|
||||||
F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
|
F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
|
||||||
Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,params%rotation_BC), &
|
Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,params%rotation_BC), &
|
||||||
timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)]))
|
timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)]))
|
||||||
|
phaseFieldDot = (phaseField - phaseField_lastInc)/timeinc_old
|
||||||
F_lastInc2 = F_lastInc
|
F_lastInc2 = F_lastInc
|
||||||
F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid(3)])
|
F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid(3)])
|
||||||
|
phaseField_lastInc = phaseField
|
||||||
|
heatSource_lastInc = reshape(materialpoint_heat(1,:),[grid(1),grid(2),grid(3)])
|
||||||
endif
|
endif
|
||||||
F_aim = F_aim + f_aimDot * timeinc
|
F_aim = F_aim + f_aimDot * timeinc
|
||||||
|
|
||||||
F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot,math_rotate_backward33(F_aim, &
|
F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot,math_rotate_backward33(F_aim, &
|
||||||
rotation_BC)),[9,grid(1),grid(2),grid(3)])
|
rotation_BC)),[9,grid(1),grid(2),grid(3)])
|
||||||
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
|
phaseField = phaseField_lastInc + phaseFieldDot*timeinc
|
||||||
|
call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr)
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! update stiffness (and gamma operator)
|
! update stiffness (and gamma operator)
|
||||||
|
@ -412,26 +435,50 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
|
||||||
wgt, &
|
wgt, &
|
||||||
field_real, &
|
field_real, &
|
||||||
field_fourier, &
|
field_fourier, &
|
||||||
|
phaseField_real, &
|
||||||
|
phaseField_fourier, &
|
||||||
Utilities_FFTforward, &
|
Utilities_FFTforward, &
|
||||||
Utilities_FFTbackward, &
|
Utilities_FFTbackward, &
|
||||||
|
utilities_scalarFFTforward, &
|
||||||
|
utilities_scalarFFTbackward, &
|
||||||
Utilities_fourierConvolution, &
|
Utilities_fourierConvolution, &
|
||||||
Utilities_inverseLaplace, &
|
Utilities_inverseLaplace, &
|
||||||
|
Utilities_diffusion, &
|
||||||
Utilities_constitutiveResponse, &
|
Utilities_constitutiveResponse, &
|
||||||
Utilities_divergenceRMS
|
Utilities_divergenceRMS
|
||||||
use IO, only: &
|
use IO, only: &
|
||||||
IO_intOut
|
IO_intOut
|
||||||
|
use crystallite, only: &
|
||||||
|
crystallite_temperature
|
||||||
|
use homogenization, only: &
|
||||||
|
materialpoint_heat
|
||||||
|
|
||||||
implicit none
|
implicit none
|
||||||
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &
|
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &
|
||||||
in
|
in
|
||||||
PetscScalar, dimension(3,3,grid(1),grid(2),grid(3)) :: &
|
PetscScalar, target, dimension(10, &
|
||||||
x_scal, &
|
XG_RANGE,YG_RANGE,ZG_RANGE) :: &
|
||||||
|
x_scal
|
||||||
|
PetscScalar, target, dimension(10, &
|
||||||
|
X_RANGE,Y_RANGE,Z_RANGE) :: &
|
||||||
f_scal
|
f_scal
|
||||||
|
PetscScalar, pointer, dimension(:,:,:,:) :: &
|
||||||
|
F, &
|
||||||
|
residual_F
|
||||||
|
PetscScalar, pointer, dimension(:,:,:) :: &
|
||||||
|
phaseField, &
|
||||||
|
residual_phaseField
|
||||||
PetscInt :: &
|
PetscInt :: &
|
||||||
PETScIter, &
|
PETScIter, &
|
||||||
nfuncs
|
nfuncs
|
||||||
PetscObject :: dummy
|
PetscObject :: dummy
|
||||||
PetscErrorCode :: ierr
|
PetscErrorCode :: ierr
|
||||||
|
real(pReal) :: mobility, diffusivity
|
||||||
|
|
||||||
|
F => x_scal(1:9,1:grid(1),1:grid(2),1:grid(3))
|
||||||
|
residual_F => f_scal(1:9,1:grid(1),1:grid(2),1:grid(3))
|
||||||
|
phaseField => x_scal(10,1:grid(1),1:grid(2),1:grid(3))
|
||||||
|
residual_phaseField => f_scal(10,1:grid(1),1:grid(2),1:grid(3))
|
||||||
|
|
||||||
call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr)
|
call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr)
|
||||||
call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr)
|
call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr)
|
||||||
|
@ -453,24 +500,26 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! evaluate inertia
|
! evaluate inertia
|
||||||
if (params%density > 0.0_pReal) then
|
dynamic: if (params%density > 0.0_pReal) then
|
||||||
f_scal = ((x_scal - F_lastInc)/params%timeinc - (F_lastInc - F_lastInc2)/params%timeincOld)/&
|
residual_F = ((F - reshape(F_lastInc,[9,grid(1),grid(2),grid(3)]))/params%timeinc - &
|
||||||
|
reshape(F_lastInc - F_lastInc2, [9,grid(1),grid(2),grid(3)])/params%timeincOld)/&
|
||||||
((params%timeinc + params%timeincOld)/2.0_pReal)
|
((params%timeinc + params%timeincOld)/2.0_pReal)
|
||||||
f_scal = params%density*product(geomSize/grid)*f_scal
|
residual_F = params%density*product(geomSize/grid)*residual_F
|
||||||
field_real = 0.0_pReal
|
field_real = 0.0_pReal
|
||||||
field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3) = reshape(f_scal,[grid(1),grid(2),grid(3),3,3],&
|
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],&
|
||||||
order=[4,5,1,2,3]) ! field real has a different order
|
order=[4,5,1,2,3]) ! field real has a different order
|
||||||
call Utilities_FFTforward()
|
call Utilities_FFTforward()
|
||||||
call Utilities_inverseLaplace()
|
call Utilities_inverseLaplace()
|
||||||
inertiaField_fourier = field_fourier
|
inertiaField_fourier = field_fourier
|
||||||
else
|
else dynamic
|
||||||
inertiaField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
|
inertiaField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
|
||||||
endif
|
endif dynamic
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! evaluate constitutive response
|
! evaluate constitutive response
|
||||||
call Utilities_constitutiveResponse(F_lastInc,x_scal,params%temperature,params%timeinc, &
|
crystallite_temperature(1,1_pInt:product(grid)) = reshape(phaseField,[product(grid)])
|
||||||
f_scal,C_volAvg,C_minmaxAvg,P_av,ForwardData,params%rotation_BC)
|
call Utilities_constitutiveResponse(F_lastInc,F,params%temperature,params%timeinc, &
|
||||||
|
residual_F,C_volAvg,C_minmaxAvg,P_av,ForwardData,params%rotation_BC)
|
||||||
ForwardData = .false.
|
ForwardData = .false.
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
|
@ -482,7 +531,7 @@ if (params%density > 0.0_pReal) then
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! updated deformation gradient using fix point algorithm of basic scheme
|
! updated deformation gradient using fix point algorithm of basic scheme
|
||||||
field_real = 0.0_pReal
|
field_real = 0.0_pReal
|
||||||
field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3) = reshape(f_scal,[grid(1),grid(2),grid(3),3,3],&
|
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],&
|
||||||
order=[4,5,1,2,3]) ! field real has a different order
|
order=[4,5,1,2,3]) ! field real has a different order
|
||||||
call Utilities_FFTforward()
|
call Utilities_FFTforward()
|
||||||
field_fourier = field_fourier + inertiaField_fourier
|
field_fourier = field_fourier + inertiaField_fourier
|
||||||
|
@ -492,7 +541,29 @@ if (params%density > 0.0_pReal) then
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! constructing residual
|
! constructing residual
|
||||||
f_scal = reshape(field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3),shape(x_scal),order=[3,4,5,1,2])
|
residual_F = reshape(field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3),&
|
||||||
|
[9,grid(1),grid(2),grid(3)],order=[2,3,4,1])
|
||||||
|
|
||||||
|
!--------------------------------------------------------------------------------------------------
|
||||||
|
! constructing phase field residual
|
||||||
|
diffusivity = 400.0; mobility = 4e6 !sample coefficients for copper... need to clear this up with better mode of input
|
||||||
|
phaseField_real = 0.0_pReal
|
||||||
|
phaseField_real(1:grid(1),1:grid(2),1:grid(3)) = phaseField_lastInc
|
||||||
|
call utilities_scalarFFTforward()
|
||||||
|
call utilities_diffusion(diffusivity/mobility,params%timeinc)
|
||||||
|
call utilities_scalarFFTbackward()
|
||||||
|
residual_phaseField = phaseField_real(1:grid(1),1:grid(2),1:grid(3))
|
||||||
|
phaseField_real = 0.0_pReal
|
||||||
|
phaseField_real(1:grid(1),1:grid(2),1:grid(3)) = &
|
||||||
|
(heatSource_lastInc + reshape(materialpoint_heat(1,1_pInt:product(grid)),[grid(1),grid(2),grid(3)]))* &
|
||||||
|
params%timeinc/mobility/2.0_pReal
|
||||||
|
call utilities_scalarFFTforward()
|
||||||
|
call utilities_diffusion(diffusivity/mobility,params%timeinc/2.0_pReal)
|
||||||
|
call utilities_scalarFFTbackward()
|
||||||
|
residual_phaseField = phaseField - &
|
||||||
|
(residual_phaseField + phaseField_real(1:grid(1),1:grid(2),1:grid(3)))
|
||||||
|
|
||||||
|
err_phaseField = maxval(abs(residual_phaseField)); phaseField_Avg = sum(phaseField)*wgt
|
||||||
|
|
||||||
end subroutine BasicPETSc_formResidual
|
end subroutine BasicPETSc_formResidual
|
||||||
|
|
||||||
|
@ -531,7 +602,8 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du
|
||||||
|
|
||||||
converged: if ((totalIter >= itmin .and. &
|
converged: if ((totalIter >= itmin .and. &
|
||||||
all([ err_div/divTol, &
|
all([ err_div/divTol, &
|
||||||
err_stress/stressTol ] < 1.0_pReal)) &
|
err_stress/stressTol, &
|
||||||
|
err_phaseField/phaseField_Avg/1.0e-3 ] < 1.0_pReal)) &
|
||||||
.or. terminallyIll) then
|
.or. terminallyIll) then
|
||||||
reason = 1
|
reason = 1
|
||||||
elseif (totalIter >= itmax) then converged
|
elseif (totalIter >= itmax) then converged
|
||||||
|
@ -547,6 +619,8 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du
|
||||||
err_div/divTol, ' (',err_div,' / m, tol =',divTol,')'
|
err_div/divTol, ' (',err_div,' / m, tol =',divTol,')'
|
||||||
write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', &
|
write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', &
|
||||||
err_stress/stressTol, ' (',err_stress, ' Pa, tol =',stressTol,')'
|
err_stress/stressTol, ' (',err_stress, ' Pa, tol =',stressTol,')'
|
||||||
|
write(6,'(a,f10.2,a,es8.2,a,es9.2,a)') ' error phase field = ', &
|
||||||
|
err_phaseField/phaseField_Avg/1.0e-3, ' (',err_phaseField/phaseField_Avg, ' Pa, tol =',1.0e-3,')'
|
||||||
write(6,'(/,a)') ' ==========================================================================='
|
write(6,'(/,a)') ' ==========================================================================='
|
||||||
flush(6)
|
flush(6)
|
||||||
|
|
||||||
|
|
|
@ -47,6 +47,8 @@ module DAMASK_spectral_utilities
|
||||||
integer(pInt), public :: grid1Red !< grid(1)/2
|
integer(pInt), public :: grid1Red !< grid(1)/2
|
||||||
real(pReal), public, dimension(:,:,:,:,:), pointer :: field_real !< real representation (some stress or deformation) of field_fourier
|
real(pReal), public, dimension(:,:,:,:,:), pointer :: field_real !< real representation (some stress or deformation) of field_fourier
|
||||||
complex(pReal),public, dimension(:,:,:,:,:), pointer :: field_fourier !< field on which the Fourier transform operates
|
complex(pReal),public, dimension(:,:,:,:,:), pointer :: field_fourier !< field on which the Fourier transform operates
|
||||||
|
real(pReal), public, dimension(:,:,:), pointer :: phaseField_real !< real representation (some stress or deformation) of field_fourier
|
||||||
|
complex(pReal),public, dimension(:,:,:), pointer :: phaseField_fourier !< field on which the Fourier transform operates
|
||||||
real(pReal), private, dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat !< gamma operator (field) for spectral method
|
real(pReal), private, dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat !< gamma operator (field) for spectral method
|
||||||
real(pReal), private, dimension(:,:,:,:), allocatable :: xi !< wave vector field for divergence and for gamma operator
|
real(pReal), private, dimension(:,:,:,:), allocatable :: xi !< wave vector field for divergence and for gamma operator
|
||||||
real(pReal), private, dimension(3,3,3,3) :: C_ref !< reference stiffness
|
real(pReal), private, dimension(3,3,3,3) :: C_ref !< reference stiffness
|
||||||
|
@ -67,6 +69,8 @@ module DAMASK_spectral_utilities
|
||||||
type(C_PTR), private :: &
|
type(C_PTR), private :: &
|
||||||
planForth, & !< FFTW plan P(x) to P(k)
|
planForth, & !< FFTW plan P(x) to P(k)
|
||||||
planBack, & !< FFTW plan F(k) to F(x)
|
planBack, & !< FFTW plan F(k) to F(x)
|
||||||
|
planPhaseFieldForth, & !< FFTW plan P(x) to P(k)
|
||||||
|
planPhaseFieldBack, & !< FFTW plan F(k) to F(x)
|
||||||
planDebugForth, & !< FFTW plan for scalar field (proof that order of usual transform is correct)
|
planDebugForth, & !< FFTW plan for scalar field (proof that order of usual transform is correct)
|
||||||
planDebugBack, & !< FFTW plan for scalar field inverse (proof that order of usual transform is correct)
|
planDebugBack, & !< FFTW plan for scalar field inverse (proof that order of usual transform is correct)
|
||||||
planDiv !< plan for FFTW in case of debugging divergence calculation
|
planDiv !< plan for FFTW in case of debugging divergence calculation
|
||||||
|
@ -101,8 +105,11 @@ module DAMASK_spectral_utilities
|
||||||
utilities_updateGamma, &
|
utilities_updateGamma, &
|
||||||
utilities_FFTforward, &
|
utilities_FFTforward, &
|
||||||
utilities_FFTbackward, &
|
utilities_FFTbackward, &
|
||||||
|
utilities_scalarFFTforward, &
|
||||||
|
utilities_scalarFFTbackward, &
|
||||||
utilities_fourierConvolution, &
|
utilities_fourierConvolution, &
|
||||||
utilities_inverseLaplace, &
|
utilities_inverseLaplace, &
|
||||||
|
utilities_diffusion, &
|
||||||
utilities_divergenceRMS, &
|
utilities_divergenceRMS, &
|
||||||
utilities_curlRMS, &
|
utilities_curlRMS, &
|
||||||
utilities_maskedCompliance, &
|
utilities_maskedCompliance, &
|
||||||
|
@ -169,6 +176,7 @@ subroutine utilities_init()
|
||||||
integer(pInt), parameter :: fileUnit = 228_pInt
|
integer(pInt), parameter :: fileUnit = 228_pInt
|
||||||
integer(pInt), dimension(3) :: k_s
|
integer(pInt), dimension(3) :: k_s
|
||||||
type(C_PTR) :: &
|
type(C_PTR) :: &
|
||||||
|
phaseFieldFFT, &
|
||||||
tensorField, & !< field cotaining data for FFTW in real and fourier space (in place)
|
tensorField, & !< field cotaining data for FFTW in real and fourier space (in place)
|
||||||
scalarField_realC, & !< field cotaining data for FFTW in real space when debugging FFTW (no in place)
|
scalarField_realC, & !< field cotaining data for FFTW in real space when debugging FFTW (no in place)
|
||||||
scalarField_fourierC, & !< field cotaining data for FFTW in fourier space when debugging FFTW (no in place)
|
scalarField_fourierC, & !< field cotaining data for FFTW in fourier space when debugging FFTW (no in place)
|
||||||
|
@ -231,8 +239,11 @@ subroutine utilities_init()
|
||||||
! allocation
|
! allocation
|
||||||
allocate (xi(3,grid1Red,grid(2),grid(3)),source = 0.0_pReal) ! frequencies, only half the size for first dimension
|
allocate (xi(3,grid1Red,grid(2),grid(3)),source = 0.0_pReal) ! frequencies, only half the size for first dimension
|
||||||
tensorField = fftw_alloc_complex(int(grid1Red*grid(2)*grid(3)*9_pInt,C_SIZE_T)) ! allocate aligned data using a C function, C_SIZE_T is of type integer(8)
|
tensorField = fftw_alloc_complex(int(grid1Red*grid(2)*grid(3)*9_pInt,C_SIZE_T)) ! allocate aligned data using a C function, C_SIZE_T is of type integer(8)
|
||||||
|
phaseFieldFFT = fftw_alloc_complex(int(grid1Red*grid(2)*grid(3),C_SIZE_T)) ! allocate aligned data using a C function, C_SIZE_T is of type integer(8)
|
||||||
call c_f_pointer(tensorField, field_real, [grid(1)+2_pInt-mod(grid(1),2_pInt),grid(2),grid(3),3,3])! place a pointer for a real representation on tensorField
|
call c_f_pointer(tensorField, field_real, [grid(1)+2_pInt-mod(grid(1),2_pInt),grid(2),grid(3),3,3])! place a pointer for a real representation on tensorField
|
||||||
call c_f_pointer(tensorField, field_fourier,[grid1Red, grid(2),grid(3),3,3])! place a pointer for a complex representation on tensorField
|
call c_f_pointer(tensorField, field_fourier,[grid1Red, grid(2),grid(3),3,3])! place a pointer for a complex representation on tensorField
|
||||||
|
call c_f_pointer(phaseFieldFFT,phaseField_real,[grid(1)+2_pInt-mod(grid(1),2_pInt),grid(2),grid(3)])! place a pointer for a real representation on tensorField
|
||||||
|
call c_f_pointer(phaseFieldFFT,phaseField_fourier,[grid1Red,grid(2),grid(3)])! place a pointer for a complex representation on tensorField
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! general initialization of FFTW (see manual on fftw.org for more details)
|
! general initialization of FFTW (see manual on fftw.org for more details)
|
||||||
|
@ -259,6 +270,21 @@ subroutine utilities_init()
|
||||||
1, grid(3)*grid(2)*(grid(1)+2_pInt-mod(grid(1),2_pInt)), & ! striding, product of physical length in the 3 dimensions
|
1, grid(3)*grid(2)*(grid(1)+2_pInt-mod(grid(1),2_pInt)), & ! striding, product of physical length in the 3 dimensions
|
||||||
fftw_planner_flag) ! planner precision
|
fftw_planner_flag) ! planner precision
|
||||||
|
|
||||||
|
!--------------------------------------------------------------------------------------------------
|
||||||
|
! creating plans for the convolution
|
||||||
|
planPhaseFieldForth = fftw_plan_many_dft_r2c(3,[grid(3),grid(2) ,grid(1)], 1, & ! dimensions, logical length in each dimension in reversed order, no. of transforms
|
||||||
|
phaseField_real,[grid(3),grid(2) ,grid(1)+2_pInt-mod(grid(1),2_pInt)], & ! input data, physical length in each dimension in reversed order
|
||||||
|
1, grid(3)*grid(2)*(grid(1)+2_pInt-mod(grid(1),2_pInt)), & ! striding, product of physical length in the 3 dimensions
|
||||||
|
phaseField_fourier,[grid(3),grid(2) ,grid1Red], & ! output data, physical length in each dimension in reversed order
|
||||||
|
1, grid(3)*grid(2)* grid1Red, fftw_planner_flag) ! striding, product of physical length in the 3 dimensions, planner precision
|
||||||
|
|
||||||
|
planPhaseFieldBack = fftw_plan_many_dft_c2r(3,[grid(3),grid(2) ,grid(1)], 1, & ! dimensions, logical length in each dimension in reversed order, no. of transforms
|
||||||
|
phaseField_fourier,[grid(3),grid(2) ,grid1Red], & ! input data, physical length in each dimension in reversed order
|
||||||
|
1, grid(3)*grid(2)* grid1Red, & ! striding, product of physical length in the 3 dimensions
|
||||||
|
phaseField_real,[grid(3),grid(2) ,grid(1)+2_pInt-mod(grid(1),2_pInt)], & ! output data, physical length in each dimension in reversed order
|
||||||
|
1, grid(3)*grid(2)*(grid(1)+2_pInt-mod(grid(1),2_pInt)), & ! striding, product of physical length in the 3 dimensions
|
||||||
|
fftw_planner_flag) ! planner precision
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! depending on debug options, allocate more memory and create additional plans
|
! depending on debug options, allocate more memory and create additional plans
|
||||||
if (debugDivergence) then
|
if (debugDivergence) then
|
||||||
|
@ -487,6 +513,59 @@ subroutine utilities_FFTbackward()
|
||||||
end subroutine utilities_FFTbackward
|
end subroutine utilities_FFTbackward
|
||||||
|
|
||||||
|
|
||||||
|
!--------------------------------------------------------------------------------------------------
|
||||||
|
!> @brief forward FFT of data in field_real to field_fourier with highest freqs. removed
|
||||||
|
!> @details Does an unweighted FFT transform from real to complex.
|
||||||
|
!> In case of debugging the FFT, also one component of the tensor (specified by row and column)
|
||||||
|
!> is independetly transformed complex to complex and compared to the whole tensor transform
|
||||||
|
!--------------------------------------------------------------------------------------------------
|
||||||
|
subroutine utilities_scalarFFTforward()
|
||||||
|
use math
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
integer(pInt), dimension(2:3,2) :: Nyquist ! highest frequencies to be removed (1 if even, 2 if odd)
|
||||||
|
|
||||||
|
!--------------------------------------------------------------------------------------------------
|
||||||
|
! doing the FFT
|
||||||
|
call fftw_execute_dft_r2c(planPhaseFieldForth,phaseField_real,phaseField_fourier)
|
||||||
|
|
||||||
|
!--------------------------------------------------------------------------------------------------
|
||||||
|
! removing highest frequencies
|
||||||
|
Nyquist(2,1:2) = [grid(2)/2_pInt + 1_pInt, grid(2)/2_pInt + 1_pInt + mod(grid(2),2_pInt)]
|
||||||
|
Nyquist(3,1:2) = [grid(3)/2_pInt + 1_pInt, grid(3)/2_pInt + 1_pInt + mod(grid(3),2_pInt)]
|
||||||
|
|
||||||
|
if(grid(1)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation
|
||||||
|
phaseField_fourier (grid1Red, 1:grid(2), 1:grid(3)) &
|
||||||
|
= cmplx(0.0_pReal,0.0_pReal,pReal)
|
||||||
|
if(grid(2)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation
|
||||||
|
phaseField_fourier (1:grid1Red,Nyquist(2,1):Nyquist(2,2),1:grid(3)) &
|
||||||
|
= cmplx(0.0_pReal,0.0_pReal,pReal)
|
||||||
|
if(grid(3)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation
|
||||||
|
phaseField_fourier (1:grid1Red,1:grid(2),Nyquist(3,1):Nyquist(3,2)) &
|
||||||
|
= cmplx(0.0_pReal,0.0_pReal,pReal)
|
||||||
|
end subroutine utilities_scalarFFTforward
|
||||||
|
|
||||||
|
|
||||||
|
!--------------------------------------------------------------------------------------------------
|
||||||
|
!> @brief backward FFT of data in field_fourier to field_real
|
||||||
|
!> @details Does an inverse FFT transform from complex to real
|
||||||
|
!> In case of debugging the FFT, also one component of the tensor (specified by row and column)
|
||||||
|
!> is independetly transformed complex to complex and compared to the whole tensor transform
|
||||||
|
!> results is weighted by number of points stored in wgt
|
||||||
|
!--------------------------------------------------------------------------------------------------
|
||||||
|
subroutine utilities_scalarFFTbackward()
|
||||||
|
use math !< must use the whole module for use of FFTW
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
!--------------------------------------------------------------------------------------------------
|
||||||
|
! doing the iFFT
|
||||||
|
call fftw_execute_dft_c2r(planPhaseFieldBack,phaseField_fourier,phaseField_real) ! back transform of fluct deformation gradient
|
||||||
|
phaseField_real = phaseField_real * wgt ! normalize the result by number of elements
|
||||||
|
|
||||||
|
end subroutine utilities_scalarFFTbackward
|
||||||
|
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
!> @brief doing convolution with inverse laplace kernel
|
!> @brief doing convolution with inverse laplace kernel
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
|
@ -522,6 +601,39 @@ field_fourier(1,1,1,1:3,1:3) = cmplx(0.0_pReal,0.0_pReal,pReal)
|
||||||
end subroutine utilities_inverseLaplace
|
end subroutine utilities_inverseLaplace
|
||||||
|
|
||||||
|
|
||||||
|
!--------------------------------------------------------------------------------------------------
|
||||||
|
!> @brief doing convolution with inverse laplace kernel
|
||||||
|
!--------------------------------------------------------------------------------------------------
|
||||||
|
subroutine utilities_diffusion(coefficient,timeinc)
|
||||||
|
use math, only: &
|
||||||
|
PI
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
real(pReal),intent(in) :: timeinc, coefficient
|
||||||
|
integer(pInt) :: i, j, k
|
||||||
|
integer(pInt), dimension(3) :: k_s
|
||||||
|
|
||||||
|
write(6,'(/,a)') ' ... doing diffusion .......................................................'
|
||||||
|
flush(6)
|
||||||
|
|
||||||
|
!--------------------------------------------------------------------------------------------------
|
||||||
|
! do the actual spectral method calculation (mechanical equilibrium)
|
||||||
|
do k = 1_pInt, grid(3)
|
||||||
|
k_s(3) = k - 1_pInt
|
||||||
|
if(k > grid(3)/2_pInt + 1_pInt) k_s(3) = k_s(3) - grid(3) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1
|
||||||
|
do j = 1_pInt, grid(2)
|
||||||
|
k_s(2) = j - 1_pInt
|
||||||
|
if(j > grid(2)/2_pInt + 1_pInt) k_s(2) = k_s(2) - grid(2) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1
|
||||||
|
do i = 1_pInt, grid1Red
|
||||||
|
k_s(1) = i - 1_pInt
|
||||||
|
phaseField_fourier(i,j,k) = phaseField_fourier(i,j,k)* &
|
||||||
|
cmplx(exp(-sum((2.0_pReal*PI*real(k_s,pReal)/geomSize)*(2.0_pReal*PI*real(k_s,pReal)/geomSize))* &
|
||||||
|
coefficient*timeinc),0.0_pReal,pReal) ! symmetry, junst running from 0,1,...,N/2,N/2+1
|
||||||
|
enddo; enddo; enddo
|
||||||
|
|
||||||
|
end subroutine utilities_diffusion
|
||||||
|
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
!> @brief doing convolution gamma_hat * field_real, ensuring that average value = fieldAim
|
!> @brief doing convolution gamma_hat * field_real, ensuring that average value = fieldAim
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
|
@ -1007,19 +1119,20 @@ end function utilities_getFilter
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
!> @brief calculates filter for fourier convolution depending on type given in numerics.config
|
!> @brief calculates filter for fourier convolution depending on type given in numerics.config
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
subroutine utilities_temperatureUpdate()
|
subroutine utilities_temperatureUpdate(timeinc)
|
||||||
use crystallite, only: &
|
use crystallite, only: &
|
||||||
crystallite_temperature
|
crystallite_temperature
|
||||||
use homogenization, only: &
|
use homogenization, only: &
|
||||||
materialpoint_heat
|
materialpoint_heat
|
||||||
|
|
||||||
implicit none
|
implicit none
|
||||||
|
real(pReal),intent(in) :: timeinc
|
||||||
integer :: &
|
integer :: &
|
||||||
x,y,z,e
|
x,y,z,e
|
||||||
real(pReal) :: &
|
real(pReal) :: &
|
||||||
a
|
a
|
||||||
forall(e=1_pInt:product(grid)) &
|
forall(e=1_pInt:product(grid)) &
|
||||||
crystallite_temperature(1,e) = crystallite_temperature(1,e) + materialpoint_heat(1,e)*0.0001_pReal
|
crystallite_temperature(1,e) = crystallite_temperature(1,e) + materialpoint_heat(1,e)*timeinc
|
||||||
e = 0_pInt
|
e = 0_pInt
|
||||||
z = 0_pInt
|
z = 0_pInt
|
||||||
y = 0_pInt
|
y = 0_pInt
|
||||||
|
|
|
@ -1303,9 +1303,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity)
|
||||||
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
|
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
|
||||||
do g = 1,myNgrains
|
do g = 1,myNgrains
|
||||||
crystallite_heat(g,i,e) = 0.98_pReal* &
|
crystallite_heat(g,i,e) = 0.98_pReal* &
|
||||||
math_mul33xx33(crystallite_P(1:3,1:3,g,i,e), &
|
abs(math_mul33xx33(math_Mandel6to33(crystallite_Tstar_v), &
|
||||||
crystallite_partionedF(1:3,1:3,g,i,e)-&
|
crystallite_Lp(1:3,1:3,g,i,e)))
|
||||||
crystallite_partionedF0(1:3,1:3,g,i,e))
|
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo elementLooping12
|
enddo elementLooping12
|
||||||
|
|
Loading…
Reference in New Issue