rolled back phase field changes

This commit is contained in:
Pratheek Shanthraj 2013-12-18 09:09:32 +00:00
parent 6fa9ed8f48
commit ff6211b78c
3 changed files with 48 additions and 337 deletions

View File

@ -76,9 +76,7 @@ program DAMASK_spectral_Driver
geomSize, &
tBoundaryCondition, &
tSolutionState, &
phaseFieldDataBin, &
cutBack, &
maxPhaseFields
cutBack
use DAMASK_spectral_SolverBasic
#ifdef PETSc
use DAMASK_spectral_SolverBasicPETSC
@ -91,17 +89,14 @@ program DAMASK_spectral_Driver
real(pReal), dimension (3,3) :: rotation = math_I3 !< rotation of BC
type(tBoundaryCondition) :: P, & !< stress BC
deformation !< deformation BC (Fdot or L)
type(phaseFieldDataBin) :: phaseFieldData(maxPhaseFields)
real(pReal) :: time = 0.0_pReal, & !< length of increment
temperature = 300.0_pReal, & !< isothermal starting condition
temperature = 300.0_pReal, & !< isothermal starting conditions
density = 0.0_pReal !< density
integer(pInt) :: incs = 0_pInt, & !< number of increments
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
thermal_active = .false., & !< activate thermal phase field
fracture_active = .false. !< activate fracture phase field
logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase
end type tLoadCase
!--------------------------------------------------------------------------------------------------
@ -117,8 +112,7 @@ program DAMASK_spectral_Driver
integer(pInt) :: &
N_t = 0_pInt, & !< # of time indicators found in load case file
N_n = 0_pInt, & !< # of increment specifiers found in load case file
N_def = 0_pInt, & !< # of rate of deformation specifiers found in load case file
nActivePhaseFields
N_def = 0_pInt !< # of rate of deformation specifiers found in load case file
character(len=65536) :: &
line
@ -236,22 +230,8 @@ program DAMASK_spectral_Driver
loadCases(currentLoadCase)%P%values = math_plain9to33(temp_valueVector)
case('t','time','delta') ! increment time
loadCases(currentLoadCase)%time = IO_floatValue(line,positions,i+1_pInt)
case('temperature')
case('temp','temperature') ! starting temperature
loadCases(currentLoadCase)%temperature = IO_floatValue(line,positions,i+1_pInt)
case('thermal') ! starting temperature, conductivity and mobility
loadCases(:)%phaseFieldData(1)%label = 'thermal'
loadCases(:)%phaseFieldData(1)%active = .true.
loadCases(:)%phaseFieldData(1)%phaseField0 = 300.0_pReal ! initialize to meaningful value if not defined
loadCases(currentLoadCase)%phaseFieldData(1)%phaseField0 = IO_floatValue(line,positions,i+1_pInt)
loadCases(currentLoadCase)%phaseFieldData(1)%diffusion = IO_floatValue(line,positions,i+2_pInt)
loadCases(currentLoadCase)%phaseFieldData(1)%mobility = IO_floatValue(line,positions,i+3_pInt)
case('fracture') ! starting damage, diffusion and mobility
loadCases(:)%phaseFieldData(2)%label = 'fracture'
loadCases(:)%phaseFieldData(2)%active = .true.
loadCases(:)%phaseFieldData(2)%phaseField0 = 1.0_pReal ! initialize to meaningful value if not defined
loadCases(currentLoadCase)%phaseFieldData(2)%phaseField0 = IO_floatValue(line,positions,i+1_pInt)
loadCases(currentLoadCase)%phaseFieldData(2)%diffusion = IO_floatValue(line,positions,i+2_pInt)
loadCases(currentLoadCase)%phaseFieldData(2)%mobility = IO_floatValue(line,positions,i+3_pInt)
case('den','density') ! starting density
loadCases(currentLoadCase)%density = IO_floatValue(line,positions,i+1_pInt)
case('n','incs','increments','steps') ! number of increments
@ -291,14 +271,6 @@ program DAMASK_spectral_Driver
end select
enddo; enddo
close(FILEUNIT)
! reorder phase field data to remove redundant non-active fields
nActivePhaseFields = 0_pInt
do i = 1, maxPhaseFields
if (loadCases(1)%phaseFieldData(i)%active) then
nActivePhaseFields = nActivePhaseFields + 1_pInt
loadCases(:)%phaseFieldData(nActivePhaseFields) = loadCases(:)%phaseFieldData(i)
endif
enddo
!--------------------------------------------------------------------------------------------------
! consistency checks and output of load case
@ -364,7 +336,7 @@ program DAMASK_spectral_Driver
call basic_init(loadCases(1)%temperature)
#ifdef PETSc
case (DAMASK_spectral_SolverBasicPETSc_label)
call basicPETSc_init(loadCases(1)%temperature,nActivePhaseFields,loadCases(1)%phaseFieldData(1:nActivePhaseFields))
call basicPETSc_init(loadCases(1)%temperature)
case (DAMASK_spectral_SolverAL_label)
if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) &
call IO_warning(42_pInt, ext_msg='debug Divergence')
@ -495,9 +467,7 @@ program DAMASK_spectral_Driver
F_BC = loadCases(currentLoadCase)%deformation, &
temperature_bc = loadCases(currentLoadCase)%temperature, &
rotation_BC = loadCases(currentLoadCase)%rotation, &
density = loadCases(currentLoadCase)%density, &
nActivePhaseFields = nActivePhaseFields, &
phaseFieldData = loadCases(1)%phaseFieldData(1:nActivePhaseFields))
density = loadCases(currentLoadCase)%density)
case (DAMASK_spectral_SolverAL_label)
solres = AL_solution (&
incInfo,guess,timeinc,timeIncOld,remainingLoadCaseTime, &

View File

@ -31,9 +31,7 @@ module DAMASK_spectral_SolverBasicPETSc
use math, only: &
math_I3
use DAMASK_spectral_Utilities, only: &
tSolutionState, &
phaseFieldDataBin, &
maxPhaseFields
tSolutionState
implicit none
private
@ -51,8 +49,6 @@ module DAMASK_spectral_SolverBasicPETSc
real(pReal) :: timeincOld
real(pReal) :: temperature
real(pReal) :: density
integer(pInt) :: nActivePhaseFields
type(phaseFieldDataBin) :: phaseFieldData(maxPhaseFields)
end type tSolutionParams
type(tSolutionParams), private :: params
@ -66,11 +62,6 @@ module DAMASK_spectral_SolverBasicPETSc
!--------------------------------------------------------------------------------------------------
! common pointwise data
real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, Fdot, F_lastInc2
real(pReal), private, dimension(:,:,:,:), allocatable :: &
phaseFieldRHS_lastInc, &
phaseField_lastInc, &
phaseFieldRHS, &
phaseFieldDot
complex(pReal), private, dimension(:,:,:,:,:), allocatable :: inertiaField_fourier
!--------------------------------------------------------------------------------------------------
@ -88,7 +79,6 @@ module DAMASK_spectral_SolverBasicPETSc
C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 stiffness
S = 0.0_pReal !< current compliance (filled up with zeros)
real(pReal), private :: err_stress, err_div, err_divPrev, err_divDummy
real(pReal), private, dimension(:), allocatable :: err_phaseField, phaseField_Avg
logical, private :: ForwardData
integer(pInt), private :: &
totalIter = 0_pInt !< total iteration in current increment
@ -121,7 +111,7 @@ contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields and fills them with data, potentially from restart info
!--------------------------------------------------------------------------------------------------
subroutine basicPETSc_init(temperature,nActivePhaseFields,phaseFieldData)
subroutine basicPETSc_init(temperature)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment)
use IO, only: &
IO_intOut, &
@ -150,14 +140,13 @@ subroutine basicPETSc_init(temperature,nActivePhaseFields,phaseFieldData)
math_invSym3333
implicit none
integer(pInt), intent(in) :: nActivePhaseFields
type(phaseFieldDataBin), intent(in) :: phaseFieldData(nActivePhaseFields)
real(pReal), intent(inOut) :: temperature
real(pReal), intent(inout) :: &
temperature
#include <finclude/petscdmda.h90>
#include <finclude/petscsnes.h90>
#include <finclude/petscvec.h>
real(pReal), dimension(:,:,:,:,:), allocatable :: P
PetscScalar, dimension(:,:,:,:), pointer :: xx_psc, F
PetscScalar, dimension(:,:,:,:), pointer :: F
PetscErrorCode :: ierr
PetscObject :: dummy
real(pReal), dimension(3,3) :: &
@ -165,7 +154,6 @@ subroutine basicPETSc_init(temperature,nActivePhaseFields,phaseFieldData)
real(pReal), dimension(3,3,3,3) :: &
temp3333_Real = 0.0_pReal
KSP :: ksp
integer(pInt) :: i
call Utilities_init()
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasicPETSc init -+>>>'
@ -173,19 +161,13 @@ subroutine basicPETSc_init(temperature,nActivePhaseFields,phaseFieldData)
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
allocate (P (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! allocate global fields
allocate (P (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
allocate (F_lastInc (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 (inertiaField_fourier (grid1Red,grid(2),grid(3),3,3),source = cmplx(0.0_pReal,0.0_pReal,pReal))
allocate (phaseFieldRHS_lastInc (nActivePhaseFields,grid(1),grid(2),grid(3)),source = 0.0_pReal)
allocate (phaseField_lastInc (nActivePhaseFields,grid(1),grid(2),grid(3)),source = 0.0_pReal)
allocate (phaseFieldDot (nActivePhaseFields,grid(1),grid(2),grid(3)),source = 0.0_pReal)
allocate (phaseFieldRHS (nActivePhaseFields,grid(1),grid(2),grid(3)),source = 0.0_pReal)
allocate (err_phaseField(nActivePhaseFields), source = 0.0_pReal)
allocate (phaseField_Avg(nActivePhaseFields), source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
@ -193,7 +175,7 @@ subroutine basicPETSc_init(temperature,nActivePhaseFields,phaseFieldData)
call DMDACreate3d(PETSC_COMM_WORLD, &
DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, &
DMDA_STENCIL_BOX,grid(1),grid(2),grid(3),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, &
9+nActivePhaseFields,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)
9,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)
CHKERRQ(ierr)
call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr)
call DMDASNESSetFunctionLocal(da,INSERT_VALUES,BasicPETSC_formResidual,dummy,ierr)
@ -208,16 +190,12 @@ subroutine basicPETSc_init(temperature,nActivePhaseFields,phaseFieldData)
!--------------------------------------------------------------------------------------------------
! init fields
call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with
F => xx_psc(0:8,:,:,:)
if (restartInc == 1_pInt) then ! no deformation (no restart)
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with
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
xx_psc(0:8,:,:,:) = 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
do i = 1, nActivePhaseFields
xx_psc(8+i,:,:,:) = phaseFieldData(i)%phaseField0
phaseField_lastInc(i,:,:,:) = phaseFieldData(i)%phaseField0
enddo
elseif (restartInc > 1_pInt) then ! using old values from file
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) &
write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
@ -251,10 +229,10 @@ subroutine basicPETSc_init(temperature,nActivePhaseFields,phaseFieldData)
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(&
F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
call Utilities_constitutiveResponse(&
reshape(F,[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(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)]),&
temperature,0.0_pReal,P,C_volAvg,C_minmaxAvg,temp33_Real,.false.,math_I3)
call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! write data back into PETSc
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! write data back into PETSc
if (restartInc == 1_pInt) then ! use initial stiffness as reference stiffness
temp3333_Real = C_minMaxAvg
endif
@ -267,8 +245,7 @@ end subroutine basicPETSc_init
!> @brief solution for the Basic PETSC scheme with internal iterations
!--------------------------------------------------------------------------------------------------
type(tSolutionState) function basicPETSc_solution( &
incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,temperature_bc,rotation_BC,density, &
nActivePhaseFields,phaseFieldData)
incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,temperature_bc,rotation_BC,density)
use numerics, only: &
update_gamma, &
itmax
@ -292,8 +269,6 @@ type(tSolutionState) function basicPETSc_solution( &
use FEsolving, only: &
restartWrite, &
terminallyIll
use homogenization, only: &
materialpoint_heat
implicit none
#include <finclude/petscdmda.h90>
@ -301,8 +276,6 @@ type(tSolutionState) function basicPETSc_solution( &
!--------------------------------------------------------------------------------------------------
! input data for solution
integer(pInt), intent(in) :: nActivePhaseFields
type(phaseFieldDataBin), intent(in) :: phaseFieldData(nActivePhaseFields)
real(pReal), intent(in) :: &
timeinc, & !< increment in time for current solution
timeinc_old, & !< increment in time of last increment
@ -317,17 +290,15 @@ type(tSolutionState) function basicPETSc_solution( &
character(len=*), intent(in) :: &
incInfoIn
real(pReal), dimension(3,3), intent(in) :: rotation_BC
integer(pInt) :: i
!--------------------------------------------------------------------------------------------------
! PETSc Data
PetscScalar, pointer :: xx_psc(:,:,:,:), F(:,:,:,:)
PetscScalar, pointer :: F(:,:,:,:)
PetscErrorCode :: ierr
SNESConvergedReason :: reason
incInfo = incInfoIn
call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with
F => xx_psc(0:8,:,:,:)
call DMDAVecGetArrayF90(da,solution_vec,F,ierr)
!--------------------------------------------------------------------------------------------------
! restart information for spectral solver
if (restartWrite) then
@ -356,11 +327,8 @@ type(tSolutionState) function basicPETSc_solution( &
F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
if (cutBack) then
F_aim = F_aim_lastInc
xx_psc(0:8,:,:,:) = 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
do i = 1, nActivePhaseFields
xx_psc(8+i,:,:,:) = phaseField_lastInc(i,:,:,:)
enddo
else
C_volAvgLastInc = C_volAvg
@ -383,11 +351,6 @@ type(tSolutionState) function basicPETSc_solution( &
F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
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)]))
do i = 1, nActivePhaseFields
phaseFieldDot(i,:,:,:) = (xx_psc(8+i,:,:,:) - phaseField_lastInc(i,:,:,:))/timeinc_old
phaseField_lastInc(i,:,:,:) = xx_psc(8+i,:,:,:)
phaseFieldRHS_lastInc(i,:,:,:) = phaseFieldRHS(i,:,:,:)
enddo
F_lastInc2 = F_lastInc
F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid(3)])
endif
@ -395,10 +358,7 @@ type(tSolutionState) function basicPETSc_solution( &
F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot,math_rotate_backward33(F_aim, &
rotation_BC)),[9,grid(1),grid(2),grid(3)])
do i = 1, nActivePhaseFields
xx_psc(8+i,:,:,:) = phaseField_lastInc(i,:,:,:) + phaseFieldDot(i,:,:,:)*timeinc
enddo
call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr)
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
@ -416,8 +376,6 @@ type(tSolutionState) function basicPETSc_solution( &
params%timeincOld = timeinc_old
params%temperature = temperature_BC
params%density = density
params%nActivePhaseFields = nActivePhaseFields
params%phaseFieldData(1:nActivePhaseFields) = phaseFieldData(1:nActivePhaseFields)
call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr); CHKERRQ(ierr)
call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr)
@ -456,46 +414,26 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
wgt, &
field_real, &
field_fourier, &
phaseField_real, &
phaseField_fourier, &
Utilities_FFTforward, &
Utilities_FFTbackward, &
utilities_scalarFFTforward, &
utilities_scalarFFTbackward, &
Utilities_fourierConvolution, &
Utilities_inverseLaplace, &
Utilities_diffusion, &
Utilities_constitutiveResponse, &
Utilities_divergenceRMS
use IO, only: &
IO_intOut
use crystallite, only: &
crystallite_temperature
use homogenization, only: &
materialpoint_heat, &
materialpoint_P
implicit none
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &
in
PetscScalar, target, dimension(9+params%nActivePhaseFields, &
XG_RANGE,YG_RANGE,ZG_RANGE) :: &
x_scal
PetscScalar, target, dimension(9+params%nActivePhaseFields, &
X_RANGE,Y_RANGE,Z_RANGE) :: &
PetscScalar, dimension(3,3,grid(1),grid(2),grid(3)) :: &
x_scal, &
f_scal
PetscScalar, pointer, dimension(:,:,:,:) :: &
F, &
residual_F
PetscInt :: &
PETScIter, &
nfuncs
PetscObject :: dummy
PetscErrorCode :: ierr
integer(pInt) :: i
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))
call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr)
call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr)
@ -517,38 +455,26 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! evaluate inertia
dynamic: if (params%density > 0.0_pReal) then
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)
residual_F = params%density*product(geomSize/grid)*residual_F
if (params%density > 0.0_pReal) then
f_scal = ((x_scal - F_lastInc)/params%timeinc - (F_lastInc - F_lastInc2)/params%timeincOld)/&
((params%timeinc + params%timeincOld)/2.0_pReal)
f_scal = params%density*product(geomSize/grid)*f_scal
field_real = 0.0_pReal
field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3) = reshape(residual_F,[grid(1),grid(2),grid(3),3,3],&
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],&
order=[4,5,1,2,3]) ! field real has a different order
call Utilities_FFTforward()
call Utilities_inverseLaplace()
inertiaField_fourier = field_fourier
else dynamic
else
inertiaField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
endif dynamic
endif
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
do i = 1, params%nActivePhaseFields
if(params%phaseFieldData(i)%label == 'thermal') &
crystallite_temperature(1,1_pInt:product(grid)) = &
reshape(x_scal(9+i,1:grid(1),1:grid(2),1:grid(3)),[product(grid)])
enddo
call Utilities_constitutiveResponse(F_lastInc,F,params%temperature,params%timeinc, &
residual_F,C_volAvg,C_minmaxAvg,P_av,ForwardData,params%rotation_BC)
call Utilities_constitutiveResponse(F_lastInc,x_scal,params%temperature,params%timeinc, &
f_scal,C_volAvg,C_minmaxAvg,P_av,ForwardData,params%rotation_BC)
ForwardData = .false.
do i = 1, params%nActivePhaseFields
if(params%phaseFieldData(i)%label == 'fracture') &
residual_F = residual_F * spread(x_scal(9+i,1:grid(1),1:grid(2),1:grid(3)),dim=1,ncopies=9)
enddo
!--------------------------------------------------------------------------------------------------
! stress BC handling
F_aim_lastIter = F_aim
@ -558,7 +484,7 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! updated deformation gradient using fix point algorithm of basic scheme
field_real = 0.0_pReal
field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3) = reshape(residual_F,[grid(1),grid(2),grid(3),3,3],&
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],&
order=[4,5,1,2,3]) ! field real has a different order
call Utilities_FFTforward()
field_fourier = field_fourier + inertiaField_fourier
@ -566,78 +492,9 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
call Utilities_fourierConvolution(math_rotate_backward33(F_aim_lastIter-F_aim,params%rotation_BC))
call Utilities_FFTbackward()
!--------------------------------------------------------------------------------------------------
! constructing phase field residual
do i = 1, params%nActivePhaseFields
select case (params%phaseFieldData(i)%label)
case ('thermal')
phaseField_real = 0.0_pReal
phaseField_real(1:grid(1),1:grid(2),1:grid(3)) = &
phaseField_lastInc(i,1:grid(1),1:grid(2),1:grid(3))
call utilities_scalarFFTforward()
call utilities_diffusion(params%phaseFieldData(i)%diffusion,params%timeinc)
call utilities_scalarFFTbackward()
f_scal(9+i,1:grid(1),1:grid(2),1:grid(3)) = &
phaseField_real(1:grid(1),1:grid(2),1:grid(3))
phaseFieldRHS(i,1:grid(1),1:grid(2),1:grid(3)) = &
reshape(materialpoint_heat(1,1_pInt:product(grid)),[grid(1),grid(2),grid(3)])
phaseField_real = 0.0_pReal
phaseField_real(1:grid(1),1:grid(2),1:grid(3)) = &
params%timeinc*params%phaseFieldData(i)%mobility* &
(phaseFieldRHS_lastInc(i,1:grid(1),1:grid(2),1:grid(3)) + &
phaseFieldRHS (i,1:grid(1),1:grid(2),1:grid(3)))/2.0_pReal
call utilities_scalarFFTforward()
call utilities_diffusion(params%phaseFieldData(i)%diffusion,params%timeinc/2.0_pReal)
call utilities_scalarFFTbackward()
f_scal(9+i,1:grid(1),1:grid(2),1:grid(3)) = &
x_scal(9+i,1:grid(1),1:grid(2),1:grid(3)) - &
f_scal(9+i,1:grid(1),1:grid(2),1:grid(3)) - &
phaseField_real(1:grid(1),1:grid(2),1:grid(3))
err_phaseField(i) = maxval(abs(f_scal(9+i,1:grid(1),1:grid(2),1:grid(3))))
phaseField_Avg(i) = sum(x_scal(9+i,1:grid(1),1:grid(2),1:grid(3)))*wgt
case ('fracture')
phaseField_real = 0.0_pReal
phaseField_real(1:grid(1),1:grid(2),1:grid(3)) = &
phaseField_lastInc(i,1:grid(1),1:grid(2),1:grid(3))
call utilities_scalarFFTforward()
call utilities_diffusion(2.0_pReal*maxval(geomSize/real(grid,pReal))* &
params%phaseFieldData(i)%diffusion,params%timeinc)
call utilities_scalarFFTbackward()
f_scal(9+i,1:grid(1),1:grid(2),1:grid(3)) = &
phaseField_real(1:grid(1),1:grid(2),1:grid(3))
phaseFieldRHS(i,1:grid(1),1:grid(2),1:grid(3)) = &
- params%phaseFieldData(i)%mobility* &
sum(residual_F* &
(F-reshape(spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid(3)),[9,grid(1),grid(2),grid(3)])),dim=1) &
- params%phaseFieldData(i)%diffusion*(x_scal(9+i,1:grid(1),1:grid(2),1:grid(3)) - 1.0_pReal)/ &
8.0_pReal/maxval(geomSize/real(grid,pReal))
phaseField_real = 0.0_pReal
phaseField_real(1:grid(1),1:grid(2),1:grid(3)) = &
params%timeinc*params%phaseFieldData(i)%mobility* &
(phaseFieldRHS_lastInc(i,1:grid(1),1:grid(2),1:grid(3)) + &
phaseFieldRHS (i,1:grid(1),1:grid(2),1:grid(3)))/2.0_pReal
call utilities_scalarFFTforward()
call utilities_diffusion(2.0_pReal*maxval(geomSize/real(grid,pReal))* &
params%phaseFieldData(i)%diffusion,params%timeinc/2.0_pReal)
call utilities_scalarFFTbackward()
f_scal(9+i,1:grid(1),1:grid(2),1:grid(3)) = &
x_scal(9+i,1:grid(1),1:grid(2),1:grid(3)) - &
f_scal(9+i,1:grid(1),1:grid(2),1:grid(3)) - &
phaseField_real(1:grid(1),1:grid(2),1:grid(3))
err_phaseField(i) = maxval(abs(f_scal(9+i,1:grid(1),1:grid(2),1:grid(3))))
phaseField_Avg(i) = sum(x_scal(9+i,1:grid(1),1:grid(2),1:grid(3)))*wgt
end select
enddo
!--------------------------------------------------------------------------------------------------
! constructing residual
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])
f_scal = reshape(field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3),shape(f_scal),order=[3,4,5,1,2])
end subroutine BasicPETSc_formResidual
@ -668,17 +525,15 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du
PetscErrorCode :: ierr
real(pReal) :: &
divTol, &
stressTol, &
phaseField_err = 0.0_pReal
stressTol
divTol = max(maxval(abs(P_av))*err_div_tolRel,err_div_tolAbs)
stressTol = max(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs)
err_divPrev = err_div; err_div = err_divDummy
if (params%nActivePhaseFields .ne. 0_pInt) phaseField_err = maxval(err_phaseField/phaseField_Avg)
converged: if ((totalIter >= itmin .and. &
all([ err_div/divTol, err_stress/stressTol] < 1.0_pReal) .and. &
phaseField_err < 1.0e-3_pReal) &
all([ err_div/divTol, &
err_stress/stressTol ] < 1.0_pReal)) &
.or. terminallyIll) then
reason = 1
elseif (totalIter >= itmax) then converged
@ -694,16 +549,11 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du
err_div/divTol, ' (',err_div,' / m, tol =',divTol,')'
write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', &
err_stress/stressTol, ' (',err_stress, ' Pa, tol =',stressTol,')'
if (params%nActivePhaseFields .ne. 0_pInt) then
write(6,'(a,f10.2,a,es8.2,a,es9.2,a)') ' error phase field = ', &
phaseField_err/1.0e-3, ' (',phaseField_err, ' Pa, tol =',1.0e-3,')'
endif
write(6,'(/,a)') ' ==========================================================================='
flush(6)
end subroutine BasicPETSc_converged
!--------------------------------------------------------------------------------------------------
!> @brief convergence check
!--------------------------------------------------------------------------------------------------

View File

@ -48,8 +48,6 @@ module DAMASK_spectral_utilities
integer(pInt), public :: grid1Red !< grid(1)/2
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
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 :: xi !< wave vector field for divergence and for gamma operator
real(pReal), private, dimension(3,3,3,3) :: C_ref !< reference stiffness
@ -70,8 +68,6 @@ module DAMASK_spectral_utilities
type(C_PTR), private :: &
planForth, & !< FFTW plan P(x) to P(k)
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)
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
@ -114,11 +110,8 @@ module DAMASK_spectral_utilities
utilities_updateGamma, &
utilities_FFTforward, &
utilities_FFTbackward, &
utilities_scalarFFTforward, &
utilities_scalarFFTbackward, &
utilities_fourierConvolution, &
utilities_inverseLaplace, &
utilities_diffusion, &
utilities_divergenceRMS, &
utilities_curlRMS, &
utilities_maskedCompliance, &
@ -184,7 +177,6 @@ subroutine utilities_init()
integer(pInt), parameter :: fileUnit = 228_pInt
integer(pInt), dimension(3) :: k_s
type(C_PTR) :: &
phaseFieldFFT, &
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_fourierC, & !< field cotaining data for FFTW in fourier space when debugging FFTW (no in place)
@ -247,11 +239,8 @@ subroutine utilities_init()
! allocation
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)
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_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)
@ -278,21 +267,6 @@ 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
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
if (debugDivergence) then
@ -308,7 +282,7 @@ subroutine utilities_init()
endif
if (debugFFTW) then
scalarField_realC = fftw_alloc_complex(int(product(grid),C_SIZE_T)) ! allllocate data for real representation (no in place transform)
scalarField_realC = fftw_alloc_complex(int(product(grid),C_SIZE_T)) ! allocate data for real representation (no in place transform)
scalarField_fourierC = fftw_alloc_complex(int(product(grid),C_SIZE_T)) ! allocate data for fourier representation (no in place transform)
call c_f_pointer(scalarField_realC, scalarField_real, grid) ! place a pointer for a real representation
call c_f_pointer(scalarField_fourierC, scalarField_fourier, grid) ! place a pointer for a fourier representation
@ -521,59 +495,6 @@ 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
!--------------------------------------------------------------------------------------------------
@ -609,36 +530,6 @@ field_fourier(1,1,1,1:3,1:3) = cmplx(0.0_pReal,0.0_pReal,pReal)
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
!--------------------------------------------------------------------------------------------------
! 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
!--------------------------------------------------------------------------------------------------