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

reworked phase field interface to damask spectral solvers. now specify 'thermal a b c' or 'fracture a b c' to activate either phase field where a b c are the initial value, diffusion coefficient and mobility respectively. 

Right now only thermal and fracture phase fields implemented and only in the basic petsc solver
This commit is contained in:
Pratheek Shanthraj 2013-11-13 19:21:35 +00:00
parent d2cbca35b2
commit 50db944c0c
3 changed files with 189 additions and 145 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

43
code/DAMASK_spectral_driver.f90
View File

@ -75,8 +75,9 @@ program DAMASK_spectral_Driver
geomSize, & geomSize, &
tBoundaryCondition, & tBoundaryCondition, &
tSolutionState, & tSolutionState, &
phaseFieldDataBin, &
cutBack, & cutBack, &
utilities_temperatureUpdate maxPhaseFields
use DAMASK_spectral_SolverBasic use DAMASK_spectral_SolverBasic
#ifdef PETSc #ifdef PETSc
use DAMASK_spectral_SolverBasicPETSC use DAMASK_spectral_SolverBasicPETSC
@ -89,14 +90,17 @@ program DAMASK_spectral_Driver
real(pReal), dimension (3,3) :: rotation = math_I3 !< rotation of BC real(pReal), dimension (3,3) :: rotation = math_I3 !< rotation of BC
type(tBoundaryCondition) :: P, & !< stress BC type(tBoundaryCondition) :: P, & !< stress BC
deformation !< deformation BC (Fdot or L) deformation !< deformation BC (Fdot or L)
type(phaseFieldDataBin) :: phaseFieldData(maxPhaseFields)
real(pReal) :: time = 0.0_pReal, & !< length of increment real(pReal) :: time = 0.0_pReal, & !< length of increment
temperature = 300.0_pReal, & !< isothermal starting conditions temperature = 300.0_pReal, & !< isothermal starting condition
density = 0.0_pReal !< density density = 0.0_pReal !< density
integer(pInt) :: incs = 0_pInt, & !< number of increments integer(pInt) :: incs = 0_pInt, & !< number of increments
outputfrequency = 1_pInt, & !< frequency of result writes outputfrequency = 1_pInt, & !< frequency of result writes
restartfrequency = 0_pInt, & !< frequency of restart writes restartfrequency = 0_pInt, & !< frequency of restart writes
logscale = 0_pInt !< linear/logarithmic time inc flag logscale = 0_pInt !< linear/logarithmic time inc flag
logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase logical :: followFormerTrajectory = .true., & !< follow trajectory of former loadcase
thermal_active = .false., & !< activate thermal phase field
fracture_active = .false. !< activate fracture phase field
end type tLoadCase end type tLoadCase
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -112,7 +116,8 @@ program DAMASK_spectral_Driver
integer(pInt) :: & integer(pInt) :: &
N_t = 0_pInt, & !< # of time indicators found in load case file 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_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 N_def = 0_pInt, & !< # of rate of deformation specifiers found in load case file
nActivePhaseFields
character(len=65536) :: & character(len=65536) :: &
line line
@ -230,8 +235,22 @@ program DAMASK_spectral_Driver
loadCases(currentLoadCase)%P%values = math_plain9to33(temp_valueVector) loadCases(currentLoadCase)%P%values = math_plain9to33(temp_valueVector)
case('t','time','delta') ! increment time case('t','time','delta') ! increment time
loadCases(currentLoadCase)%time = IO_floatValue(line,positions,i+1_pInt) loadCases(currentLoadCase)%time = IO_floatValue(line,positions,i+1_pInt)
case('temp','temperature') ! starting temperature case('temperature')
loadCases(currentLoadCase)%temperature = IO_floatValue(line,positions,i+1_pInt) 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 case('den','density') ! starting density
loadCases(currentLoadCase)%density = IO_floatValue(line,positions,i+1_pInt) loadCases(currentLoadCase)%density = IO_floatValue(line,positions,i+1_pInt)
case('n','incs','increments','steps') ! number of increments case('n','incs','increments','steps') ! number of increments
@ -271,6 +290,14 @@ program DAMASK_spectral_Driver
end select end select
enddo; enddo enddo; enddo
close(myUnit) close(myUnit)
! 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 ! consistency checks and output of load case
@ -336,7 +363,7 @@ program DAMASK_spectral_Driver
call basic_init(loadCases(1)%temperature) call basic_init(loadCases(1)%temperature)
#ifdef PETSc #ifdef PETSc
case (DAMASK_spectral_SolverBasicPETSc_label) case (DAMASK_spectral_SolverBasicPETSc_label)
call basicPETSc_init(loadCases(1)%temperature) call basicPETSc_init(loadCases(1)%temperature,nActivePhaseFields,loadCases(1)%phaseFieldData(1:nActivePhaseFields))
case (DAMASK_spectral_SolverAL_label) case (DAMASK_spectral_SolverAL_label)
if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) & if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) &
call IO_warning(42_pInt, ext_msg='debug Divergence') call IO_warning(42_pInt, ext_msg='debug Divergence')
@ -467,7 +494,9 @@ program DAMASK_spectral_Driver
F_BC = loadCases(currentLoadCase)%deformation, & F_BC = loadCases(currentLoadCase)%deformation, &
temperature_bc = loadCases(currentLoadCase)%temperature, & temperature_bc = loadCases(currentLoadCase)%temperature, &
rotation_BC = loadCases(currentLoadCase)%rotation, & rotation_BC = loadCases(currentLoadCase)%rotation, &
density = loadCases(currentLoadCase)%density) density = loadCases(currentLoadCase)%density, &
nActivePhaseFields = nActivePhaseFields, &
phaseFieldData = loadCases(1)%phaseFieldData(1:nActivePhaseFields))
case (DAMASK_spectral_SolverAL_label) case (DAMASK_spectral_SolverAL_label)
solres = AL_solution (& solres = AL_solution (&
incInfo,guess,timeinc,timeIncOld,remainingLoadCaseTime, & incInfo,guess,timeinc,timeIncOld,remainingLoadCaseTime, &

210
code/DAMASK_spectral_solverBasicPETSc.f90
View File

@ -31,7 +31,9 @@ module DAMASK_spectral_SolverBasicPETSc
use math, only: & use math, only: &
math_I3 math_I3
use DAMASK_spectral_Utilities, only: & use DAMASK_spectral_Utilities, only: &
tSolutionState tSolutionState, &
phaseFieldDataBin, &
maxPhaseFields
implicit none implicit none
private private
@ -49,6 +51,8 @@ module DAMASK_spectral_SolverBasicPETSc
real(pReal) :: timeincOld real(pReal) :: timeincOld
real(pReal) :: temperature real(pReal) :: temperature
real(pReal) :: density real(pReal) :: density
integer(pInt) :: nActivePhaseFields
type(phaseFieldDataBin) :: phaseFieldData(maxPhaseFields)
end type tSolutionParams end type tSolutionParams
type(tSolutionParams), private :: params type(tSolutionParams), private :: params
@ -62,9 +66,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 :: & real(pReal), private, dimension(:,:,:,:), allocatable :: &
heatSource_lastInc, & phaseFieldRHS_lastInc, &
phaseField_lastInc, & phaseField_lastInc, &
phaseFieldRHS, &
phaseFieldDot phaseFieldDot
complex(pReal), private, dimension(:,:,:,:,:), allocatable :: inertiaField_fourier complex(pReal), private, dimension(:,:,:,:,:), allocatable :: inertiaField_fourier
@ -82,7 +87,8 @@ 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, err_phaseField, phaseField_Avg real(pReal), private :: err_stress, err_div, err_divPrev, err_divDummy
real(pReal), private, dimension(:), allocatable :: 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
@ -115,7 +121,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields and fills them with data, potentially from restart info !> @brief allocates all neccessary fields and fills them with data, potentially from restart info
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine basicPETSc_init(temperature) subroutine basicPETSc_init(temperature,nActivePhaseFields,phaseFieldData)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment)
use IO, only: & use IO, only: &
IO_intOut, & IO_intOut, &
@ -144,14 +150,14 @@ subroutine basicPETSc_init(temperature)
math_invSym3333 math_invSym3333
implicit none implicit none
real(pReal), intent(inout) :: & integer(pInt), intent(in) :: nActivePhaseFields
temperature type(phaseFieldDataBin), intent(in) :: phaseFieldData(nActivePhaseFields)
real(pReal), intent(inOut) :: temperature
#include <finclude/petscdmda.h90> #include <finclude/petscdmda.h90>
#include <finclude/petscsnes.h90> #include <finclude/petscsnes.h90>
#include <finclude/petscvec.h> #include <finclude/petscvec.h>
real(pReal), dimension(:,:,:,:,:), allocatable :: P real(pReal), dimension(:,:,:,:,:), allocatable :: P
PetscScalar, dimension(:,:,:,:), pointer :: xx_psc, 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) :: &
@ -159,6 +165,7 @@ subroutine basicPETSc_init(temperature)
real(pReal), dimension(3,3,3,3) :: & real(pReal), dimension(3,3,3,3) :: &
temp3333_Real = 0.0_pReal temp3333_Real = 0.0_pReal
KSP :: ksp KSP :: ksp
integer(pInt) :: i
call Utilities_init() call Utilities_init()
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasicPETSc init -+>>>' write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasicPETSc init -+>>>'
@ -166,16 +173,19 @@ subroutine basicPETSc_init(temperature)
write(6,'(a15,a)') ' Current time: ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
allocate (P (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate global fields ! 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_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 (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 (phaseFieldRHS_lastInc (nActivePhaseFields,grid(1),grid(2),grid(3)),source = 0.0_pReal)
allocate (phaseField_lastInc(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 (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 ! initialize solver specific parts of PETSc
@ -183,31 +193,29 @@ 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, &
10,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr) 9+nActivePhaseFields,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 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); CHKERRQ(ierr)
CHKERRQ(ierr)
call SNESGetKSP(snes,ksp,ierr); CHKERRQ(ierr) call SNESGetKSP(snes,ksp,ierr); CHKERRQ(ierr)
call KSPSetConvergenceTest(ksp,BasicPETSC_convergedKSP,dummy,PETSC_NULL_FUNCTION,ierr) call KSPSetConvergenceTest(ksp,BasicPETSC_convergedKSP,dummy,PETSC_NULL_FUNCTION,ierr); CHKERRQ(ierr)
CHKERRQ(ierr)
call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr) call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init fields ! init fields
call DMDAVecGetArrayF90(da,solution_vec,xx_psc,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,:,:,:) 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)]) xx_psc(0:8,:,:,:) = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)])
F_lastInc2 = F_lastInc F_lastInc2 = F_lastInc
phaseField = temperature do i = 1, nActivePhaseFields
phaseField_lastInc = temperature xx_psc(8+i,:,:,:) = phaseFieldData(i)%phaseField0
phaseField_lastInc(i,:,:,:) = phaseFieldData(i)%phaseField0
enddo
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)') &
@ -257,7 +265,8 @@ end subroutine basicPETSc_init
!> @brief solution for the Basic PETSC scheme with internal iterations !> @brief solution for the Basic PETSC scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
type(tSolutionState) function basicPETSc_solution( & type(tSolutionState) function basicPETSc_solution( &
incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,temperature_bc,rotation_BC,density) incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,temperature_bc,rotation_BC,density, &
nActivePhaseFields,phaseFieldData)
use numerics, only: & use numerics, only: &
update_gamma, & update_gamma, &
itmax itmax
@ -290,6 +299,8 @@ type(tSolutionState) function basicPETSc_solution( &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! input data for solution ! input data for solution
integer(pInt), intent(in) :: nActivePhaseFields
type(phaseFieldDataBin), intent(in) :: phaseFieldData(nActivePhaseFields)
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
timeinc, & !< increment in time for current solution timeinc, & !< increment in time for current solution
timeinc_old, & !< increment in time of last increment timeinc_old, & !< increment in time of last increment
@ -304,17 +315,17 @@ type(tSolutionState) function basicPETSc_solution( &
character(len=*), intent(in) :: & character(len=*), intent(in) :: &
incInfoIn incInfoIn
real(pReal), dimension(3,3), intent(in) :: rotation_BC real(pReal), dimension(3,3), intent(in) :: rotation_BC
integer(pInt) :: i
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc Data ! PETSc Data
PetscScalar, pointer :: xx_psc(:,:,:,:), F(:,:,:,:), phaseField(:,:,:) PetscScalar, pointer :: xx_psc(:,:,:,:), F(:,:,:,:)
PetscErrorCode :: ierr PetscErrorCode :: ierr
SNESConvergedReason :: reason SNESConvergedReason :: reason
incInfo = incInfoIn incInfo = incInfoIn
call DMDAVecGetArrayF90(da,solution_vec,xx_psc,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,:,:,:) 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
@ -343,9 +354,11 @@ 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)])
if (cutBack) then if (cutBack) then
F_aim = F_aim_lastInc F_aim = F_aim_lastInc
F = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)]) xx_psc(0:8,:,:,:) = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)])
C_volAvg = C_volAvgLastInc C_volAvg = C_volAvgLastInc
phaseField = phaseField_lastInc do i = 1, nActivePhaseFields
xx_psc(8+i,:,:,:) = phaseField_lastInc(i,:,:,:)
enddo
else else
C_volAvgLastInc = C_volAvg C_volAvgLastInc = C_volAvg
@ -368,17 +381,21 @@ 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 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_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)])
phaseField = phaseField_lastInc + phaseFieldDot*timeinc 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,xx_psc,ierr); CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -397,6 +414,8 @@ type(tSolutionState) function basicPETSc_solution( &
params%timeincOld = timeinc_old params%timeincOld = timeinc_old
params%temperature = temperature_BC params%temperature = temperature_BC
params%density = density 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 SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr); CHKERRQ(ierr)
call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr) call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr)
@ -451,34 +470,32 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
use crystallite, only: & use crystallite, only: &
crystallite_temperature crystallite_temperature
use homogenization, only: & use homogenization, only: &
materialpoint_heat materialpoint_heat, &
materialpoint_P
use constitutive, only: &
constitutive_damage
implicit none implicit none
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &
in in
PetscScalar, target, dimension(10, & PetscScalar, target, dimension(9+params%nActivePhaseFields, &
XG_RANGE,YG_RANGE,ZG_RANGE) :: & XG_RANGE,YG_RANGE,ZG_RANGE) :: &
x_scal x_scal
PetscScalar, target, dimension(10, & PetscScalar, target, dimension(9+params%nActivePhaseFields, &
X_RANGE,Y_RANGE,Z_RANGE) :: & X_RANGE,Y_RANGE,Z_RANGE) :: &
f_scal f_scal
PetscScalar, pointer, dimension(:,:,:,:) :: & PetscScalar, pointer, dimension(:,:,:,:) :: &
F, & F, &
residual_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 integer(pInt) :: i
F => x_scal(1:9,1:grid(1),1:grid(2),1:grid(3)) 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)) 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)
@ -517,7 +534,15 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate constitutive response ! evaluate constitutive response
crystallite_temperature(1,1_pInt:product(grid)) = reshape(phaseField,[product(grid)]) 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)])
if(params%phaseFieldData(i)%label == 'fracture') &
constitutive_damage(1,1,1: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, & call Utilities_constitutiveResponse(F_lastInc,F,params%temperature,params%timeinc, &
residual_F,C_volAvg,C_minmaxAvg,P_av,ForwardData,params%rotation_BC) residual_F,C_volAvg,C_minmaxAvg,P_av,ForwardData,params%rotation_BC)
ForwardData = .false. ForwardData = .false.
@ -539,32 +564,78 @@ 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_fourierConvolution(math_rotate_backward33(F_aim_lastIter-F_aim,params%rotation_BC))
call Utilities_FFTbackward() 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(0.5_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)/ &
2.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 ! constructing residual
residual_F = reshape(field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3),& 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]) [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
@ -601,9 +672,8 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du
err_divPrev = err_div; err_div = err_divDummy err_divPrev = err_div; err_div = err_divDummy
converged: if ((totalIter >= itmin .and. & converged: if ((totalIter >= itmin .and. &
all([ err_div/divTol, & all([ err_div/divTol, err_stress/stressTol] < 1.0_pReal) .and. &
err_stress/stressTol, & maxval(err_phaseField/phaseField_Avg) < 1.0e-3_pReal) &
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
@ -620,7 +690,7 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du
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 = ', & 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,')' maxval(err_phaseField/phaseField_Avg)/1.0e-3, ' (',maxval(err_phaseField/phaseField_Avg), ' Pa, tol =',1.0e-3,')'
write(6,'(/,a)') ' ===========================================================================' write(6,'(/,a)') ' ==========================================================================='
flush(6) flush(6)

75
code/DAMASK_spectral_utilities.f90
View File

@ -36,6 +36,7 @@ module DAMASK_spectral_utilities
#include <finclude/petscsys.h> #include <finclude/petscsys.h>
#endif #endif
logical, public :: cutBack =.false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill 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
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! grid related information information ! grid related information information
integer(pInt), public, dimension(3) :: grid !< grid points as specified in geometry file integer(pInt), public, dimension(3) :: grid !< grid points as specified in geometry file
@ -100,6 +101,14 @@ module DAMASK_spectral_utilities
character(len=64) :: myType = 'None' character(len=64) :: myType = 'None'
end type tBoundaryCondition end type tBoundaryCondition
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
logical :: active = .false.
character(len=64) :: label = ''
end type phaseFieldDataBin
public :: & public :: &
utilities_init, & utilities_init, &
utilities_updateGamma, & utilities_updateGamma, &
@ -116,8 +125,7 @@ module DAMASK_spectral_utilities
utilities_constitutiveResponse, & utilities_constitutiveResponse, &
utilities_calculateRate, & utilities_calculateRate, &
utilities_forwardField, & utilities_forwardField, &
utilities_destroy, & utilities_destroy
utilities_temperatureUpdate
private :: & private :: &
utilities_getFilter utilities_getFilter
@ -613,9 +621,6 @@ subroutine utilities_diffusion(coefficient,timeinc)
integer(pInt) :: i, j, k integer(pInt) :: i, j, k
integer(pInt), dimension(3) :: k_s integer(pInt), dimension(3) :: k_s
write(6,'(/,a)') ' ... doing diffusion .......................................................'
flush(6)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! do the actual spectral method calculation (mechanical equilibrium) ! do the actual spectral method calculation (mechanical equilibrium)
do k = 1_pInt, grid(3) do k = 1_pInt, grid(3)
@ -1116,66 +1121,6 @@ real(pReal) function utilities_getFilter(k)
end function utilities_getFilter end function utilities_getFilter
!--------------------------------------------------------------------------------------------------
!> @brief calculates filter for fourier convolution depending on type given in numerics.config
!--------------------------------------------------------------------------------------------------
subroutine utilities_temperatureUpdate(timeinc)
use crystallite, only: &
crystallite_temperature
use homogenization, only: &
materialpoint_heat
implicit none
real(pReal),intent(in) :: timeinc
integer :: &
x,y,z,e
real(pReal) :: &
a
forall(e=1_pInt:product(grid)) &
crystallite_temperature(1,e) = crystallite_temperature(1,e) + materialpoint_heat(1,e)*timeinc
e = 0_pInt
z = 0_pInt
y = 0_pInt
x = 0_pInt
!< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me]
do z = 0_pInt,grid(3)-1_pInt
do y = 0_pInt,grid(2)-1_pInt
do x = 0_pInt,grid(1)-1_pInt
e = e + 1_pInt
a = 0.0_pReal
a = a+ crystallite_temperature(1, z * grid(1) * grid(2) &
+ y * grid(1) &
+ modulo(x+1_pInt,grid(1)) &
+ 1_pInt)
a = a+ crystallite_temperature(1,z * grid(1) * grid(2) &
+ y * grid(1) &
+ modulo(x-1_pInt,grid(1)) &
+ 1_pInt)
a = a+ crystallite_temperature(1,z * grid(1) * grid(2) &
+ modulo(y+1_pInt,grid(2)) * grid(1) &
+ x &
+ 1_pInt)
a = a+ crystallite_temperature(1,z * grid(1) * grid(2) &
+ modulo(y-1_pInt,grid(2)) * grid(1) &
+ x &
+ 1_pInt)
a = a+ crystallite_temperature(1, modulo(z+1_pInt,grid(3)) * grid(1) * grid(2) &
+ y * grid(1) &
+ x &
+ 1_pInt)
a = a+ crystallite_temperature(1,modulo(z-1_pInt,grid(3)) * grid(1) * grid(2) &
+ y * grid(1) &
+ x &
+ 1_pInt)
crystallite_temperature(1,e) = (crystallite_temperature(1,e)+a)/7.0_pReal
enddo
enddo
enddo
end subroutine utilities_temperatureUpdate
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief cleans up !> @brief cleans up
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------