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

Merge branch 'stress-ramp-loadcase' into development

This commit is contained in:
Sharan Roongta 2020-10-02 11:53:49 +02:00
parent 4a7701d86a 20b393ac06
commit 42d3ac045b
10 changed files with 201 additions and 238 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 25ce39dd0f5bf49dc5c2bec20767d93d2b76d353
Subproject commit 05024077bc89350d0313fc33d695a25612ac5130

45
src/grid/DAMASK_grid.f90
View File

@ -181,7 +181,7 @@ program DAMASK_grid
if ((N_def /= N_n) .or. (N_n /= N_t) .or. N_n < 1) & ! sanity check
call IO_error(error_ID=837,el=currentLoadCase,ext_msg = trim(interface_loadFile)) ! error message for incomplete loadcase
newLoadCase%stress%myType='stress'
newLoadCase%stress%myType='p'
field = 1
newLoadCase%ID(field) = FIELD_MECH_ID ! mechanical active by default
thermalActive: if (any(thermal_type == THERMAL_conduction_ID)) then
@ -210,18 +210,16 @@ program DAMASK_grid
temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not a *
if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable
enddo
newLoadCase%deformation%maskLogical = transpose(reshape(temp_maskVector,[ 3,3])) ! logical mask in 3x3 notation
newLoadCase%deformation%maskFloat = merge(ones,zeros,newLoadCase%deformation%maskLogical) ! float (1.0/0.0) mask in 3x3 notation
newLoadCase%deformation%values = math_9to33(temp_valueVector) ! values in 3x3 notation
newLoadCase%deformation%mask = transpose(reshape(temp_maskVector,[ 3,3])) ! mask in 3x3 notation
newLoadCase%deformation%values = math_9to33(temp_valueVector) ! values in 3x3 notation
case('p','stress', 's')
temp_valueVector = 0.0_pReal
do j = 1, 9
temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not an asterisk
if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable
enddo
newLoadCase%stress%maskLogical = transpose(reshape(temp_maskVector,[ 3,3]))
newLoadCase%stress%maskFloat = merge(ones,zeros,newLoadCase%stress%maskLogical)
newLoadCase%stress%values = math_9to33(temp_valueVector)
newLoadCase%stress%mask = transpose(reshape(temp_maskVector,[ 3,3]))
newLoadCase%stress%values = math_9to33(temp_valueVector)
case('t','time','delta') ! increment time
newLoadCase%time = IO_floatValue(line,chunkPos,i+1)
case('n','incs','increments') ! number of increments
@ -268,8 +266,8 @@ program DAMASK_grid
print*, ' drop guessing along trajectory'
if (newLoadCase%deformation%myType == 'l') then
do j = 1, 3
if (any(newLoadCase%deformation%maskLogical(j,1:3) .eqv. .true.) .and. &
any(newLoadCase%deformation%maskLogical(j,1:3) .eqv. .false.)) errorID = 832 ! each row should be either fully or not at all defined
if (any(newLoadCase%deformation%mask(j,1:3) .eqv. .true.) .and. &
any(newLoadCase%deformation%mask(j,1:3) .eqv. .false.)) errorID = 832 ! each row should be either fully or not at all defined
enddo
print*, ' velocity gradient:'
else if (newLoadCase%deformation%myType == 'f') then
@ -278,20 +276,19 @@ program DAMASK_grid
print*, ' deformation gradient rate:'
endif
do i = 1, 3; do j = 1, 3
if(newLoadCase%deformation%maskLogical(i,j)) then
if(newLoadCase%deformation%mask(i,j)) then
write(IO_STDOUT,'(2x,f12.7)',advance='no') newLoadCase%deformation%values(i,j)
else
write(IO_STDOUT,'(2x,12a)',advance='no') ' * '
endif
enddo; write(IO_STDOUT,'(/)',advance='no')
enddo
if (any(newLoadCase%stress%maskLogical .eqv. &
newLoadCase%deformation%maskLogical)) errorID = 831 ! exclusive or masking only
if (any(newLoadCase%stress%maskLogical .and. transpose(newLoadCase%stress%maskLogical) &
.and. (math_I3<1))) errorID = 838 ! no rotation is allowed by stress BC
if (any(newLoadCase%stress%mask .eqv. newLoadCase%deformation%mask)) errorID = 831 ! exclusive or masking only
if (any(newLoadCase%stress%mask .and. transpose(newLoadCase%stress%mask) .and. (math_I3<1))) &
errorID = 838 ! no rotation is allowed by stress BC
print*, ' stress / GPa:'
do i = 1, 3; do j = 1, 3
if(newLoadCase%stress%maskLogical(i,j)) then
if(newLoadCase%stress%mask(i,j)) then
write(IO_STDOUT,'(2x,f12.7)',advance='no') newLoadCase%stress%values(i,j)*1e-9_pReal
else
write(IO_STDOUT,'(2x,12a)',advance='no') ' * '
@ -368,11 +365,7 @@ program DAMASK_grid
timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal)
else
if (currentLoadCase == 1) then ! 1st load case of logarithmic scale
if (inc == 1) then ! 1st inc of 1st load case of logarithmic scale
timeinc = loadCases(1)%time*(2.0_pReal**real( 1-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd
else ! not-1st inc of 1st load case of logarithmic scale
timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1-loadCases(1)%incs ,pReal))
endif
timeinc = loadCases(1)%time*(2.0_pReal**real(max(inc-1,1)-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd
else ! not-1st load case of logarithmic scale
timeinc = time0 * &
( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc ,pReal)/&
@ -432,17 +425,11 @@ program DAMASK_grid
do field = 1, nActiveFields
select case(loadCases(currentLoadCase)%ID(field))
case(FIELD_MECH_ID)
solres(field) = mech_solution (&
incInfo,timeinc,timeIncOld, &
stress_BC = loadCases(currentLoadCase)%stress, &
rotation_BC = loadCases(currentLoadCase)%rot)
solres(field) = mech_solution(incInfo)
case(FIELD_THERMAL_ID)
solres(field) = grid_thermal_spectral_solution(timeinc,timeIncOld)
solres(field) = grid_thermal_spectral_solution(timeinc)
case(FIELD_DAMAGE_ID)
solres(field) = grid_damage_spectral_solution(timeinc,timeIncOld)
solres(field) = grid_damage_spectral_solution(timeinc)
end select
if (.not. solres(field)%converged) exit ! no solution found

6
src/grid/grid_damage_spectral.f90
View File

@ -156,11 +156,10 @@ end subroutine grid_damage_spectral_init
!--------------------------------------------------------------------------------------------------
!> @brief solution for the spectral damage scheme with internal iterations
!--------------------------------------------------------------------------------------------------
function grid_damage_spectral_solution(timeinc,timeinc_old) result(solution)
function grid_damage_spectral_solution(timeinc) result(solution)
real(pReal), intent(in) :: &
timeinc, & !< increment in time for current solution
timeinc_old !< increment in time of last increment
timeinc !< increment in time for current solution
integer :: i, j, k, cell
type(tSolutionState) :: solution
PetscInt :: devNull
@ -174,7 +173,6 @@ function grid_damage_spectral_solution(timeinc,timeinc_old) result(solution)
!--------------------------------------------------------------------------------------------------
! set module wide availabe data
params%timeinc = timeinc
params%timeincOld = timeinc_old
call SNESSolve(damage_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr)
call SNESGetConvergedReason(damage_snes,reason,ierr); CHKERRQ(ierr)

128
src/grid/grid_mech_FEM.f90
View File

@ -25,63 +25,56 @@ module grid_mech_FEM
implicit none
private
!--------------------------------------------------------------------------------------------------
! derived types
type(tSolutionParams), private :: params
type(tSolutionParams) :: params
type, private :: tNumerics
type :: tNumerics
integer :: &
itmin, & !< minimum number of iterations
itmax !< maximum number of iterations
real(pReal) :: &
err_div, &
divTol, &
BCTol, &
eps_div_atol, & !< absolute tolerance for equilibrium
eps_div_rtol, & !< relative tolerance for equilibrium
eps_stress_atol, & !< absolute tolerance for fullfillment of stress BC
eps_stress_rtol !< relative tolerance for fullfillment of stress BC
end type tNumerics
type(tNumerics), private :: num
logical, private:: &
debugRotation
type(tNumerics) :: num ! numerics parameters. Better name?
logical :: debugRotation
!--------------------------------------------------------------------------------------------------
! PETSc data
DM, private :: mech_grid
SNES, private :: mech_snes
Vec, private :: solution_current, solution_lastInc, solution_rate
DM :: mech_grid
SNES :: mech_snes
Vec :: solution_current, solution_lastInc, solution_rate
!--------------------------------------------------------------------------------------------------
! common pointwise data
real(pReal), private, dimension(:,:,:,:,:), allocatable :: F, P_current, F_lastInc
real(pReal), private :: detJ
real(pReal), private, dimension(3) :: delta
real(pReal), private, dimension(3,8) :: BMat
real(pReal), private, dimension(8,8) :: HGMat
PetscInt, private :: xstart,ystart,zstart,xend,yend,zend
real(pReal), dimension(:,:,:,:,:), allocatable :: F, P_current, F_lastInc
real(pReal) :: detJ
real(pReal), dimension(3) :: delta
real(pReal), dimension(3,8) :: BMat
real(pReal), dimension(8,8) :: HGMat
PetscInt :: xstart,ystart,zstart,xend,yend,zend
!--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc.
real(pReal), private, dimension(3,3) :: &
real(pReal), dimension(3,3) :: &
F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient
F_aim = math_I3, & !< current prescribed deformation gradient
F_aim_lastIter = math_I3, &
F_aim_lastInc = math_I3, & !< previous average deformation gradient
P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress
character(len=:), allocatable, private :: incInfo !< time and increment information
real(pReal), private, dimension(3,3,3,3) :: &
P_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress
P_aim = 0.0_pReal
character(len=:), allocatable :: incInfo !< time and increment information
real(pReal), dimension(3,3,3,3) :: &
C_volAvg = 0.0_pReal, & !< current volume average stiffness
C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness
S = 0.0_pReal !< current compliance (filled up with zeros)
real(pReal), private :: &
real(pReal) :: &
err_BC !< deviation from stress BC
integer, private :: &
integer :: &
totalIter = 0 !< total iteration in current increment
public :: &
@ -99,7 +92,6 @@ contains
subroutine grid_mech_FEM_init
real(pReal) :: HGCoeff = 0.0e-2_pReal
PetscInt, dimension(0:worldsize-1) :: localK
real(pReal), dimension(3,3) :: &
temp33_Real = 0.0_pReal
real(pReal), dimension(4,8) :: &
@ -111,33 +103,34 @@ subroutine grid_mech_FEM_init
-1.0_pReal, 1.0_pReal,-1.0_pReal,-1.0_pReal, &
1.0_pReal,-1.0_pReal,-1.0_pReal,-1.0_pReal, &
1.0_pReal, 1.0_pReal, 1.0_pReal, 1.0_pReal], [4,8])
real(pReal), dimension(3,3,3,3) :: devNull
PetscErrorCode :: ierr
PetscScalar, pointer, dimension(:,:,:,:) :: &
u_current,u_lastInc
PetscInt, dimension(0:worldsize-1) :: localK
integer(HID_T) :: fileHandle, groupHandle
character(len=pStringLen) :: &
fileName
class(tNode), pointer :: &
num_grid, &
debug_grid
real(pReal), dimension(3,3,3,3) :: devNull
PetscScalar, pointer, dimension(:,:,:,:) :: &
u_current,u_lastInc
print'(/,a)', ' <<<+- grid_mech_FEM init -+>>>'; flush(IO_STDOUT)
!-----------------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------------------
! debugging options
debug_grid => config_debug%get('grid', defaultVal=emptyList)
debugRotation = debug_grid%contains('rotation')
!-------------------------------------------------------------------------------------------------
! read numerical parameters and do sanity checks
num_grid => config_numerics%get('grid',defaultVal=emptyDict)
num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal)
num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal)
num%eps_stress_atol = num_grid%get_asFloat ('eps_stress_atol', defaultVal=1.0e3_pReal)
num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol', defaultVal=0.01_pReal)
num%itmin = num_grid%get_asInt ('itmin',defaultVal=1)
num%itmax = num_grid%get_asInt ('itmax',defaultVal=250)
num%eps_div_atol = num_grid%get_asFloat('eps_div_atol', defaultVal=1.0e-4_pReal)
num%eps_div_rtol = num_grid%get_asFloat('eps_div_rtol', defaultVal=5.0e-4_pReal)
num%eps_stress_atol = num_grid%get_asFloat('eps_stress_atol',defaultVal=1.0e3_pReal)
num%eps_stress_rtol = num_grid%get_asFloat('eps_stress_rtol',defaultVal=1.0e-3_pReal)
num%itmin = num_grid%get_asInt ('itmin', defaultVal=1)
num%itmax = num_grid%get_asInt ('itmax', defaultVal=250)
if (num%eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol')
if (num%eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol')
@ -242,6 +235,7 @@ subroutine grid_mech_FEM_init
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity
F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3)
endif restartRead
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
call utilities_updateCoords(F)
call utilities_constitutiveResponse(P_current,temp33_Real,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2
@ -268,19 +262,12 @@ end subroutine grid_mech_FEM_init
!--------------------------------------------------------------------------------------------------
!> @brief solution for the FEM scheme with internal iterations
!--------------------------------------------------------------------------------------------------
function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution)
function grid_mech_FEM_solution(incInfoIn) result(solution)
!--------------------------------------------------------------------------------------------------
! input data for solution
character(len=*), intent(in) :: &
incInfoIn
real(pReal), intent(in) :: &
timeinc, & !< time increment of current solution
timeinc_old !< time increment of last successful increment
type(tBoundaryCondition), intent(in) :: &
stress_BC
type(rotation), intent(in) :: &
rotation_BC
type(tSolutionState) :: &
solution
!--------------------------------------------------------------------------------------------------
@ -292,14 +279,7 @@ function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg)
!--------------------------------------------------------------------------------------------------
! set module wide available data
params%stress_mask = stress_BC%maskFloat
params%stress_BC = stress_BC%values
params%rotation_BC = rotation_BC
params%timeinc = timeinc
params%timeincOld = timeinc_old
S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg)
!--------------------------------------------------------------------------------------------------
! solve BVP
@ -341,6 +321,14 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,
PetscScalar, pointer, dimension(:,:,:,:) :: &
u_current,u_lastInc
!--------------------------------------------------------------------------------------------------
! set module wide available data
params%stress_mask = stress_BC%mask
params%rotation_BC = rotation_BC
params%timeinc = timeinc
params%timeincOld = timeinc_old
call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr)
call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr)
@ -349,20 +337,20 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,
else
C_volAvgLastInc = C_volAvg
F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess)
F_aimDot = merge(merge((F_aim-F_aim_lastInc)/timeinc_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess)
F_aim_lastInc = F_aim
!-----------------------------------------------------------------------------------------------
! calculate rate for aim
if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc)
F_aimDot = F_aimDot &
+ merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask)
elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * deformation_BC%values
F_aimDot = F_aimDot &
+ merge(deformation_BC%values,.0_pReal,deformation_BC%mask)
elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime
F_aimDot = F_aimDot &
+ merge((deformation_BC%values - F_aim_lastInc)/loadCaseTime,.0_pReal,deformation_BC%mask)
endif
if (guess) then
@ -382,6 +370,12 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,
!--------------------------------------------------------------------------------------------------
! update average and local deformation gradients
F_aim = F_aim_lastInc + F_aimDot * timeinc
if (stress_BC%myType=='p') then
P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask)
elseif (stress_BC%myType=='pdot') then !UNTESTED
P_aim = P_aim + merge(stress_BC%values*timeinc,.0_pReal,stress_BC%mask)
endif
call VecAXPY(solution_current,timeinc,solution_rate,ierr); CHKERRQ(ierr)
call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr);CHKERRQ(ierr)
@ -497,8 +491,6 @@ subroutine formResidual(da_local,x_local, &
PetscScalar, pointer,dimension(:,:,:,:) :: x_scal, f_scal
PetscScalar, dimension(8,3) :: x_elem, f_elem
PetscInt :: i, ii, j, jj, k, kk, ctr, ele
real(pReal), dimension(3,3) :: &
deltaF_aim
PetscInt :: &
PETScIter, &
nfuncs
@ -547,10 +539,8 @@ subroutine formResidual(da_local,x_local, &
!--------------------------------------------------------------------------------------------------
! stress BC handling
F_aim_lastIter = F_aim
deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC)
F_aim = F_aim - deltaF_aim
err_BC = maxval(abs(params%stress_mask * (P_av - params%stress_BC))) ! mask = 0.0 when no stress bc
F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc
err_BC = maxval(abs(merge(P_av - P_aim,.0_pReal,params%stress_mask)))
!--------------------------------------------------------------------------------------------------
! constructing residual

100
src/grid/grid_mech_spectral_basic.f90
View File

@ -24,11 +24,9 @@ module grid_mech_spectral_basic
implicit none
private
!--------------------------------------------------------------------------------------------------
! derived types
type(tSolutionParams) :: params
type, private :: tNumerics
type :: tNumerics
logical :: update_gamma !< update gamma operator with current stiffness
integer :: &
itmin, & !< minimum number of iterations
@ -42,7 +40,7 @@ module grid_mech_spectral_basic
type(tNumerics) :: num ! numerics parameters. Better name?
logical, private :: debugRotation
logical :: debugRotation
!--------------------------------------------------------------------------------------------------
! PETSc data
@ -62,10 +60,10 @@ module grid_mech_spectral_basic
F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient
F_aim = math_I3, & !< current prescribed deformation gradient
F_aim_lastInc = math_I3, & !< previous average deformation gradient
P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress
P_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress
P_aim = 0.0_pReal
character(len=:), allocatable :: incInfo !< time and increment information
real(pReal), private, dimension(3,3,3,3) :: &
real(pReal), dimension(3,3,3,3) :: &
C_volAvg = 0.0_pReal, & !< current volume average stiffness
C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness
C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 stiffness
@ -96,18 +94,17 @@ subroutine grid_mech_spectral_basic_init
real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P
real(pReal), dimension(3,3) :: &
temp33_Real = 0.0_pReal
class (tNode), pointer :: &
num_grid, &
debug_grid
PetscErrorCode :: ierr
PetscScalar, pointer, dimension(:,:,:,:) :: &
F ! pointer to solution data
PetscInt, dimension(worldsize) :: localK
PetscInt, dimension(0:worldsize-1) :: localK
integer(HID_T) :: fileHandle, groupHandle
integer :: fileUnit
character(len=pStringLen) :: &
fileName
class (tNode), pointer :: &
num_grid, &
debug_grid
print'(/,a)', ' <<<+- grid_mech_spectral_basic init -+>>>'; flush(IO_STDOUT)
@ -126,13 +123,13 @@ subroutine grid_mech_spectral_basic_init
! read numerical parameters and do sanity checks
num_grid => config_numerics%get('grid',defaultVal=emptyDict)
num%update_gamma = num_grid%get_asBool ('update_gamma', defaultVal=.false.)
num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal)
num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal)
num%eps_stress_atol = num_grid%get_asFloat ('eps_stress_atol',defaultVal=1.0e3_pReal)
num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol',defaultVal=0.01_pReal)
num%itmin = num_grid%get_asInt ('itmin',defaultVal=1)
num%itmax = num_grid%get_asInt ('itmax',defaultVal=250)
num%update_gamma = num_grid%get_asBool ('update_gamma', defaultVal=.false.)
num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal)
num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal)
num%eps_stress_atol = num_grid%get_asFloat ('eps_stress_atol',defaultVal=1.0e3_pReal)
num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol',defaultVal=1.0e-3_pReal)
num%itmin = num_grid%get_asInt ('itmin',defaultVal=1)
num%itmax = num_grid%get_asInt ('itmax',defaultVal=250)
if (num%eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol')
if (num%eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol')
@ -158,7 +155,7 @@ subroutine grid_mech_spectral_basic_init
call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr)
localK = 0
localK(worldrank+1) = grid3
localK(worldrank) = grid3
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr)
call DMDACreate3d(PETSC_COMM_WORLD, &
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary
@ -202,8 +199,8 @@ subroutine grid_mech_spectral_basic_init
endif restartRead
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
call Utilities_updateCoords(reshape(F,shape(F_lastInc)))
call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
call utilities_updateCoords(reshape(F,shape(F_lastInc)))
call utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
reshape(F,shape(F_lastInc)), & ! target F
0.0_pReal) ! time increment
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! deassociate pointer
@ -231,19 +228,12 @@ end subroutine grid_mech_spectral_basic_init
!--------------------------------------------------------------------------------------------------
!> @brief solution for the basic scheme with internal iterations
!--------------------------------------------------------------------------------------------------
function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution)
function grid_mech_spectral_basic_solution(incInfoIn) result(solution)
!--------------------------------------------------------------------------------------------------
! input data for solution
character(len=*), intent(in) :: &
incInfoIn
real(pReal), intent(in) :: &
timeinc, & !< time increment of current solution
timeinc_old !< time increment of last successful increment
type(tBoundaryCondition), intent(in) :: &
stress_BC
type(rotation), intent(in) :: &
rotation_BC
type(tSolutionState) :: &
solution
!--------------------------------------------------------------------------------------------------
@ -255,17 +245,9 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg)
S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg)
if(num%update_gamma) call utilities_updateGamma(C_minMaxAvg)
!--------------------------------------------------------------------------------------------------
! set module wide available data
params%stress_mask = stress_BC%maskFloat
params%stress_BC = stress_BC%values
params%rotation_BC = rotation_BC
params%timeinc = timeinc
params%timeincOld = timeinc_old
!--------------------------------------------------------------------------------------------------
! solve BVP
call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr)
@ -303,7 +285,14 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo
type(rotation), intent(in) :: &
rotation_BC
PetscErrorCode :: ierr
PetscScalar, dimension(:,:,:,:), pointer :: F
PetscScalar, pointer, dimension(:,:,:,:) :: F
!--------------------------------------------------------------------------------------------------
! set module wide available data
params%stress_mask = stress_BC%mask
params%rotation_BC = rotation_BC
params%timeinc = timeinc
params%timeincOld = timeinc_old
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
@ -314,20 +303,20 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo
C_volAvgLastInc = C_volAvg
C_minMaxAvgLastInc = C_minMaxAvg
F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess)
F_aimDot = merge(merge((F_aim-F_aim_lastInc)/timeinc_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess)
F_aim_lastInc = F_aim
!-----------------------------------------------------------------------------------------------
! calculate rate for aim
if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc)
F_aimDot = F_aimDot &
+ merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask)
elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * deformation_BC%values
F_aimDot = F_aimDot &
+ merge(deformation_BC%values,.0_pReal,deformation_BC%mask)
elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime
F_aimDot = F_aimDot &
+ merge((deformation_BC%values - F_aim_lastInc)/loadCaseTime,.0_pReal,deformation_BC%mask)
endif
Fdot = utilities_calculateRate(guess, &
@ -341,7 +330,13 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo
!--------------------------------------------------------------------------------------------------
! update average and local deformation gradients
F_aim = F_aim_lastInc + F_aimDot * timeinc
F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average
if (stress_BC%myType=='p') then
P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask)
elseif (stress_BC%myType=='pdot') then !UNTESTED
P_aim = P_aim + merge(stress_BC%values*timeinc,.0_pReal,stress_BC%mask)
endif
F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average
rotation_BC%rotate(F_aim,active=.true.)),[9,grid(1),grid(2),grid3])
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
@ -375,7 +370,7 @@ subroutine grid_mech_spectral_basic_restartWrite
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
print'(a)', ' writing solver data required for restart to file'; flush(IO_STDOUT)
print*, 'writing solver data required for restart to file'; flush(IO_STDOUT)
write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5'
fileHandle = HDF5_openFile(fileName,'w')
@ -469,6 +464,7 @@ subroutine formResidual(in, F, &
call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr)
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment
!--------------------------------------------------------------------------------------------------
! begin of new iteration
newIteration: if (totalIter <= PETScIter) then
@ -491,16 +487,16 @@ subroutine formResidual(in, F, &
!--------------------------------------------------------------------------------------------------
! stress BC handling
deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC)
deltaF_aim = math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc
F_aim = F_aim - deltaF_aim
err_BC = maxval(abs(params%stress_mask * (P_av - params%stress_BC))) ! mask = 0.0 when no stress bc
err_BC = maxval(abs(merge(P_av - P_aim,.0_pReal,params%stress_mask)))
!--------------------------------------------------------------------------------------------------
! updated deformation gradient using fix point algorithm of basic scheme
tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residuum ! store fPK field for subsequent FFT forward transform
call utilities_FFTtensorForward ! FFT forward of global "tensorField_real"
err_div = Utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use
err_div = utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use
call utilities_fourierGammaConvolution(params%rotation_BC%rotate(deltaF_aim,active=.true.)) ! convolution of Gamma and tensorField_fourier
call utilities_FFTtensorBackward ! FFT backward of global tensorField_fourier

107
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -15,7 +15,6 @@ module grid_mech_spectral_polarisation
use DAMASK_interface
use HDF5_utilities
use math
use rotations
use spectral_utilities
use FEsolving
use config
@ -25,8 +24,6 @@ module grid_mech_spectral_polarisation
implicit none
private
!--------------------------------------------------------------------------------------------------
! derived types
type(tSolutionParams) :: params
type :: tNumerics
@ -48,7 +45,7 @@ module grid_mech_spectral_polarisation
type(tNumerics) :: num ! numerics parameters. Better name?
logical, private :: debugRotation
logical :: debugRotation
!--------------------------------------------------------------------------------------------------
! PETSc data
@ -71,8 +68,8 @@ module grid_mech_spectral_polarisation
F_aim = math_I3, & !< current prescribed deformation gradient
F_aim_lastInc = math_I3, & !< previous average deformation gradient
F_av = 0.0_pReal, & !< average incompatible def grad field
P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress
P_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress
P_aim = 0.0_pReal
character(len=:), allocatable :: incInfo !< time and increment information
real(pReal), dimension(3,3,3,3) :: &
C_volAvg = 0.0_pReal, & !< current volume average stiffness
@ -108,10 +105,6 @@ subroutine grid_mech_spectral_polarisation_init
real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P
real(pReal), dimension(3,3) :: &
temp33_Real = 0.0_pReal
class (tNode), pointer :: &
num_grid, &
debug_grid
PetscErrorCode :: ierr
PetscScalar, pointer, dimension(:,:,:,:) :: &
FandF_tau, & ! overall pointer to solution data
@ -122,13 +115,16 @@ subroutine grid_mech_spectral_polarisation_init
integer :: fileUnit
character(len=pStringLen) :: &
fileName
class (tNode), pointer :: &
num_grid, &
debug_grid
print'(/,a)', ' <<<+- grid_mech_spectral_polarisation init -+>>>'; flush(IO_STDOUT)
print*, 'Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006'
!------------------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------------------
! debugging options
debug_grid => config_debug%get('grid',defaultVal=emptyList)
debugRotation = debug_grid%contains('rotation')
@ -136,17 +132,18 @@ subroutine grid_mech_spectral_polarisation_init
!-------------------------------------------------------------------------------------------------
! read numerical parameters and do sanity checks
num_grid => config_numerics%get('grid',defaultVal=emptyDict)
num%update_gamma = num_grid%get_asBool ('update_gamma', defaultVal=.false.)
num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal)
num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal)
num%eps_curl_atol = num_grid%get_asFloat ('eps_curl_atol', defaultVal=1.0e-10_pReal)
num%eps_curl_rtol = num_grid%get_asFloat ('eps_curl_rtol', defaultVal=5.0e-4_pReal)
num%eps_stress_atol = num_grid%get_asFloat ('eps_stress_atol', defaultVal=1.0e3_pReal)
num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol', defaultVal=0.01_pReal)
num%itmin = num_grid%get_asInt ('itmin', defaultVal=1)
num%itmax = num_grid%get_asInt ('itmax', defaultVal=250)
num%alpha = num_grid%get_asFloat ('alpha', defaultVal=1.0_pReal)
num%beta = num_grid%get_asFloat ('beta', defaultVal=1.0_pReal)
num%update_gamma = num_grid%get_asBool ('update_gamma', defaultVal=.false.)
num%eps_div_atol = num_grid%get_asFloat('eps_div_atol', defaultVal=1.0e-4_pReal)
num%eps_div_rtol = num_grid%get_asFloat('eps_div_rtol', defaultVal=5.0e-4_pReal)
num%eps_curl_atol = num_grid%get_asFloat('eps_curl_atol', defaultVal=1.0e-10_pReal)
num%eps_curl_rtol = num_grid%get_asFloat('eps_curl_rtol', defaultVal=5.0e-4_pReal)
num%eps_stress_atol = num_grid%get_asFloat('eps_stress_atol',defaultVal=1.0e3_pReal)
num%eps_stress_rtol = num_grid%get_asFloat('eps_stress_rtol',defaultVal=1.0e-3_pReal)
num%itmin = num_grid%get_asInt ('itmin', defaultVal=1)
num%itmax = num_grid%get_asInt ('itmax', defaultVal=250)
num%alpha = num_grid%get_asFloat('alpha', defaultVal=1.0_pReal)
num%beta = num_grid%get_asFloat('beta', defaultVal=1.0_pReal)
if (num%eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol')
if (num%eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol')
@ -228,8 +225,8 @@ subroutine grid_mech_spectral_polarisation_init
endif restartRead
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
call Utilities_updateCoords(reshape(F,shape(F_lastInc)))
call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
call utilities_updateCoords(reshape(F,shape(F_lastInc)))
call utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
reshape(F,shape(F_lastInc)), & ! target F
0.0_pReal) ! time increment
call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! deassociate pointer
@ -259,19 +256,12 @@ end subroutine grid_mech_spectral_polarisation_init
!--------------------------------------------------------------------------------------------------
!> @brief solution for the Polarisation scheme with internal iterations
!--------------------------------------------------------------------------------------------------
function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution)
function grid_mech_spectral_polarisation_solution(incInfoIn) result(solution)
!--------------------------------------------------------------------------------------------------
! input data for solution
character(len=*), intent(in) :: &
incInfoIn
real(pReal), intent(in) :: &
timeinc, & !< time increment of current solution
timeinc_old !< time increment of last successful increment
type(tBoundaryCondition), intent(in) :: &
stress_BC
type(rotation), intent(in) :: &
rotation_BC
type(tSolutionState) :: &
solution
!--------------------------------------------------------------------------------------------------
@ -283,21 +273,13 @@ function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old,
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg)
if (num%update_gamma) then
S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg)
if(num%update_gamma) then
call utilities_updateGamma(C_minMaxAvg)
C_scale = C_minMaxAvg
S_scale = math_invSym3333(C_minMaxAvg)
endif
!--------------------------------------------------------------------------------------------------
! set module wide available data
params%stress_mask = stress_BC%maskFloat
params%stress_BC = stress_BC%values
params%rotation_BC = rotation_BC
params%timeinc = timeinc
params%timeincOld = timeinc_old
!--------------------------------------------------------------------------------------------------
! solve BVP
call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr)
@ -335,10 +317,17 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc
type(rotation), intent(in) :: &
rotation_BC
PetscErrorCode :: ierr
PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau, F, F_tau
PetscScalar, pointer, dimension(:,:,:,:) :: FandF_tau, F, F_tau
integer :: i, j, k
real(pReal), dimension(3,3) :: F_lambda33
!--------------------------------------------------------------------------------------------------
! set module wide available data
params%stress_mask = stress_BC%mask
params%rotation_BC = rotation_BC
params%timeinc = timeinc
params%timeincOld = timeinc_old
call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
F => FandF_tau(0: 8,:,:,:)
F_tau => FandF_tau(9:17,:,:,:)
@ -350,20 +339,20 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc
C_volAvgLastInc = C_volAvg
C_minMaxAvgLastInc = C_minMaxAvg
F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess)
F_aimDot = merge(merge((F_aim-F_aim_lastInc)/timeinc_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess)
F_aim_lastInc = F_aim
!-----------------------------------------------------------------------------------------------
! calculate rate for aim
if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc)
F_aimDot = F_aimDot &
+ merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask)
elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * deformation_BC%values
F_aimDot = F_aimDot &
+ merge(deformation_BC%values,.0_pReal,deformation_BC%mask)
elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime
F_aimDot = F_aimDot &
+ merge((deformation_BC%values - F_aim_lastInc)/loadCaseTime,.0_pReal,deformation_BC%mask)
endif
Fdot = utilities_calculateRate(guess, &
@ -381,6 +370,12 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc
!--------------------------------------------------------------------------------------------------
! update average and local deformation gradients
F_aim = F_aim_lastInc + F_aimDot * timeinc
if (stress_BC%myType=='p') then
P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask)
elseif (stress_BC%myType=='pdot') then !UNTESTED
P_aim = P_aim + merge(stress_BC%values*timeinc,.0_pReal,stress_BC%mask)
endif
F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average
rotation_BC%rotate(F_aim,active=.true.)),&
[9,grid(1),grid(2),grid3])
@ -595,15 +590,15 @@ subroutine formResidual(in, FandF_tau, &
!--------------------------------------------------------------------------------------------------
! stress BC handling
F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc
err_BC = maxval(abs((1.0_pReal-params%stress_mask) * math_mul3333xx33(C_scale,F_aim &
-params%rotation_BC%rotate(F_av)) + &
params%stress_mask * (P_av-params%stress_BC))) ! mask = 0.0 for no bc
F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc
err_BC = maxval(abs(merge(P_av-P_aim, &
math_mul3333xx33(C_scale,F_aim-params%rotation_BC%rotate(F_av)),&
params%stress_mask)))
! calculate divergence
tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F !< stress field in disguise
call utilities_FFTtensorForward
err_div = Utilities_divergenceRMS() !< root mean squared error in divergence of stress
err_div = utilities_divergenceRMS() !< root mean squared error in divergence of stress
!--------------------------------------------------------------------------------------------------
! constructing residual
@ -622,7 +617,7 @@ subroutine formResidual(in, FandF_tau, &
tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F
call utilities_FFTtensorForward
err_curl = Utilities_curlRMS()
err_curl = utilities_curlRMS()
end subroutine formResidual

6
src/grid/grid_thermal_spectral.f90
View File

@ -148,11 +148,10 @@ end subroutine grid_thermal_spectral_init
!--------------------------------------------------------------------------------------------------
!> @brief solution for the spectral thermal scheme with internal iterations
!--------------------------------------------------------------------------------------------------
function grid_thermal_spectral_solution(timeinc,timeinc_old) result(solution)
function grid_thermal_spectral_solution(timeinc) result(solution)
real(pReal), intent(in) :: &
timeinc, & !< increment in time for current solution
timeinc_old !< increment in time of last increment
timeinc !< increment in time for current solution
integer :: i, j, k, cell
type(tSolutionState) :: solution
PetscInt :: devNull
@ -166,7 +165,6 @@ function grid_thermal_spectral_solution(timeinc,timeinc_old) result(solution)
!--------------------------------------------------------------------------------------------------
! set module wide availabe data
params%timeinc = timeinc
params%timeincOld = timeinc_old
call SNESSolve(thermal_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr)
call SNESGetConvergedReason(thermal_snes,reason,ierr); CHKERRQ(ierr)

31
src/grid/spectral_utilities.f90
View File

@ -47,15 +47,15 @@ module spectral_utilities
complex(C_DOUBLE_COMPLEX),public, dimension(:,:,:,:), pointer :: vectorField_fourier !< vector field fourier representation for fftw
real(C_DOUBLE), public, dimension(:,:,:), pointer :: scalarField_real !< scalar field real representation for fftw
complex(C_DOUBLE_COMPLEX),public, dimension(:,:,:), pointer :: scalarField_fourier !< scalar field fourier representation for fftw
complex(pReal), private, dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat !< gamma operator (field) for spectral method
complex(pReal), private, dimension(:,:,:,:), allocatable :: xi1st !< wave vector field for first derivatives
complex(pReal), private, dimension(:,:,:,:), allocatable :: xi2nd !< wave vector field for second derivatives
real(pReal), private, dimension(3,3,3,3) :: C_ref !< mechanic reference stiffness
complex(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat !< gamma operator (field) for spectral method
complex(pReal), dimension(:,:,:,:), allocatable :: xi1st !< wave vector field for first derivatives
complex(pReal), dimension(:,:,:,:), allocatable :: xi2nd !< wave vector field for second derivatives
real(pReal), dimension(3,3,3,3) :: C_ref !< mechanic reference stiffness
!--------------------------------------------------------------------------------------------------
! plans for FFTW
type(C_PTR), private :: &
type(C_PTR) :: &
planTensorForth, & !< FFTW MPI plan P(x) to P(k)
planTensorBack, & !< FFTW MPI plan F(k) to F(x)
planVectorForth, & !< FFTW MPI plan v(x) to v(k)
@ -65,7 +65,7 @@ module spectral_utilities
!--------------------------------------------------------------------------------------------------
! variables controlling debugging
logical, private :: &
logical :: &
debugGeneral, & !< general debugging of spectral solver
debugRotation, & !< also printing out results in lab frame
debugPETSc !< use some in debug defined options for more verbose PETSc solution
@ -82,10 +82,9 @@ module spectral_utilities
end type tSolutionState
type, public :: tBoundaryCondition !< set of parameters defining a boundary condition
real(pReal), dimension(3,3) :: values = 0.0_pReal, &
maskFloat = 0.0_pReal
logical, dimension(3,3) :: maskLogical = .false.
character(len=pStringLen) :: myType = 'None'
real(pReal), dimension(3,3) :: values = 0.0_pReal
logical, dimension(3,3) :: mask = .false.
character(len=pStringLen) :: myType = 'None'
end type tBoundaryCondition
type, public :: tLoadCase
@ -101,21 +100,21 @@ module spectral_utilities
integer(kind(FIELD_UNDEFINED_ID)), allocatable :: ID(:)
end type tLoadCase
type, public :: tSolutionParams !< @todo use here the type definition for a full loadcase
real(pReal), dimension(3,3) :: stress_mask, stress_BC
type, public :: tSolutionParams
real(pReal), dimension(3,3) :: stress_BC
logical, dimension(3,3) :: stress_mask
type(rotation) :: rotation_BC
real(pReal) :: timeinc
real(pReal) :: timeincOld
real(pReal) :: timeinc, timeincOld
end type tSolutionParams
type, private :: tNumerics
type :: tNumerics
integer :: &
divergence_correction !< scale divergence/curl calculation: [0: no correction, 1: size scaled to 1, 2: size scaled to Npoints]
logical :: &
memory_efficient !< calculate gamma operator on the fly
end type tNumerics
type(tNumerics), private :: num ! numerics parameters. Better name?
type(tNumerics) :: num ! numerics parameters. Better name?
enum, bind(c); enumerator :: &
DERIVATIVE_CONTINUOUS_ID, &

12
src/material.f90
View File

@ -2,7 +2,7 @@
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Parses material config file, either solverJobName.materialConfig or material.config
!> @brief Defines phase and homogenization
!--------------------------------------------------------------------------------------------------
module material
use prec
@ -83,7 +83,7 @@ module material
material_homogenizationAt !< homogenization ID of each element
integer, dimension(:,:), allocatable, public, target :: & ! (ip,elem) ToDo: ugly target for mapping hack
material_homogenizationMemberAt !< position of the element within its homogenization instance
integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem)
integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem)
material_phaseAt !< phase ID of each element
integer, dimension(:,:,:), allocatable, public, protected :: & ! (constituent,elem)
material_phaseMemberAt !< position of the element within its phase instance
@ -326,7 +326,7 @@ subroutine material_parseMicrostructure
constituents, & !> list of constituents
constituent, & !> constituent definition
phases, &
homogenization
homogenizations
integer, dimension(:), allocatable :: &
counterPhase, &
@ -362,14 +362,14 @@ subroutine material_parseMicrostructure
phases => config_material%get('phase')
allocate(counterPhase(phases%length),source=0)
homogenization => config_material%get('homogenization')
allocate(counterHomogenization(homogenization%length),source=0)
homogenizations => config_material%get('homogenization')
allocate(counterHomogenization(homogenizations%length),source=0)
do e = 1, discretization_nElem
microstructure => microstructures%get(discretization_microstructureAt(e))
constituents => microstructure%get('constituents')
material_homogenizationAt(e) = homogenization%getIndex(microstructure%get_asString('homogenization'))
material_homogenizationAt(e) = homogenizations%getIndex(microstructure%get_asString('homogenization'))
do i = 1, discretization_nIP
counterHomogenization(material_homogenizationAt(e)) = counterHomogenization(material_homogenizationAt(e)) + 1
material_homogenizationMemberAt(i,e) = counterHomogenization(material_homogenizationAt(e))

2
src/rotations.f90
View File

@ -1452,7 +1452,7 @@ subroutine selfTest
contains
function quaternion_equal(qu1,qu2) result(ok)
pure recursive function quaternion_equal(qu1,qu2) result(ok)
real(pReal), intent(in), dimension(4) :: qu1,qu2
logical :: ok