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

introduced dynamics for high strain rate simulations. so far only works for materials with homogeneous density. the density is read in as an optional argument in the load file using the keyword 'den' or 'density'. if density > 0 then inertial terms are added to the stress equilibrium. an implicit euler time discretisation is used to calculate the inertial terms. 

dynamic problems are length scale dependent so make sure the geometry size is meaningful (1 m^3 or N m^3 for example is too big for a RVE!).

not specifying a density (density = 0) will perform a quasi static simulation as before.
This commit is contained in:
Pratheek Shanthraj 2013-07-26 16:25:37 +00:00
parent 187a70569e
commit 03d8f14a98
5 changed files with 160 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
code/DAMASK_spectral_driver.f90
View File

@ -91,7 +91,8 @@ program DAMASK_spectral_Driver
type(tBoundaryCondition) :: P, & !< stress BC type(tBoundaryCondition) :: P, & !< stress BC
deformation !< deformation BC (Fdot or L) deformation !< deformation BC (Fdot or L)
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 conditions
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
@ -232,6 +233,8 @@ program DAMASK_spectral_Driver
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('temp','temperature') ! starting temperature
loadCases(currentLoadCase)%temperature = IO_floatValue(line,positions,i+1_pInt) loadCases(currentLoadCase)%temperature = IO_floatValue(line,positions,i+1_pInt)
case('den','density') ! starting density
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
loadCases(currentLoadCase)%incs = IO_intValue(line,positions,i+1_pInt) loadCases(currentLoadCase)%incs = IO_intValue(line,positions,i+1_pInt)
case('logincs','logincrements','logsteps') ! number of increments (switch to log time scaling) case('logincs','logincrements','logsteps') ! number of increments (switch to log time scaling)
@ -314,6 +317,7 @@ program DAMASK_spectral_Driver
write(6,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'rotation of loadframe:',& write(6,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'rotation of loadframe:',&
math_transpose33(loadCases(currentLoadCase)%rotation) math_transpose33(loadCases(currentLoadCase)%rotation)
write(6,'(2x,a,f12.6)') 'temperature:', loadCases(currentLoadCase)%temperature write(6,'(2x,a,f12.6)') 'temperature:', loadCases(currentLoadCase)%temperature
write(6,'(2x,a,f12.6)') 'density: ', loadCases(currentLoadCase)%density
if (loadCases(currentLoadCase)%time < 0.0_pReal) errorID = 834_pInt ! negative time increment if (loadCases(currentLoadCase)%time < 0.0_pReal) errorID = 834_pInt ! negative time increment
write(6,'(2x,a,f12.6)') 'time: ', loadCases(currentLoadCase)%time write(6,'(2x,a,f12.6)') 'time: ', loadCases(currentLoadCase)%time
if (loadCases(currentLoadCase)%incs < 1_pInt) errorID = 835_pInt ! non-positive incs count if (loadCases(currentLoadCase)%incs < 1_pInt) errorID = 835_pInt ! non-positive incs count
@ -463,21 +467,24 @@ program DAMASK_spectral_Driver
P_BC = loadCases(currentLoadCase)%P, & P_BC = loadCases(currentLoadCase)%P, &
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)
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, &
P_BC = loadCases(currentLoadCase)%P, & P_BC = loadCases(currentLoadCase)%P, &
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)
case (DAMASK_spectral_SolverPolarisation_label) case (DAMASK_spectral_SolverPolarisation_label)
solres = Polarisation_solution (& solres = Polarisation_solution (&
incInfo,guess,timeinc,timeIncOld,remainingLoadCaseTime, & incInfo,guess,timeinc,timeIncOld,remainingLoadCaseTime, &
P_BC = loadCases(currentLoadCase)%P, & P_BC = loadCases(currentLoadCase)%P, &
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)
#endif #endif
end select end select

45
code/DAMASK_spectral_solverAL.f90
View File

@ -47,7 +47,9 @@ module DAMASK_spectral_solverAL
type tSolutionParams !< @todo use here the type definition for a full loadcase including mask type tSolutionParams !< @todo use here the type definition for a full loadcase including mask
real(pReal), dimension(3,3) :: P_BC, rotation_BC real(pReal), dimension(3,3) :: P_BC, rotation_BC
real(pReal) :: timeinc real(pReal) :: timeinc
real(pReal) :: timeincOld
real(pReal) :: temperature real(pReal) :: temperature
real(pReal) :: density
end type tSolutionParams end type tSolutionParams
type(tSolutionParams), private :: params type(tSolutionParams), private :: params
@ -63,9 +65,11 @@ module DAMASK_spectral_solverAL
! common pointwise data ! common pointwise data
real(pReal), private, dimension(:,:,:,:,:), allocatable :: & real(pReal), private, dimension(:,:,:,:,:), allocatable :: &
F_lastInc, & !< field of previous compatible deformation gradients F_lastInc, & !< field of previous compatible deformation gradients
F_lambda_lastInc, & !< field of previous incompatible deformation gradient F_lastInc2, & !< field of 2nd previous compatible deformation gradients
F_lambda_lastInc, & !< field of previous incompatible deformation gradient
Fdot, & !< field of assumed rate of compatible deformation gradient Fdot, & !< field of assumed rate of compatible deformation gradient
F_lambdaDot !< field of assumed rate of incopatible deformation gradient F_lambdaDot !< field of assumed rate of incopatible deformation gradient
complex(pReal),private, dimension(:,:,:,:,:), allocatable :: inertiaField_fourier
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc. ! stress, stiffness and compliance average etc.
@ -140,6 +144,7 @@ subroutine AL_init(temperature)
Utilities_constitutiveResponse, & Utilities_constitutiveResponse, &
Utilities_updateGamma, & Utilities_updateGamma, &
grid, & grid, &
grid1Red, &
geomSize, & geomSize, &
wgt wgt
use mesh, only: & use mesh, only: &
@ -174,9 +179,11 @@ subroutine AL_init(temperature)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate global fields ! allocate global fields
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 (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 (F_lambda_lastInc(3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal) allocate (F_lambda_lastInc(3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
allocate (F_lambdaDot (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal) allocate (F_lambdaDot (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))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc Init ! PETSc Init
@ -200,6 +207,7 @@ subroutine AL_init(temperature)
F_lambda => xx_psc(9:17,:,:,:) F_lambda => xx_psc(9:17,:,:,:)
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_lastInc2 = F_lastInc
F_lambda_lastInc = F_lastInc F_lambda_lastInc = F_lastInc
F = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)]) F = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)])
F_lambda = F F_lambda = F
@ -216,6 +224,10 @@ subroutine AL_init(temperature)
trim(getSolverJobName()),size(F_lastInc)) trim(getSolverJobName()),size(F_lastInc))
read (777,rec=1) F_lastInc read (777,rec=1) F_lastInc
close (777) close (777)
call IO_read_jobBinaryFile(777,'F_lastInc2',&
trim(getSolverJobName()),size(F_lastInc2))
read (777,rec=1) F_lastInc2
close (777)
F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F
F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc
call IO_read_jobBinaryFile(777,'F_lambda',& call IO_read_jobBinaryFile(777,'F_lambda',&
@ -263,7 +275,7 @@ end subroutine AL_init
!> @brief solution for the AL scheme with internal iterations !> @brief solution for the AL scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
type(tSolutionState) function & type(tSolutionState) function &
AL_solution(incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,temperature_bc,rotation_BC) AL_solution(incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,temperature_bc,rotation_BC,density)
use numerics, only: & use numerics, only: &
update_gamma, & update_gamma, &
itmax, & itmax, &
@ -311,7 +323,7 @@ use mesh, only: &
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
real(pReal), intent(in) :: density
real(pReal) :: err_stress_tol real(pReal) :: err_stress_tol
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc Data ! PETSc Data
@ -380,7 +392,7 @@ use mesh, only: &
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)]))
F_lambdaDot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & F_lambdaDot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
timeinc_old,guess,F_lambda_lastInc,reshape(F_lambda,[3,3,grid(1),grid(2),grid(3)])) timeinc_old,guess,F_lambda_lastInc,reshape(F_lambda,[3,3,grid(1),grid(2),grid(3)]))
F_lastInc2 = F_lastInc
F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid(3)]) F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid(3)])
F_lambda_lastInc = reshape(F_lambda,[3,3,grid(1),grid(2),grid(3)]) F_lambda_lastInc = reshape(F_lambda,[3,3,grid(1),grid(2),grid(3)])
endif endif
@ -408,7 +420,9 @@ use mesh, only: &
params%P_BC = P_BC%values params%P_BC = P_BC%values
params%rotation_BC = rotation_BC params%rotation_BC = rotation_BC
params%timeinc = timeinc params%timeinc = timeinc
params%timeincOld = timeinc_old
params%temperature = temperature_bc params%temperature = temperature_bc
params%density = density
err_stress_tol = 1.0_pReal; err_stress = 2.0_pReal * err_stress_tol ! ensure to start loop err_stress_tol = 1.0_pReal; err_stress = 2.0_pReal * err_stress_tol ! ensure to start loop
totalIter = -1_pInt totalIter = -1_pInt
@ -461,10 +475,13 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
math_invSym3333 math_invSym3333
use DAMASK_spectral_Utilities, only: & use DAMASK_spectral_Utilities, only: &
grid, & grid, &
geomSize, &
wgt, & wgt, &
field_real, & field_real, &
field_fourier, &
Utilities_FFTforward, & Utilities_FFTforward, &
Utilities_fourierConvolution, & Utilities_fourierConvolution, &
Utilities_inverseLaplace, &
Utilities_FFTbackward, & Utilities_FFTbackward, &
Utilities_constitutiveResponse Utilities_constitutiveResponse
use debug, only: & use debug, only: &
@ -528,7 +545,22 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
math_transpose33(F_aim) math_transpose33(F_aim)
flush(6) flush(6)
endif endif
if (params%density > 0.0_pReal) then
!--------------------------------------------------------------------------------------------------
! evaluate inertia
residual_F = ((F - F_lastInc)/params%timeinc - (F_lastInc - F_lastInc2)/params%timeincOld)/((params%timeinc + params%timeincOld)/2.0_pReal)
residual_F = params%density*product(geomSize/grid)*residual_F
field_real = 0.0_pReal
field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3) = reshape(residual_F,[grid(1),grid(2),grid(3),3,3],&
order=[4,5,1,2,3]) ! field real has a different order
call Utilities_FFTforward()
call Utilities_inverseLaplace()
inertiaField_fourier = field_fourier
else
inertiaField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! !
field_real = 0.0_pReal field_real = 0.0_pReal
@ -540,6 +572,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! doing convolution in Fourier space ! doing convolution in Fourier space
call Utilities_FFTforward() call Utilities_FFTforward()
field_fourier = field_fourier + polarAlpha*inertiaField_fourier
call Utilities_fourierConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) call Utilities_fourierConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC))
call Utilities_FFTbackward() call Utilities_FFTbackward()

42
code/DAMASK_spectral_solverBasicPETSc.f90
View File

@ -46,7 +46,9 @@ module DAMASK_spectral_SolverBasicPETSc
type tSolutionParams type tSolutionParams
real(pReal), dimension(3,3) :: P_BC, rotation_BC real(pReal), dimension(3,3) :: P_BC, rotation_BC
real(pReal) :: timeinc real(pReal) :: timeinc
real(pReal) :: timeincOld
real(pReal) :: temperature real(pReal) :: temperature
real(pReal) :: density
end type tSolutionParams end type tSolutionParams
type(tSolutionParams), private :: params type(tSolutionParams), private :: params
@ -59,7 +61,8 @@ module DAMASK_spectral_SolverBasicPETSc
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! common pointwise data ! common pointwise data
real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, Fdot real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, Fdot, F_lastInc2
complex(pReal), private, dimension(:,:,:,:,:), allocatable :: inertiaField_fourier
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc. ! stress, stiffness and compliance average etc.
@ -126,6 +129,7 @@ subroutine basicPETSc_init(temperature)
Utilities_constitutiveResponse, & Utilities_constitutiveResponse, &
Utilities_updateGamma, & Utilities_updateGamma, &
grid, & grid, &
grid1Red, &
wgt, & wgt, &
geomSize geomSize
use mesh, only: & use mesh, only: &
@ -157,8 +161,10 @@ subroutine basicPETSc_init(temperature)
allocate (P (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal) allocate (P (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate global fields ! allocate global fields
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 (Fdot (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))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc ! initialize solver specific parts of PETSc
@ -182,6 +188,7 @@ subroutine basicPETSc_init(temperature)
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
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)') &
@ -195,6 +202,10 @@ subroutine basicPETSc_init(temperature)
trim(getSolverJobName()),size(F_lastInc)) trim(getSolverJobName()),size(F_lastInc))
read (777,rec=1) F_lastInc read (777,rec=1) F_lastInc
close (777) close (777)
call IO_read_jobBinaryFile(777,'F_lastInc2',&
trim(getSolverJobName()),size(F_lastInc2))
read (777,rec=1) F_lastInc2
close (777)
F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F
F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc
@ -230,7 +241,7 @@ 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) incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,temperature_bc,rotation_BC,density)
use numerics, only: & use numerics, only: &
update_gamma, & update_gamma, &
itmax itmax
@ -269,6 +280,7 @@ type(tSolutionState) function basicPETSc_solution( &
type(tBoundaryCondition), intent(in) :: P_BC,F_BC type(tBoundaryCondition), intent(in) :: P_BC,F_BC
real(pReal), dimension(3,3), intent(in) :: rotation_BC real(pReal), dimension(3,3), intent(in) :: rotation_BC
character(len=*), intent(in) :: incInfoIn character(len=*), intent(in) :: incInfoIn
real(pReal), intent(in) :: density
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! !
@ -326,6 +338,7 @@ type(tSolutionState) function basicPETSc_solution( &
! update rate and forward last inc ! update rate and forward last inc
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)]))
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)])
endif endif
F_aim = F_aim + f_aimDot * timeinc F_aim = F_aim + f_aimDot * timeinc
@ -347,7 +360,9 @@ type(tSolutionState) function basicPETSc_solution( &
params%P_BC = P_BC%values params%P_BC = P_BC%values
params%rotation_BC = rotation_BC params%rotation_BC = rotation_BC
params%timeinc = timeinc params%timeinc = timeinc
params%timeincOld = timeinc_old
params%temperature = temperature_BC params%temperature = temperature_BC
params%density = density
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)
@ -382,11 +397,14 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
debug_spectralRotation debug_spectralRotation
use DAMASK_spectral_Utilities, only: & use DAMASK_spectral_Utilities, only: &
grid, & grid, &
geomSize, &
wgt, & wgt, &
field_real, & field_real, &
field_fourier, &
Utilities_FFTforward, & Utilities_FFTforward, &
Utilities_FFTbackward, & Utilities_FFTbackward, &
Utilities_fourierConvolution, & Utilities_fourierConvolution, &
Utilities_inverseLaplace, &
Utilities_constitutiveResponse, & Utilities_constitutiveResponse, &
Utilities_divergenceRMS Utilities_divergenceRMS
use IO, only: & use IO, only: &
@ -424,6 +442,21 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
flush(6) flush(6)
endif endif
if (params%density > 0.0_pReal) then
!--------------------------------------------------------------------------------------------------
! evaluate inertia
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(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
inertiaField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate constitutive response ! evaluate constitutive response
call Utilities_constitutiveResponse(F_lastInc,x_scal,params%temperature,params%timeinc, & call Utilities_constitutiveResponse(F_lastInc,x_scal,params%temperature,params%timeinc, &
@ -442,6 +475,7 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
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(f_scal,[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
err_div = Utilities_divergenceRMS() err_div = Utilities_divergenceRMS()
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()

38
code/DAMASK_spectral_solverPolarisation.f90
View File

@ -47,7 +47,9 @@ module DAMASK_spectral_solverPolarisation
type tSolutionParams !< @todo use here the type definition for a full loadcase including mask type tSolutionParams !< @todo use here the type definition for a full loadcase including mask
real(pReal), dimension(3,3) :: P_BC, rotation_BC real(pReal), dimension(3,3) :: P_BC, rotation_BC
real(pReal) :: timeinc real(pReal) :: timeinc
real(pReal) :: timeincOld
real(pReal) :: temperature real(pReal) :: temperature
real(pReal) :: density
end type tSolutionParams end type tSolutionParams
type(tSolutionParams), private :: params type(tSolutionParams), private :: params
@ -63,9 +65,11 @@ module DAMASK_spectral_solverPolarisation
! common pointwise data ! common pointwise data
real(pReal), private, dimension(:,:,:,:,:), allocatable :: & real(pReal), private, dimension(:,:,:,:,:), allocatable :: &
F_lastInc, & !< field of previous compatible deformation gradients F_lastInc, & !< field of previous compatible deformation gradients
F_lastInc2, & !< field of 2nd previous compatible deformation gradients
F_tau_lastInc, & !< field of previous incompatible deformation gradient F_tau_lastInc, & !< field of previous incompatible deformation gradient
Fdot, & !< field of assumed rate of compatible deformation gradient Fdot, & !< field of assumed rate of compatible deformation gradient
F_tauDot !< field of assumed rate of incopatible deformation gradient F_tauDot !< field of assumed rate of incopatible deformation gradient
complex(pReal),private, dimension(:,:,:,:,:), allocatable :: inertiaField_fourier
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc. ! stress, stiffness and compliance average etc.
@ -140,6 +144,7 @@ subroutine Polarisation_init(temperature)
Utilities_constitutiveResponse, & Utilities_constitutiveResponse, &
Utilities_updateGamma, & Utilities_updateGamma, &
grid, & grid, &
grid1Red, &
geomSize, & geomSize, &
wgt wgt
use mesh, only: & use mesh, only: &
@ -174,9 +179,11 @@ subroutine Polarisation_init(temperature)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate global fields ! allocate global fields
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 (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 (F_tau_lastInc(3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal) allocate (F_tau_lastInc(3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
allocate (F_tauDot (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal) allocate (F_tauDot (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))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc Init ! PETSc Init
@ -200,6 +207,7 @@ subroutine Polarisation_init(temperature)
F_tau => xx_psc(9:17,:,:,:) F_tau => xx_psc(9:17,:,:,:)
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_lastInc2 = F_lastInc
F_tau_lastInc = F_lastInc F_tau_lastInc = F_lastInc
F = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)]) F = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)])
F_tau = F F_tau = F
@ -216,6 +224,10 @@ subroutine Polarisation_init(temperature)
trim(getSolverJobName()),size(F_lastInc)) trim(getSolverJobName()),size(F_lastInc))
read (777,rec=1) F_lastInc read (777,rec=1) F_lastInc
close (777) close (777)
call IO_read_jobBinaryFile(777,'F_lastInc2',&
trim(getSolverJobName()),size(F_lastInc2))
read (777,rec=1) F_lastInc2
close (777)
F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F
F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc
call IO_read_jobBinaryFile(777,'F_tau',& call IO_read_jobBinaryFile(777,'F_tau',&
@ -263,7 +275,7 @@ end subroutine Polarisation_init
!> @brief solution for the Polarisation scheme with internal iterations !> @brief solution for the Polarisation scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
type(tSolutionState) function & type(tSolutionState) function &
Polarisation_solution(incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,temperature_bc,rotation_BC) Polarisation_solution(incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,temperature_bc,rotation_BC,density)
use numerics, only: & use numerics, only: &
update_gamma, & update_gamma, &
itmax, & itmax, &
@ -311,6 +323,7 @@ type(tSolutionState) function &
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
real(pReal), intent(in) :: density
real(pReal) :: err_stress_tol real(pReal) :: err_stress_tol
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc Data ! PETSc Data
@ -379,7 +392,7 @@ type(tSolutionState) function &
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)]))
F_tauDot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & F_tauDot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
timeinc_old,guess,F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid(3)])) timeinc_old,guess,F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid(3)]))
F_lastInc2 = F_lastInc
F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid(3)]) F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid(3)])
F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid(3)]) F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid(3)])
endif endif
@ -407,7 +420,9 @@ type(tSolutionState) function &
params%P_BC = P_BC%values params%P_BC = P_BC%values
params%rotation_BC = rotation_BC params%rotation_BC = rotation_BC
params%timeinc = timeinc params%timeinc = timeinc
params%timeincOld = timeinc_old
params%temperature = temperature_bc params%temperature = temperature_bc
params%density = density
err_stress_tol = 1.0_pReal; err_stress = 2.0_pReal * err_stress_tol ! ensure to start loop err_stress_tol = 1.0_pReal; err_stress = 2.0_pReal * err_stress_tol ! ensure to start loop
totalIter = -1_pInt totalIter = -1_pInt
@ -461,10 +476,13 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
math_invSym3333 math_invSym3333
use DAMASK_spectral_Utilities, only: & use DAMASK_spectral_Utilities, only: &
grid, & grid, &
geomSize, &
wgt, & wgt, &
field_real, & field_real, &
field_fourier, &
Utilities_FFTforward, & Utilities_FFTforward, &
Utilities_fourierConvolution, & Utilities_fourierConvolution, &
Utilities_inverseLaplace, &
Utilities_FFTbackward, & Utilities_FFTbackward, &
Utilities_constitutiveResponse Utilities_constitutiveResponse
use debug, only: & use debug, only: &
@ -528,6 +546,21 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
math_transpose33(F_aim) math_transpose33(F_aim)
flush(6) flush(6)
endif endif
if (params%density > 0.0_pReal) then
!--------------------------------------------------------------------------------------------------
! evaluate inertia
residual_F = ((F - F_lastInc)/params%timeinc - (F_lastInc - F_lastInc2)/params%timeincOld)/((params%timeinc + params%timeincOld)/2.0_pReal)
residual_F = params%density*product(geomSize/grid)*residual_F
field_real = 0.0_pReal
field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3) = reshape(residual_F,[grid(1),grid(2),grid(3),3,3],&
order=[4,5,1,2,3]) ! field real has a different order
call Utilities_FFTforward()
call Utilities_inverseLaplace()
inertiaField_fourier = field_fourier
else
inertiaField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! !
@ -540,6 +573,7 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! doing convolution in Fourier space ! doing convolution in Fourier space
call Utilities_FFTforward() call Utilities_FFTforward()
field_fourier = field_fourier + polarAlpha*inertiaField_fourier
call Utilities_fourierConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) call Utilities_fourierConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC))
call Utilities_FFTbackward() call Utilities_FFTbackward()

37
code/DAMASK_spectral_utilities.f90
View File

@ -46,7 +46,7 @@ module DAMASK_spectral_utilities
! variables storing information for spectral method and FFTW ! variables storing information for spectral method and FFTW
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),private, 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), 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
@ -102,6 +102,7 @@ module DAMASK_spectral_utilities
utilities_FFTforward, & utilities_FFTforward, &
utilities_FFTbackward, & utilities_FFTbackward, &
utilities_fourierConvolution, & utilities_fourierConvolution, &
utilities_inverseLaplace, &
utilities_divergenceRMS, & utilities_divergenceRMS, &
utilities_curlRMS, & utilities_curlRMS, &
utilities_maskedCompliance, & utilities_maskedCompliance, &
@ -485,6 +486,40 @@ subroutine utilities_FFTbackward()
end subroutine utilities_FFTbackward end subroutine utilities_FFTbackward
!--------------------------------------------------------------------------------------------------
!> @brief doing convolution with inverse laplace kernel
!--------------------------------------------------------------------------------------------------
subroutine utilities_inverseLaplace()
use math, only: &
math_inv33, &
PI
implicit none
integer(pInt) :: i, j, k
integer(pInt), dimension(3) :: k_s
write(6,'(/,a)') ' ... doing inverse laplace .................................................'
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
if (any(k_s /= 0_pInt)) field_fourier(i,j,k, 1:3,1:3) = field_fourier(i,j,k, 1:3,1:3)/&
cmplx(-sum((2.0_pReal*PI*k_s/geomSize)*&
(2.0_pReal*PI*k_s/geomSize)),0.0_pReal,pReal) ! symmetry, junst running from 0,1,...,N/2,N/2+1
enddo; enddo; enddo
field_fourier(1,1,1,1:3,1:3) = cmplx(0.0_pReal,0.0_pReal,pReal)
end subroutine utilities_inverseLaplace
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief doing convolution gamma_hat * field_real, ensuring that average value = fieldAim !> @brief doing convolution gamma_hat * field_real, ensuring that average value = fieldAim
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------