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

Merge branch 'grid-solver-polishing' into 'development' 

Grid solver polishing

See merge request damask/DAMASK!781
This commit is contained in:
Philip Eisenlohr 2023-07-25 08:51:15 +00:00
parent 040b58f795 039d5e0fce
commit ee8876abdb
22 changed files with 496 additions and 576 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/IO.f90
View File

@ -548,6 +548,8 @@ subroutine IO_error(error_ID,ext_msg,label1,ID1,label2,ID2)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! user errors ! user errors
case (600)
msg = 'only one source entry allowed'
case (603) case (603)
msg = 'invalid data for table' msg = 'invalid data for table'
case (610) case (610)

5
src/Marc/DAMASK_Marc.f90
View File

@ -153,11 +153,10 @@ end module DAMASK_interface
#include "../phase_mechanical_plastic_dislotungsten.f90" #include "../phase_mechanical_plastic_dislotungsten.f90"
#include "../phase_mechanical_plastic_nonlocal.f90" #include "../phase_mechanical_plastic_nonlocal.f90"
#include "../phase_mechanical_eigen.f90" #include "../phase_mechanical_eigen.f90"
#include "../phase_mechanical_eigen_cleavageopening.f90"
#include "../phase_mechanical_eigen_thermalexpansion.f90" #include "../phase_mechanical_eigen_thermalexpansion.f90"
#include "../phase_thermal.f90" #include "../phase_thermal.f90"
#include "../phase_thermal_dissipation.f90" #include "../phase_thermal_source_dissipation.f90"
#include "../phase_thermal_externalheat.f90" #include "../phase_thermal_source_externalheat.f90"
#include "../phase_damage.f90" #include "../phase_damage.f90"
#include "../phase_damage_isobrittle.f90" #include "../phase_damage_isobrittle.f90"
#include "../phase_damage_anisobrittle.f90" #include "../phase_damage_anisobrittle.f90"

477
src/grid/DAMASK_grid.f90
View File

@ -75,7 +75,7 @@ program DAMASK_grid
cutBack = .false.,& cutBack = .false.,&
sig sig
integer :: & integer :: &
i, j, m, field, & i, j, field, &
errorID = 0, & errorID = 0, &
cutBackLevel = 0, & !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$ cutBackLevel = 0, & !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$
stepFraction = 0, & !< fraction of current time interval stepFraction = 0, & !< fraction of current time interval
@ -110,15 +110,8 @@ program DAMASK_grid
load, & load, &
num_solver, & num_solver, &
num_grid, & num_grid, &
load_step, & solver
solver, &
step_bc, &
step_mech, &
step_discretization
type(tList), pointer :: & type(tList), pointer :: &
#ifdef __INTEL_LLVM_COMPILER
tensor, &
#endif
load_steps load_steps
character(len=:), allocatable :: & character(len=:), allocatable :: &
fileContent, fname fileContent, fname
@ -210,13 +203,251 @@ program DAMASK_grid
ID(field) = FIELD_DAMAGE_ID ID(field) = FIELD_DAMAGE_ID
end if damageActive end if damageActive
!--------------------------------------------------------------------------------------------------
! doing initialization depending on active solvers
call spectral_utilities_init()
do field = 2, nActiveFields
select case (ID(field))
case (FIELD_THERMAL_ID)
call grid_thermal_spectral_init(num_grid)
case (FIELD_DAMAGE_ID)
call grid_damage_spectral_init(num_grid)
end select
end do
call mechanical_init(num_grid)
call config_numerics_deallocate()
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
load_steps => load%get_list('loadstep') ! write header of output file
allocate(loadCases(load_steps%length)) ! array of load cases if (worldrank == 0) then
writeHeader: if (CLI_restartInc < 1) then
open(newunit=statUnit,file=trim(getSolverJobName())//'.sta',form='FORMATTED',status='REPLACE')
write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded StagIterationsNeeded' ! statistics file
else writeHeader
open(newunit=statUnit,file=trim(getSolverJobName())//&
'.sta',form='FORMATTED', position='APPEND', status='OLD')
end if writeHeader
end if
writeUndeformed: if (CLI_restartInc < 1) then
print'(/,1x,a)', '... saving initial configuration ..........................................'
flush(IO_STDOUT)
call materialpoint_result(0,0.0_pREAL)
end if writeUndeformed
loadCases = parseLoadSteps(load%get_list('loadstep'))
loadCaseLooping: do l = 1, size(loadCases)
t_0 = t ! load case start time
guess = loadCases(l)%estimate_rate ! change of load case? homogeneous guess for the first inc
incLooping: do inc = 1, loadCases(l)%N
totalIncsCounter = totalIncsCounter + 1
!--------------------------------------------------------------------------------------------------
! forwarding time
Delta_t_prev = Delta_t ! last time intervall that brought former inc to an end
if (dEq(loadCases(l)%r,1.0_pREAL,1.e-9_pREAL)) then ! linear scale
Delta_t = loadCases(l)%t/real(loadCases(l)%N,pREAL)
else
Delta_t = loadCases(l)%t * (loadCases(l)%r**(inc-1)-loadCases(l)%r**inc) &
/ (1.0_pREAL-loadCases(l)%r**loadCases(l)%N)
end if
Delta_t = Delta_t * real(subStepFactor,pREAL)**real(-cutBackLevel,pREAL) ! depending on cut back level, decrease time step
skipping: if (totalIncsCounter <= CLI_restartInc) then ! not yet at restart inc?
t = t + Delta_t ! just advance time, skip already performed calculation
guess = .true. ! QUESTION:why forced guessing instead of inheriting loadcase preference
else skipping
stepFraction = 0 ! fraction scaled by stepFactor**cutLevel
subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel)
t_remaining = loadCases(l)%t + t_0 - t
t = t + Delta_t ! forward target time
stepFraction = stepFraction + 1 ! count step
!--------------------------------------------------------------------------------------------------
! report beginning of new step
print'(/,1x,a)', '###########################################################################'
print'(1x,a,1x,es12.5,6(a,i0))', &
'Time', t, &
's: Increment ', inc,'/',loadCases(l)%N,&
'-', stepFraction,'/',subStepFactor**cutBackLevel,&
' of load case ', l,'/',size(loadCases)
write(incInfo,'(4(a,i0))') &
'Increment ',totalIncsCounter,'/',sum(loadCases%N),&
'-', stepFraction,'/',subStepFactor**cutBackLevel
flush(IO_STDOUT)
!--------------------------------------------------------------------------------------------------
! forward fields
do field = 1, nActiveFields
select case(ID(field))
case(FIELD_MECH_ID)
call mechanical_forward (&
cutBack,guess,Delta_t,Delta_t_prev,t_remaining, &
deformation_BC = loadCases(l)%deformation, &
stress_BC = loadCases(l)%stress, &
rotation_BC = loadCases(l)%rot)
case(FIELD_THERMAL_ID); call grid_thermal_spectral_forward(cutBack)
case(FIELD_DAMAGE_ID); call grid_damage_spectral_forward(cutBack)
end select
end do
if (.not. cutBack) call materialpoint_forward
!--------------------------------------------------------------------------------------------------
! solve fields
stagIter = 1
stagIterate = .true.
do while (stagIterate)
if (nActiveFields > 1) print'(/,1x,a,i0)', 'Staggered Iteration ',stagIter
do field = 1, nActiveFields
select case(ID(field))
case(FIELD_MECH_ID)
solres(field) = mechanical_solution(incInfo)
case(FIELD_THERMAL_ID)
solres(field) = grid_thermal_spectral_solution(Delta_t)
case(FIELD_DAMAGE_ID)
solres(field) = grid_damage_spectral_solution(Delta_t)
end select
if (.not. solres(field)%converged) exit ! no solution found
end do
stagIter = stagIter + 1
stagIterate = stagIter <= stagItMax &
.and. all(solres(:)%converged) &
.and. .not. all(solres(:)%stagConverged) ! stationary with respect to staggered iteration
end do
!--------------------------------------------------------------------------------------------------
! check solution and either advance or retry with smaller timestep
if ( (all(solres(:)%converged .and. solres(:)%stagConverged)) & ! converged
.and. .not. solres(1)%termIll) then ! and acceptable solution found
call mechanical_updateCoords()
Delta_t_prev = Delta_t
cutBack = .false.
guess = .true. ! start guessing after first converged (sub)inc
if (worldrank == 0) then
write(statUnit,*) totalIncsCounter, t, cutBackLevel, &
solres(1)%converged, solres(1)%iterationsNeeded, StagIter
flush(statUnit)
end if
elseif (cutBackLevel < maxCutBack) then ! further cutbacking tolerated?
cutBack = .true.
stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator
cutBackLevel = cutBackLevel + 1
t = t - Delta_t
Delta_t = Delta_t/real(subStepFactor,pREAL) ! cut timestep
print'(/,1x,a)', 'cutting back '
else ! no more options to continue
if (worldrank == 0) close(statUnit)
call IO_error(950)
end if
end do subStepLooping
cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc
if (all(solres(:)%converged)) then
print'(/,1x,a,1x,i0,1x,a)', 'increment', totalIncsCounter, 'converged'
else
print'(/,1x,a,1x,i0,1x,a)', 'increment', totalIncsCounter, 'NOT converged'
end if; flush(IO_STDOUT)
call MPI_Allreduce(signal_SIGUSR1,sig,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
if (mod(inc,loadCases(l)%f_out) == 0 .or. sig) then
print'(/,1x,a)', '... saving results ........................................................'
flush(IO_STDOUT)
call materialpoint_result(totalIncsCounter,t)
end if
if (sig) call signal_setSIGUSR1(.false.)
call MPI_Allreduce(signal_SIGUSR2,sig,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
if (mod(inc,loadCases(l)%f_restart) == 0 .or. sig) then
do field = 1, nActiveFields
select case (ID(field))
case(FIELD_MECH_ID)
call mechanical_restartWrite()
case(FIELD_THERMAL_ID)
call grid_thermal_spectral_restartWrite()
case(FIELD_DAMAGE_ID)
call grid_damage_spectral_restartWrite()
end select
end do
call materialpoint_restartWrite()
end if
if (sig) call signal_setSIGUSR2(.false.)
call MPI_Allreduce(signal_SIGINT,sig,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
if (sig) exit loadCaseLooping
end if skipping
end do incLooping
end do loadCaseLooping
!--------------------------------------------------------------------------------------------------
! report summary of whole calculation
print'(/,1x,a)', '###########################################################################'
if (worldrank == 0) close(statUnit)
call quit(0) ! no complains ;)
contains
subroutine getMaskedTensor(values,mask,tensor)
real(pREAL), intent(out), dimension(3,3) :: values
logical, intent(out), dimension(3,3) :: mask
type(tList), pointer :: tensor
type(tList), pointer :: row
integer :: i,j
values = 0.0_pREAL
do i = 1,3
row => tensor%get_list(i)
do j = 1,3
mask(i,j) = row%get_asStr(j) == 'x'
if (.not. mask(i,j)) values(i,j) = row%get_asReal(j)
end do
end do
end subroutine getMaskedTensor
function parseLoadsteps(load_steps) result(loadCases)
type(tList), intent(in), target :: load_steps
type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases
integer :: l,m
type(tDict), pointer :: &
load_step, &
step_bc, &
step_mech, &
step_discretization
#ifdef __INTEL_LLVM_COMPILER
type(tList), pointer :: &
tensor
#endif
allocate(loadCases(load_steps%length))
do l = 1, load_steps%length do l = 1, load_steps%length
load_step => load_steps%get_dict(l) load_step => load_steps%get_dict(l)
step_bc => load_step%get_dict('boundary_conditions') step_bc => load_step%get_dict('boundary_conditions')
step_mech => step_bc%get_dict('mechanical') step_mech => step_bc%get_dict('mechanical')
@ -316,226 +547,6 @@ program DAMASK_grid
end if reportAndCheck end if reportAndCheck
end do end do
end function parseLoadsteps
!--------------------------------------------------------------------------------------------------
! doing initialization depending on active solvers
call spectral_utilities_init()
do field = 2, nActiveFields
select case (ID(field))
case (FIELD_THERMAL_ID)
call grid_thermal_spectral_init(num_grid)
case (FIELD_DAMAGE_ID)
call grid_damage_spectral_init(num_grid)
end select
end do
call mechanical_init(num_grid)
call config_numerics_deallocate()
!--------------------------------------------------------------------------------------------------
! write header of output file
if (worldrank == 0) then
writeHeader: if (CLI_restartInc < 1) then
open(newunit=statUnit,file=trim(getSolverJobName())//'.sta',form='FORMATTED',status='REPLACE')
write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded' ! statistics file
else writeHeader
open(newunit=statUnit,file=trim(getSolverJobName())//&
'.sta',form='FORMATTED', position='APPEND', status='OLD')
end if writeHeader
end if
writeUndeformed: if (CLI_restartInc < 1) then
print'(/,1x,a)', '... saving initial configuration ..........................................'
flush(IO_STDOUT)
call materialpoint_result(0,0.0_pREAL)
end if writeUndeformed
loadCaseLooping: do l = 1, size(loadCases)
t_0 = t ! load case start time
guess = loadCases(l)%estimate_rate ! change of load case? homogeneous guess for the first inc
incLooping: do inc = 1, loadCases(l)%N
totalIncsCounter = totalIncsCounter + 1
!--------------------------------------------------------------------------------------------------
! forwarding time
Delta_t_prev = Delta_t ! last time intervall that brought former inc to an end
if (dEq(loadCases(l)%r,1.0_pREAL,1.e-9_pREAL)) then ! linear scale
Delta_t = loadCases(l)%t/real(loadCases(l)%N,pREAL)
else
Delta_t = loadCases(l)%t * (loadCases(l)%r**(inc-1)-loadCases(l)%r**inc) &
/ (1.0_pREAL-loadCases(l)%r**loadCases(l)%N)
end if
Delta_t = Delta_t * real(subStepFactor,pREAL)**real(-cutBackLevel,pREAL) ! depending on cut back level, decrease time step
skipping: if (totalIncsCounter <= CLI_restartInc) then ! not yet at restart inc?
t = t + Delta_t ! just advance time, skip already performed calculation
guess = .true. ! QUESTION:why forced guessing instead of inheriting loadcase preference
else skipping
stepFraction = 0 ! fraction scaled by stepFactor**cutLevel
subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel)
t_remaining = loadCases(l)%t + t_0 - t
t = t + Delta_t ! forward target time
stepFraction = stepFraction + 1 ! count step
!--------------------------------------------------------------------------------------------------
! report beginning of new step
print'(/,1x,a)', '###########################################################################'
print'(1x,a,1x,es12.5,6(a,i0))', &
'Time', t, &
's: Increment ', inc,'/',loadCases(l)%N,&
'-', stepFraction,'/',subStepFactor**cutBackLevel,&
' of load case ', l,'/',size(loadCases)
write(incInfo,'(4(a,i0))') &
'Increment ',totalIncsCounter,'/',sum(loadCases%N),&
'-', stepFraction,'/',subStepFactor**cutBackLevel
flush(IO_STDOUT)
!--------------------------------------------------------------------------------------------------
! forward fields
do field = 1, nActiveFields
select case(ID(field))
case(FIELD_MECH_ID)
call mechanical_forward (&
cutBack,guess,Delta_t,Delta_t_prev,t_remaining, &
deformation_BC = loadCases(l)%deformation, &
stress_BC = loadCases(l)%stress, &
rotation_BC = loadCases(l)%rot)
case(FIELD_THERMAL_ID); call grid_thermal_spectral_forward(cutBack)
case(FIELD_DAMAGE_ID); call grid_damage_spectral_forward(cutBack)
end select
end do
if (.not. cutBack) call materialpoint_forward
!--------------------------------------------------------------------------------------------------
! solve fields
stagIter = 0
stagIterate = .true.
do while (stagIterate)
do field = 1, nActiveFields
select case(ID(field))
case(FIELD_MECH_ID)
solres(field) = mechanical_solution(incInfo)
case(FIELD_THERMAL_ID)
solres(field) = grid_thermal_spectral_solution(Delta_t)
case(FIELD_DAMAGE_ID)
solres(field) = grid_damage_spectral_solution(Delta_t)
end select
if (.not. solres(field)%converged) exit ! no solution found
end do
stagIter = stagIter + 1
stagIterate = stagIter < stagItMax &
.and. all(solres(:)%converged) &
.and. .not. all(solres(:)%stagConverged) ! stationary with respect to staggered iteration
end do
!--------------------------------------------------------------------------------------------------
! check solution and either advance or retry with smaller timestep
if ( (all(solres(:)%converged .and. solres(:)%stagConverged)) & ! converged
.and. .not. solres(1)%termIll) then ! and acceptable solution found
call mechanical_updateCoords()
Delta_t_prev = Delta_t
cutBack = .false.
guess = .true. ! start guessing after first converged (sub)inc
if (worldrank == 0) then
write(statUnit,*) totalIncsCounter, t, cutBackLevel, &
solres(1)%converged, solres(1)%iterationsNeeded
flush(statUnit)
end if
elseif (cutBackLevel < maxCutBack) then ! further cutbacking tolerated?
cutBack = .true.
stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator
cutBackLevel = cutBackLevel + 1
t = t - Delta_t
Delta_t = Delta_t/real(subStepFactor,pREAL) ! cut timestep
print'(/,1x,a)', 'cutting back '
else ! no more options to continue
if (worldrank == 0) close(statUnit)
call IO_error(950)
end if
end do subStepLooping
cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc
if (all(solres(:)%converged)) then
print'(/,1x,a,1x,i0,1x,a)', 'increment', totalIncsCounter, 'converged'
else
print'(/,1x,a,1x,i0,1x,a)', 'increment', totalIncsCounter, 'NOT converged'
end if; flush(IO_STDOUT)
call MPI_Allreduce(signal_SIGUSR1,sig,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
if (mod(inc,loadCases(l)%f_out) == 0 .or. sig) then
print'(/,1x,a)', '... saving results ........................................................'
flush(IO_STDOUT)
call materialpoint_result(totalIncsCounter,t)
end if
if (sig) call signal_setSIGUSR1(.false.)
call MPI_Allreduce(signal_SIGUSR2,sig,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
if (mod(inc,loadCases(l)%f_restart) == 0 .or. sig) then
do field = 1, nActiveFields
select case (ID(field))
case(FIELD_MECH_ID)
call mechanical_restartWrite()
case(FIELD_THERMAL_ID)
call grid_thermal_spectral_restartWrite()
case(FIELD_DAMAGE_ID)
call grid_damage_spectral_restartWrite()
end select
end do
call materialpoint_restartWrite()
end if
if (sig) call signal_setSIGUSR2(.false.)
call MPI_Allreduce(signal_SIGINT,sig,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
if (sig) exit loadCaseLooping
end if skipping
end do incLooping
end do loadCaseLooping
!--------------------------------------------------------------------------------------------------
! report summary of whole calculation
print'(/,1x,a)', '###########################################################################'
if (worldrank == 0) close(statUnit)
call quit(0) ! no complains ;)
contains
subroutine getMaskedTensor(values,mask,tensor)
real(pREAL), intent(out), dimension(3,3) :: values
logical, intent(out), dimension(3,3) :: mask
type(tList), pointer :: tensor
type(tList), pointer :: row
integer :: i,j
values = 0.0_pREAL
do i = 1,3
row => tensor%get_list(i)
do j = 1,3
mask(i,j) = row%get_asStr(j) == 'x'
if (.not. mask(i,j)) values(i,j) = row%get_asReal(j)
end do
end do
end subroutine getMaskedTensor
end program DAMASK_grid end program DAMASK_grid

28
src/grid/grid_damage_spectral.f90
View File

@ -75,7 +75,6 @@ subroutine grid_damage_spectral_init(num_grid)
type(tDict), pointer, intent(in) :: num_grid type(tDict), pointer, intent(in) :: num_grid
integer(MPI_INTEGER_KIND), dimension(0:worldsize-1) :: cells3_global integer(MPI_INTEGER_KIND), dimension(0:worldsize-1) :: cells3_global
integer :: i, j, k, ce
DM :: DM_damage DM :: DM_damage
real(pREAL), dimension(:,:,:), pointer :: phi ! 0-indexed real(pREAL), dimension(:,:,:), pointer :: phi ! 0-indexed
Vec :: uBound, lBound Vec :: uBound, lBound
@ -92,7 +91,6 @@ subroutine grid_damage_spectral_init(num_grid)
petsc_options petsc_options
print'(/,1x,a)', '<<<+- grid_spectral_damage init -+>>>' print'(/,1x,a)', '<<<+- grid_spectral_damage init -+>>>'
print'(/,1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019' print'(/,1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019'
@ -140,7 +138,7 @@ subroutine grid_damage_spectral_init(num_grid)
1_pPETSCINT, 1_pPETSCINT, int(worldsize,pPETSCINT), & 1_pPETSCINT, 1_pPETSCINT, int(worldsize,pPETSCINT), &
1_pPETSCINT, 0_pPETSCINT, & ! #dof (phi, scalar), ghost boundary width (domain overlap) 1_pPETSCINT, 0_pPETSCINT, & ! #dof (phi, scalar), ghost boundary width (domain overlap)
[int(cells(1),pPetscInt)],[int(cells(2),pPetscInt)],int(cells3_global,pPETSCINT), & ! local cells [int(cells(1),pPetscInt)],[int(cells(2),pPetscInt)],int(cells3_global,pPETSCINT), & ! local cells
DM_damage,err_PETSc) ! handle, error DM_damage,err_PETSc) ! handle, error
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call DMsetFromOptions(DM_damage,err_PETSc) call DMsetFromOptions(DM_damage,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
@ -194,11 +192,7 @@ subroutine grid_damage_spectral_init(num_grid)
phi_stagInc = phi_lastInc phi_stagInc = phi_lastInc
end if restartRead end if restartRead
ce = 0 call homogenization_set_phi(reshape(phi,[product(cells(1:2))*cells3]))
do k = 0, cells3-1; do j = 0, cells(2)-1; do i = 0, cells(1)-1
ce = ce + 1
call homogenization_set_phi(phi(i,j,k),ce)
end do; end do; end do
call DMDAVecRestoreArrayF90(DM_damage,phi_PETSc,phi,err_PETSc) call DMDAVecRestoreArrayF90(DM_damage,phi_PETSc,phi,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
@ -216,7 +210,6 @@ function grid_damage_spectral_solution(Delta_t) result(solution)
real(pREAL), intent(in) :: & real(pREAL), intent(in) :: &
Delta_t !< increment in time for current solution Delta_t !< increment in time for current solution
integer :: i, j, k, ce
type(tSolutionState) :: solution type(tSolutionState) :: solution
PetscInt :: devNull PetscInt :: devNull
PetscReal :: phi_min, phi_max, stagNorm PetscReal :: phi_min, phi_max, stagNorm
@ -227,8 +220,6 @@ function grid_damage_spectral_solution(Delta_t) result(solution)
SNESConvergedReason :: reason SNESConvergedReason :: reason
solution%converged = .false.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set module wide availabe data ! set module wide availabe data
params%Delta_t = Delta_t params%Delta_t = Delta_t
@ -256,13 +247,7 @@ function grid_damage_spectral_solution(Delta_t) result(solution)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
phi_stagInc = phi phi_stagInc = phi
!-------------------------------------------------------------------------------------------------- call homogenization_set_phi(reshape(phi,[product(cells(1:2))*cells3]))
! updating damage state
ce = 0
do k = 0, cells3-1; do j = 0, cells(2)-1; do i = 0,cells(1)-1
ce = ce + 1
call homogenization_set_phi(phi(i,j,k),ce)
end do; end do; end do
call DMDAVecRestoreArrayF90(DM_damage,phi_PETSc,phi,err_PETSc) call DMDAVecRestoreArrayF90(DM_damage,phi_PETSc,phi,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
@ -283,7 +268,6 @@ subroutine grid_damage_spectral_forward(cutBack)
logical, intent(in) :: cutBack logical, intent(in) :: cutBack
integer :: i, j, k, ce
DM :: DM_damage DM :: DM_damage
real(pREAL), dimension(:,:,:), pointer :: phi ! 0-indexed real(pREAL), dimension(:,:,:), pointer :: phi ! 0-indexed
PetscErrorCode :: err_PETSc PetscErrorCode :: err_PETSc
@ -295,11 +279,7 @@ subroutine grid_damage_spectral_forward(cutBack)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
if (cutBack) then if (cutBack) then
ce = 0 call homogenization_set_phi(reshape(phi_lastInc,[product(cells(1:2))*cells3]))
do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1)
ce = ce + 1
call homogenization_set_phi(phi_lastInc(i,j,k),ce)
end do; end do; end do
phi = phi_lastInc phi = phi_lastInc
phi_stagInc = phi_lastInc phi_stagInc = phi_lastInc
else else

2
src/grid/grid_mech_spectral_basic.f90
View File

@ -365,7 +365,7 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_
if (stress_BC%myType=='dot_P') P_aim = P_aim & if (stress_BC%myType=='dot_P') P_aim = P_aim &
+ merge(.0_pREAL,stress_BC%values,stress_BC%mask)*Delta_t + merge(.0_pREAL,stress_BC%values,stress_BC%mask)*Delta_t
F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average F = reshape(utilities_forwardTensorField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average
rotation_BC%rotate(F_aim,active=.true.)),[9,cells(1),cells(2),cells3]) rotation_BC%rotate(F_aim,active=.true.)),[9,cells(1),cells(2),cells3])
call DMDAVecRestoreArrayF90(DM_mech,F_PETSc,F,err_PETSc) call DMDAVecRestoreArrayF90(DM_mech,F_PETSc,F,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)

6
src/grid/grid_mech_spectral_polarization.f90
View File

@ -408,11 +408,11 @@ subroutine grid_mechanical_spectral_polarization_forward(cutBack,guess,Delta_t,D
if (stress_BC%myType=='dot_P') P_aim = P_aim & if (stress_BC%myType=='dot_P') P_aim = P_aim &
+ merge(.0_pREAL,stress_BC%values,stress_BC%mask)*Delta_t + merge(.0_pREAL,stress_BC%values,stress_BC%mask)*Delta_t
F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average F = reshape(utilities_forwardTensorField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average
rotation_BC%rotate(F_aim,active=.true.)),& rotation_BC%rotate(F_aim,active=.true.)),&
[9,cells(1),cells(2),cells3]) [9,cells(1),cells(2),cells3])
if (guess) then if (guess) then
F_tau = reshape(Utilities_forwardField(Delta_t,F_tau_lastInc,F_taudot), & F_tau = reshape(Utilities_forwardTensorField(Delta_t,F_tau_lastInc,F_taudot), &
[9,cells(1),cells(2),cells3]) ! does not have any average value as boundary condition [9,cells(1),cells(2),cells3]) ! does not have any average value as boundary condition
else else
do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1) do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1)

29
src/grid/grid_thermal_spectral.f90
View File

@ -74,7 +74,7 @@ subroutine grid_thermal_spectral_init(num_grid)
type(tDict), pointer, intent(in) :: num_grid type(tDict), pointer, intent(in) :: num_grid
integer(MPI_INTEGER_KIND), dimension(0:worldsize-1) :: cells3_global integer(MPI_INTEGER_KIND), dimension(0:worldsize-1) :: cells3_global
integer :: i, j, k, ce integer :: ce
DM :: DM_thermal DM :: DM_thermal
real(pREAL), dimension(:,:,:), pointer :: T ! 0-indexed real(pREAL), dimension(:,:,:), pointer :: T ! 0-indexed
integer(MPI_INTEGER_KIND) :: err_MPI integer(MPI_INTEGER_KIND) :: err_MPI
@ -88,7 +88,6 @@ subroutine grid_thermal_spectral_init(num_grid)
petsc_options petsc_options
print'(/,1x,a)', '<<<+- grid_thermal_spectral init -+>>>' print'(/,1x,a)', '<<<+- grid_thermal_spectral init -+>>>'
print'(/,1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019' print'(/,1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019'
@ -171,11 +170,8 @@ subroutine grid_thermal_spectral_init(num_grid)
dotT_lastInc = 0.0_pREAL * T_lastInc dotT_lastInc = 0.0_pREAL * T_lastInc
end if restartRead end if restartRead
ce = 0 call homogenization_thermal_setField(reshape(T,[product(cells(1:2))*cells3]), &
do k = 0, cells3-1; do j = 0, cells(2)-1; do i = 0, cells(1)-1 [(0.0_pReal, ce = 1,product(cells(1:2))*cells3)])
ce = ce + 1
call homogenization_thermal_setField(T(i,j,k),0.0_pREAL,ce)
end do; end do; end do
call DMDAVecRestoreArrayF90(DM_thermal,T_PETSc,T,err_PETSc) call DMDAVecRestoreArrayF90(DM_thermal,T_PETSc,T,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
@ -204,8 +200,6 @@ function grid_thermal_spectral_solution(Delta_t) result(solution)
SNESConvergedReason :: reason SNESConvergedReason :: reason
solution%converged = .false.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set module wide availabe data ! set module wide availabe data
params%Delta_t = Delta_t params%Delta_t = Delta_t
@ -233,13 +227,8 @@ function grid_thermal_spectral_solution(Delta_t) result(solution)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
T_stagInc = T T_stagInc = T
!-------------------------------------------------------------------------------------------------- call homogenization_thermal_setField(reshape(T,[product(cells(1:2))*cells3]), &
! updating thermal state reshape(T-T_lastInc,[product(cells(1:2))*cells3])/params%Delta_t)
ce = 0
do k = 0, cells3-1; do j = 0, cells(2)-1; do i = 0, cells(1)-1
ce = ce + 1
call homogenization_thermal_setField(T(i,j,k),(T(i,j,k)-T_lastInc(i+1,j+1,k+1))/params%Delta_t,ce)
end do; end do; end do
call DMDAVecRestoreArrayF90(DM_thermal,T_PETSc,T,err_PETSc) call DMDAVecRestoreArrayF90(DM_thermal,T_PETSc,T,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
@ -260,7 +249,6 @@ subroutine grid_thermal_spectral_forward(cutBack)
logical, intent(in) :: cutBack logical, intent(in) :: cutBack
integer :: i, j, k, ce
DM :: DM_thermal DM :: DM_thermal
real(pREAL), dimension(:,:,:), pointer :: T ! 0-indexed real(pREAL), dimension(:,:,:), pointer :: T ! 0-indexed
PetscErrorCode :: err_PETSc PetscErrorCode :: err_PETSc
@ -272,11 +260,8 @@ subroutine grid_thermal_spectral_forward(cutBack)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
if (cutBack) then if (cutBack) then
ce = 0 call homogenization_thermal_setField(reshape(T_lastInc,[product(cells(1:2))*cells3]), &
do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) reshape(dotT_lastInc,[product(cells(1:2))*cells3]))
ce = ce + 1
call homogenization_thermal_setField(T_lastInc(i,j,k),dotT_lastInc(i,j,k),ce)
end do; end do; end do
T = T_lastInc T = T_lastInc
T_stagInc = T_lastInc T_stagInc = T_lastInc
else else

16
src/grid/spectral_utilities.f90
View File

@ -124,7 +124,7 @@ module spectral_utilities
utilities_maskedCompliance, & utilities_maskedCompliance, &
utilities_constitutiveResponse, & utilities_constitutiveResponse, &
utilities_calculateRate, & utilities_calculateRate, &
utilities_forwardField, & utilities_forwardTensorField, &
utilities_updateCoords utilities_updateCoords
contains contains
@ -864,7 +864,7 @@ end function utilities_calculateRate
!> @brief forwards a field with a pointwise given rate, if aim is given, !> @brief forwards a field with a pointwise given rate, if aim is given,
!> ensures that the average matches the aim !> ensures that the average matches the aim
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function utilities_forwardField(Delta_t,field_lastInc,rate,aim) function utilities_forwardTensorField(Delta_t,field_lastInc,rate,aim)
real(pREAL), intent(in) :: & real(pREAL), intent(in) :: &
Delta_t !< Delta_t of current step Delta_t !< Delta_t of current step
@ -875,22 +875,22 @@ function utilities_forwardField(Delta_t,field_lastInc,rate,aim)
aim !< average field value aim aim !< average field value aim
real(pREAL), dimension(3,3,cells(1),cells(2),cells3) :: & real(pREAL), dimension(3,3,cells(1),cells(2),cells3) :: &
utilities_forwardField utilities_forwardTensorField
real(pREAL), dimension(3,3) :: fieldDiff !< <a + adot*t> - aim real(pREAL), dimension(3,3) :: fieldDiff !< <a + adot*t> - aim
integer(MPI_INTEGER_KIND) :: err_MPI integer(MPI_INTEGER_KIND) :: err_MPI
utilities_forwardField = field_lastInc + rate*Delta_t utilities_forwardTensorField = field_lastInc + rate*Delta_t
if (present(aim)) then !< correct to match average if (present(aim)) then !< correct to match average
fieldDiff = sum(sum(sum(utilities_forwardField,dim=5),dim=4),dim=3)*wgt fieldDiff = sum(sum(sum(utilities_forwardTensorField,dim=5),dim=4),dim=3)*wgt
call MPI_Allreduce(MPI_IN_PLACE,fieldDiff,9_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI) call MPI_Allreduce(MPI_IN_PLACE,fieldDiff,9_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
fieldDiff = fieldDiff - aim fieldDiff = fieldDiff - aim
utilities_forwardField = utilities_forwardField & utilities_forwardTensorField = utilities_forwardTensorField &
- spread(spread(spread(fieldDiff,3,cells(1)),4,cells(2)),5,cells3) - spread(spread(spread(fieldDiff,3,cells(1)),4,cells(2)),5,cells3)
end if end if
end function utilities_forwardField end function utilities_forwardTensorField
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

11
src/homogenization.f90
View File

@ -145,9 +145,8 @@ module homogenization
real(pREAL) :: f real(pREAL) :: f
end function homogenization_f_T end function homogenization_f_T
module subroutine homogenization_thermal_setField(T,dot_T, ce) module subroutine homogenization_thermal_setField(T,dot_T)
integer, intent(in) :: ce real(pREAL), dimension(:), intent(in) :: T, dot_T
real(pREAL), intent(in) :: T, dot_T
end subroutine homogenization_thermal_setField end subroutine homogenization_thermal_setField
module function homogenization_damage_active() result(active) module function homogenization_damage_active() result(active)
@ -170,10 +169,8 @@ module homogenization
real(pREAL) :: f real(pREAL) :: f
end function homogenization_f_phi end function homogenization_f_phi
module subroutine homogenization_set_phi(phi,ce) module subroutine homogenization_set_phi(phi)
integer, intent(in) :: ce real(pREAL), dimension(:), intent(in) :: phi
real(pREAL), intent(in) :: &
phi
end subroutine homogenization_set_phi end subroutine homogenization_set_phi
end interface end interface

19
src/homogenization_damage.f90
View File

@ -151,20 +151,19 @@ end function homogenization_f_phi
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Set damage field. !> @brief Set damage field.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine homogenization_set_phi(phi,ce) module subroutine homogenization_set_phi(phi)
integer, intent(in) :: ce real(pREAL), dimension(:), intent(in) :: phi
real(pREAL), intent(in) :: phi
integer :: & integer :: ho, en, ce
ho, &
en
ho = material_ID_homogenization(ce) do ce=lbound(phi,1), ubound(phi,1)
en = material_entry_homogenization(ce) ho = material_ID_homogenization(ce)
damagestate_h(ho)%state(1,en) = phi en = material_entry_homogenization(ce)
current(ho)%phi(en) = phi damagestate_h(ho)%state(1,en) = phi(ce)
current(ho)%phi(en) = phi(ce)
end do
end subroutine homogenization_set_phi end subroutine homogenization_set_phi

17
src/homogenization_thermal.f90
View File

@ -173,15 +173,20 @@ end function homogenization_f_T
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Set thermal field and its rate (T and dot_T). !> @brief Set thermal field and its rate (T and dot_T).
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine homogenization_thermal_setField(T,dot_T, ce) module subroutine homogenization_thermal_setField(T,dot_T)
integer, intent(in) :: ce real(pREAL), dimension(:), intent(in) :: T, dot_T
real(pREAL), intent(in) :: T, dot_T
integer :: ho, en, ce
current(material_ID_homogenization(ce))%T(material_entry_homogenization(ce)) = T do ce=max(lbound(T,1),lbound(dot_T,1)), min(ubound(T,1),ubound(dot_T,1))
current(material_ID_homogenization(ce))%dot_T(material_entry_homogenization(ce)) = dot_T ho = material_ID_homogenization(ce)
call thermal_partition(ce) en = material_entry_homogenization(ce)
current(ho)%T(en) = T(ce)
current(ho)%dot_T(en) = dot_T(ce)
call thermal_partition(ce)
end do
end subroutine homogenization_thermal_setField end subroutine homogenization_thermal_setField

23
src/phase.f90
View File

@ -49,6 +49,29 @@ module phase
type(tState), dimension(:), allocatable :: p !< tState for each active source mechanism in a phase type(tState), dimension(:), allocatable :: p !< tState for each active source mechanism in a phase
end type end type
enum, bind(c); enumerator :: &
UNDEFINED, &
MECHANICAL_PLASTICITY_NONE, &
MECHANICAL_PLASTICITY_ISOTROPIC, &
MECHANICAL_PLASTICITY_PHENOPOWERLAW, &
MECHANICAL_PLASTICITY_KINEHARDENING, &
MECHANICAL_PLASTICITY_DISLOTWIN, &
MECHANICAL_PLASTICITY_DISLOTUNGSTEN, &
MECHANICAL_PLASTICITY_NONLOCAL, &
MECHANICAL_EIGEN_THERMALEXPANSION, &
DAMAGE_ISOBRITTLE, &
DAMAGE_ANISOBRITTLE, &
THERMAL_SOURCE_DISSIPATION, &
THERMAL_SOURCE_EXTERNALHEAT
end enum
integer(kind(UNDEFINED)), dimension(:), allocatable :: &
mechanical_plasticity_type, & !< plasticity of each phase
damage_type !< damage type of each phase
integer(kind(UNDEFINED)), dimension(:,:), allocatable :: &
thermal_source_type, &
mechanical_eigen_kinematics_type
character(len=2), allocatable, dimension(:) :: phase_lattice character(len=2), allocatable, dimension(:) :: phase_lattice
real(pREAL), allocatable, dimension(:) :: phase_cOverA real(pREAL), allocatable, dimension(:) :: phase_cOverA

56
src/phase_damage.f90
View File

@ -9,23 +9,14 @@ submodule(phase) damage
l_c = 0.0_pREAL !< characteristic length l_c = 0.0_pREAL !< characteristic length
end type tDamageParameters end type tDamageParameters
enum, bind(c); enumerator :: &
DAMAGE_UNDEFINED_ID, &
DAMAGE_ISOBRITTLE_ID, &
DAMAGE_ANISOBRITTLE_ID
end enum
integer :: phase_damage_maxSizeDotState integer :: phase_damage_maxSizeDotState
type :: tFieldQuantities
type :: tDataContainer
real(pREAL), dimension(:), allocatable :: phi real(pREAL), dimension(:), allocatable :: phi
end type tDataContainer end type tFieldQuantities
integer(kind(DAMAGE_UNDEFINED_ID)), dimension(:), allocatable :: &
phase_damage !< active sources mechanisms of each phase
type(tDataContainer), dimension(:), allocatable :: current type(tFieldQuantities), dimension(:), allocatable :: current
type(tDamageParameters), dimension(:), allocatable :: param type(tDamageParameters), dimension(:), allocatable :: param
@ -114,11 +105,11 @@ module subroutine damage_init()
end do end do
allocate(phase_damage(phases%length), source = DAMAGE_UNDEFINED_ID) allocate(damage_type(phases%length), source = UNDEFINED)
if (damage_active) then if (damage_active) then
where(isobrittle_init() ) phase_damage = DAMAGE_ISOBRITTLE_ID where(isobrittle_init() ) damage_type = DAMAGE_ISOBRITTLE
where(anisobrittle_init()) phase_damage = DAMAGE_ANISOBRITTLE_ID where(anisobrittle_init()) damage_type = DAMAGE_ANISOBRITTLE
end if end if
phase_damage_maxSizeDotState = maxval(damageState%sizeDotState) phase_damage_maxSizeDotState = maxval(damageState%sizeDotState)
@ -159,8 +150,8 @@ module function phase_damage_C66(C66,ph,en) result(C66_degraded)
real(pREAL), dimension(6,6) :: C66_degraded real(pREAL), dimension(6,6) :: C66_degraded
damageType: select case (phase_damage(ph)) damageType: select case (damage_type(ph))
case (DAMAGE_ISOBRITTLE_ID) damageType case (DAMAGE_ISOBRITTLE) damageType
C66_degraded = C66 * damage_phi(ph,en)**2 C66_degraded = C66 * damage_phi(ph,en)**2
case default damageType case default damageType
C66_degraded = C66 C66_degraded = C66
@ -204,13 +195,14 @@ module function phase_f_phi(phi,co,ce) result(f)
ph, & ph, &
en en
ph = material_ID_phase(co,ce) ph = material_ID_phase(co,ce)
en = material_entry_phase(co,ce) en = material_entry_phase(co,ce)
select case(phase_damage(ph)) select case(damage_type(ph))
case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ANISOBRITTLE_ID) case(DAMAGE_ISOBRITTLE,DAMAGE_ANISOBRITTLE)
f = 1.0_pREAL & f = 1.0_pREAL &
- 2.0_pREAL * phi*damageState(ph)%state(1,en) - 2.0_pREAL * phi*damageState(ph)%state(1,en) ! ToDo: MD: seems to be phi**2
case default case default
f = 0.0_pREAL f = 0.0_pREAL
end select end select
@ -318,8 +310,8 @@ module subroutine damage_restartWrite(groupHandle,ph)
integer, intent(in) :: ph integer, intent(in) :: ph
select case(phase_damage(ph)) select case(damage_type(ph))
case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ANISOBRITTLE_ID) case(DAMAGE_ISOBRITTLE,DAMAGE_ANISOBRITTLE)
call HDF5_write(damageState(ph)%state,groupHandle,'omega_damage') call HDF5_write(damageState(ph)%state,groupHandle,'omega_damage')
end select end select
@ -332,8 +324,8 @@ module subroutine damage_restartRead(groupHandle,ph)
integer, intent(in) :: ph integer, intent(in) :: ph
select case(phase_damage(ph)) select case(damage_type(ph))
case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ANISOBRITTLE_ID) case(DAMAGE_ISOBRITTLE,DAMAGE_ANISOBRITTLE)
call HDF5_read(damageState(ph)%state0,groupHandle,'omega_damage') call HDF5_read(damageState(ph)%state0,groupHandle,'omega_damage')
end select end select
@ -350,15 +342,15 @@ module subroutine damage_result(group,ph)
integer, intent(in) :: ph integer, intent(in) :: ph
if (phase_damage(ph) /= DAMAGE_UNDEFINED_ID) & if (damage_type(ph) /= UNDEFINED) &
call result_closeGroup(result_addGroup(group//'damage')) call result_closeGroup(result_addGroup(group//'damage'))
sourceType: select case (phase_damage(ph)) sourceType: select case (damage_type(ph))
case (DAMAGE_ISOBRITTLE_ID) sourceType case (DAMAGE_ISOBRITTLE) sourceType
call isobrittle_result(ph,group//'damage/') call isobrittle_result(ph,group//'damage/')
case (DAMAGE_ANISOBRITTLE_ID) sourceType case (DAMAGE_ANISOBRITTLE) sourceType
call anisobrittle_result(ph,group//'damage/') call anisobrittle_result(ph,group//'damage/')
end select sourceType end select sourceType
@ -381,9 +373,9 @@ function phase_damage_collectDotState(ph,en) result(broken)
if (damageState(ph)%sizeState > 0) then if (damageState(ph)%sizeState > 0) then
sourceType: select case (phase_damage(ph)) sourceType: select case (damage_type(ph))
case (DAMAGE_ANISOBRITTLE_ID) sourceType case (DAMAGE_ANISOBRITTLE) sourceType
call anisobrittle_dotState(mechanical_S(ph,en), ph,en) ! ToDo: use M_d call anisobrittle_dotState(mechanical_S(ph,en), ph,en) ! ToDo: use M_d
end select sourceType end select sourceType
@ -446,9 +438,9 @@ function phase_damage_deltaState(Fe, ph, en) result(broken)
if (damageState(ph)%sizeState == 0) return if (damageState(ph)%sizeState == 0) return
sourceType: select case (phase_damage(ph)) sourceType: select case (damage_type(ph))
case (DAMAGE_ISOBRITTLE_ID) sourceType case (DAMAGE_ISOBRITTLE) sourceType
call isobrittle_deltaState(phase_homogenizedC66(ph,en), Fe, ph,en) call isobrittle_deltaState(phase_homogenizedC66(ph,en), Fe, ph,en)
broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,en))) broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,en)))
if (.not. broken) then if (.not. broken) then

34
src/phase_mechanical.f90
View File

@ -3,21 +3,6 @@
!---------------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------------
submodule(phase) mechanical submodule(phase) mechanical
enum, bind(c); enumerator :: &
PLASTIC_UNDEFINED_ID, &
PLASTIC_NONE_ID, &
PLASTIC_ISOTROPIC_ID, &
PLASTIC_PHENOPOWERLAW_ID, &
PLASTIC_KINEHARDENING_ID, &
PLASTIC_DISLOTWIN_ID, &
PLASTIC_DISLOTUNGSTEN_ID, &
PLASTIC_NONLOCAL_ID, &
EIGEN_UNDEFINED_ID, &
EIGEN_CLEAVAGE_OPENING_ID, &
EIGEN_THERMAL_EXPANSION_ID
end enum
type(tTensorContainer), dimension(:), allocatable :: & type(tTensorContainer), dimension(:), allocatable :: &
! current value ! current value
phase_mechanical_Fe, & phase_mechanical_Fe, &
@ -37,9 +22,6 @@ submodule(phase) mechanical
phase_mechanical_S0 phase_mechanical_S0
integer(kind(PLASTIC_undefined_ID)), dimension(:), allocatable :: &
phase_plasticity !< plasticity of each phase
interface interface
module subroutine eigen_init(phases) module subroutine eigen_init(phases)
@ -283,7 +265,7 @@ module subroutine mechanical_init(phases)
call elastic_init(phases) call elastic_init(phases)
allocate(plasticState(phases%length)) allocate(plasticState(phases%length))
allocate(phase_plasticity(phases%length),source = PLASTIC_UNDEFINED_ID) allocate(mechanical_plasticity_type(phases%length),source = UNDEFINED)
call plastic_init() call plastic_init()
do ph = 1,phases%length do ph = 1,phases%length
plasticState(ph)%state0 = plasticState(ph)%state plasticState(ph)%state0 = plasticState(ph)%state
@ -327,24 +309,24 @@ module subroutine mechanical_result(group,ph)
call results(group,ph) call results(group,ph)
select case(phase_plasticity(ph)) select case(mechanical_plasticity_type(ph))
case(PLASTIC_ISOTROPIC_ID) case(MECHANICAL_PLASTICITY_ISOTROPIC)
call plastic_isotropic_result(ph,group//'mechanical/') call plastic_isotropic_result(ph,group//'mechanical/')
case(PLASTIC_PHENOPOWERLAW_ID) case(MECHANICAL_PLASTICITY_PHENOPOWERLAW)
call plastic_phenopowerlaw_result(ph,group//'mechanical/') call plastic_phenopowerlaw_result(ph,group//'mechanical/')
case(PLASTIC_KINEHARDENING_ID) case(MECHANICAL_PLASTICITY_KINEHARDENING)
call plastic_kinehardening_result(ph,group//'mechanical/') call plastic_kinehardening_result(ph,group//'mechanical/')
case(PLASTIC_DISLOTWIN_ID) case(MECHANICAL_PLASTICITY_DISLOTWIN)
call plastic_dislotwin_result(ph,group//'mechanical/') call plastic_dislotwin_result(ph,group//'mechanical/')
case(PLASTIC_DISLOTUNGSTEN_ID) case(MECHANICAL_PLASTICITY_DISLOTUNGSTEN)
call plastic_dislotungsten_result(ph,group//'mechanical/') call plastic_dislotungsten_result(ph,group//'mechanical/')
case(PLASTIC_NONLOCAL_ID) case(MECHANICAL_PLASTICITY_NONLOCAL)
call plastic_nonlocal_result(ph,group//'mechanical/') call plastic_nonlocal_result(ph,group//'mechanical/')
end select end select

43
src/phase_mechanical_eigen.f90
View File

@ -3,15 +3,7 @@ submodule(phase:mechanical) eigen
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
Nmodels Nmodels
integer(kind(EIGEN_UNDEFINED_ID)), dimension(:,:), allocatable :: &
model
integer(kind(EIGEN_UNDEFINED_ID)), dimension(:), allocatable :: &
model_damage
interface interface
module function damage_anisobrittle_init() result(myKinematics)
logical, dimension(:), allocatable :: myKinematics
end function damage_anisobrittle_init
module function thermalexpansion_init(kinematics_length) result(myKinematics) module function thermalexpansion_init(kinematics_length) result(myKinematics)
integer, intent(in) :: kinematics_length integer, intent(in) :: kinematics_length
@ -60,17 +52,12 @@ module subroutine eigen_init(phases)
Nmodels(ph) = kinematics%length Nmodels(ph) = kinematics%length
end do end do
allocate(model(maxval(Nmodels),phases%length), source = EIGEN_undefined_ID) allocate(mechanical_eigen_kinematics_type(maxval(Nmodels),phases%length), source = UNDEFINED)
if (maxval(Nmodels) /= 0) then if (maxval(Nmodels) /= 0) then
where(thermalexpansion_init(maxval(Nmodels))) model = EIGEN_thermal_expansion_ID where(thermalexpansion_init(maxval(Nmodels))) mechanical_eigen_kinematics_type = MECHANICAL_EIGEN_THERMALEXPANSION
end if end if
allocate(model_damage(phases%length), source = EIGEN_UNDEFINED_ID)
where(damage_anisobrittle_init()) model_damage = EIGEN_cleavage_opening_ID
end subroutine eigen_init end subroutine eigen_init
@ -173,17 +160,9 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
dLi_dFi = 0.0_pREAL dLi_dFi = 0.0_pREAL
plasticType: select case (phase_plasticity(ph))
case (PLASTIC_isotropic_ID) plasticType
call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,ph,en)
Li = Li + my_Li
dLi_dS = dLi_dS + my_dLi_dS
active = .true.
end select plasticType
KinematicsLoop: do k = 1, Nmodels(ph) KinematicsLoop: do k = 1, Nmodels(ph)
kinematicsType: select case (model(k,ph)) kinematicsType: select case (mechanical_eigen_kinematics_type(k,ph))
case (EIGEN_thermal_expansion_ID) kinematicsType case (MECHANICAL_EIGEN_THERMALEXPANSION) kinematicsType
call thermalexpansion_LiAndItsTangent(my_Li, my_dLi_dS, ph,en) call thermalexpansion_LiAndItsTangent(my_Li, my_dLi_dS, ph,en)
Li = Li + my_Li Li = Li + my_Li
dLi_dS = dLi_dS + my_dLi_dS dLi_dS = dLi_dS + my_dLi_dS
@ -191,13 +170,21 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
end select kinematicsType end select kinematicsType
end do KinematicsLoop end do KinematicsLoop
select case (model_damage(ph)) plasticType: select case (mechanical_plasticity_type(ph))
case (EIGEN_cleavage_opening_ID) case (MECHANICAL_PLASTICITY_ISOTROPIC) plasticType
call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,ph,en)
Li = Li + my_Li
dLi_dS = dLi_dS + my_dLi_dS
active = .true.
end select plasticType
damageType: select case (damage_type(ph))
case (DAMAGE_ANISOBRITTLE)
call damage_anisobrittle_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, en) call damage_anisobrittle_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, en)
Li = Li + my_Li Li = Li + my_Li
dLi_dS = dLi_dS + my_dLi_dS dLi_dS = dLi_dS + my_dLi_dS
active = .true. active = .true.
end select end select damageType
if (.not. active) return if (.not. active) return

30
src/phase_mechanical_eigen_cleavageopening.f90
View File

@ -1,30 +0,0 @@
!--------------------------------------------------------------------------------------------------
!> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine incorporating kinematics resulting from opening of cleavage planes
!> @details to be done
!--------------------------------------------------------------------------------------------------
submodule(phase:eigen) cleavageopening
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
module function damage_anisobrittle_init() result(myKinematics)
logical, dimension(:), allocatable :: myKinematics
myKinematics = kinematics_active2('anisobrittle')
if (count(myKinematics) == 0) return
print'(/,1x,a)', '<<<+- phase:mechanical:eigen:cleavageopening init -+>>>'
print'(/,a,i2)', ' # phases: ',count(myKinematics); flush(IO_STDOUT)
end function damage_anisobrittle_init
end submodule cleavageopening

4
src/phase_mechanical_elastic.f90
View File

@ -199,8 +199,8 @@ module function phase_homogenizedC66(ph,en) result(C)
integer, intent(in) :: ph, en integer, intent(in) :: ph, en
plasticType: select case (phase_plasticity(ph)) plasticType: select case (mechanical_plasticity_type(ph))
case (PLASTIC_DISLOTWIN_ID) plasticType case (MECHANICAL_PLASTICITY_DISLOTWIN) plasticType
C = plastic_dislotwin_homogenizedC(ph,en) C = plastic_dislotwin_homogenizedC(ph,en)
case default plasticType case default plasticType
C = elastic_C66(ph,en) C = elastic_C66(ph,en)

68
src/phase_mechanical_plastic.f90
View File

@ -211,17 +211,17 @@ contains
module subroutine plastic_init module subroutine plastic_init
print'(/,1x,a)', '<<<+- phase:mechanical:plastic init -+>>>' print'(/,1x,a)', '<<<+- phase:mechanical:plasticity init -+>>>'
where(plastic_none_init()) phase_plasticity = PLASTIC_NONE_ID where(plastic_none_init()) mechanical_plasticity_type = MECHANICAL_PLASTICITY_NONE
where(plastic_isotropic_init()) phase_plasticity = PLASTIC_ISOTROPIC_ID where(plastic_isotropic_init()) mechanical_plasticity_type = MECHANICAL_PLASTICITY_ISOTROPIC
where(plastic_phenopowerlaw_init()) phase_plasticity = PLASTIC_PHENOPOWERLAW_ID where(plastic_phenopowerlaw_init()) mechanical_plasticity_type = MECHANICAL_PLASTICITY_PHENOPOWERLAW
where(plastic_kinehardening_init()) phase_plasticity = PLASTIC_KINEHARDENING_ID where(plastic_kinehardening_init()) mechanical_plasticity_type = MECHANICAL_PLASTICITY_KINEHARDENING
where(plastic_dislotwin_init()) phase_plasticity = PLASTIC_DISLOTWIN_ID where(plastic_dislotwin_init()) mechanical_plasticity_type = MECHANICAL_PLASTICITY_DISLOTWIN
where(plastic_dislotungsten_init()) phase_plasticity = PLASTIC_DISLOTUNGSTEN_ID where(plastic_dislotungsten_init()) mechanical_plasticity_type = MECHANICAL_PLASTICITY_DISLOTUNGSTEN
where(plastic_nonlocal_init()) phase_plasticity = PLASTIC_NONLOCAL_ID where(plastic_nonlocal_init()) mechanical_plasticity_type = MECHANICAL_PLASTICITY_NONLOCAL
if (any(phase_plasticity == PLASTIC_undefined_ID)) call IO_error(201) if (any(mechanical_plasticity_type == UNDEFINED)) call IO_error(201)
end subroutine plastic_init end subroutine plastic_init
@ -251,7 +251,7 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
i, j i, j
if (phase_plasticity(ph) == PLASTIC_NONE_ID) then if (mechanical_plasticity_type(ph) == MECHANICAL_PLASTICITY_NONE) then
Lp = 0.0_pREAL Lp = 0.0_pREAL
dLp_dFi = 0.0_pREAL dLp_dFi = 0.0_pREAL
dLp_dS = 0.0_pREAL dLp_dS = 0.0_pREAL
@ -259,24 +259,24 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
Mp = matmul(matmul(transpose(Fi),Fi),S) Mp = matmul(matmul(transpose(Fi),Fi),S)
plasticType: select case (phase_plasticity(ph)) plasticType: select case (mechanical_plasticity_type(ph))
case (PLASTIC_ISOTROPIC_ID) plasticType case (MECHANICAL_PLASTICITY_ISOTROPIC) plasticType
call isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) call isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
case (PLASTIC_PHENOPOWERLAW_ID) plasticType case (MECHANICAL_PLASTICITY_PHENOPOWERLAW) plasticType
call phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) call phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
case (PLASTIC_KINEHARDENING_ID) plasticType case (MECHANICAL_PLASTICITY_KINEHARDENING) plasticType
call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
case (PLASTIC_NONLOCAL_ID) plasticType case (MECHANICAL_PLASTICITY_NONLOCAL) plasticType
call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
case (PLASTIC_DISLOTWIN_ID) plasticType case (MECHANICAL_PLASTICITY_DISLOTWIN) plasticType
call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
case (PLASTIC_DISLOTUNGSTEN_ID) plasticType case (MECHANICAL_PLASTICITY_DISLOTUNGSTEN) plasticType
call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
end select plasticType end select plasticType
@ -308,28 +308,28 @@ module function plastic_dotState(subdt,ph,en) result(dotState)
dotState dotState
if (phase_plasticity(ph) /= PLASTIC_NONE_ID) then if (mechanical_plasticity_type(ph) /= MECHANICAL_PLASTICITY_NONE) then
Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,en)),& Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,en)),&
phase_mechanical_Fi(ph)%data(1:3,1:3,en)),phase_mechanical_S(ph)%data(1:3,1:3,en)) phase_mechanical_Fi(ph)%data(1:3,1:3,en)),phase_mechanical_S(ph)%data(1:3,1:3,en))
plasticType: select case (phase_plasticity(ph)) plasticType: select case (mechanical_plasticity_type(ph))
case (PLASTIC_ISOTROPIC_ID) plasticType case (MECHANICAL_PLASTICITY_ISOTROPIC) plasticType
dotState = isotropic_dotState(Mp,ph,en) dotState = isotropic_dotState(Mp,ph,en)
case (PLASTIC_PHENOPOWERLAW_ID) plasticType case (MECHANICAL_PLASTICITY_PHENOPOWERLAW) plasticType
dotState = phenopowerlaw_dotState(Mp,ph,en) dotState = phenopowerlaw_dotState(Mp,ph,en)
case (PLASTIC_KINEHARDENING_ID) plasticType case (MECHANICAL_PLASTICITY_KINEHARDENING) plasticType
dotState = plastic_kinehardening_dotState(Mp,ph,en) dotState = plastic_kinehardening_dotState(Mp,ph,en)
case (PLASTIC_DISLOTWIN_ID) plasticType case (MECHANICAL_PLASTICITY_DISLOTWIN) plasticType
dotState = dislotwin_dotState(Mp,ph,en) dotState = dislotwin_dotState(Mp,ph,en)
case (PLASTIC_DISLOTUNGSTEN_ID) plasticType case (MECHANICAL_PLASTICITY_DISLOTUNGSTEN) plasticType
dotState = dislotungsten_dotState(Mp,ph,en) dotState = dislotungsten_dotState(Mp,ph,en)
case (PLASTIC_NONLOCAL_ID) plasticType case (MECHANICAL_PLASTICITY_NONLOCAL) plasticType
call nonlocal_dotState(Mp,subdt,ph,en) call nonlocal_dotState(Mp,subdt,ph,en)
dotState = plasticState(ph)%dotState(:,en) dotState = plasticState(ph)%dotState(:,en)
@ -349,15 +349,15 @@ module subroutine plastic_dependentState(ph,en)
en en
plasticType: select case (phase_plasticity(ph)) plasticType: select case (mechanical_plasticity_type(ph))
case (PLASTIC_DISLOTWIN_ID) plasticType case (MECHANICAL_PLASTICITY_DISLOTWIN) plasticType
call dislotwin_dependentState(ph,en) call dislotwin_dependentState(ph,en)
case (PLASTIC_DISLOTUNGSTEN_ID) plasticType case (MECHANICAL_PLASTICITY_DISLOTUNGSTEN) plasticType
call dislotungsten_dependentState(ph,en) call dislotungsten_dependentState(ph,en)
case (PLASTIC_NONLOCAL_ID) plasticType case (MECHANICAL_PLASTICITY_NONLOCAL) plasticType
call nonlocal_dependentState(ph,en) call nonlocal_dependentState(ph,en)
end select plasticType end select plasticType
@ -384,19 +384,19 @@ module function plastic_deltaState(ph, en) result(broken)
broken = .false. broken = .false.
select case (phase_plasticity(ph)) select case (mechanical_plasticity_type(ph))
case (PLASTIC_NONLOCAL_ID,PLASTIC_KINEHARDENING_ID) case (MECHANICAL_PLASTICITY_NONLOCAL,MECHANICAL_PLASTICITY_KINEHARDENING)
Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,en)),& Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,en)),&
phase_mechanical_Fi(ph)%data(1:3,1:3,en)),& phase_mechanical_Fi(ph)%data(1:3,1:3,en)),&
phase_mechanical_S(ph)%data(1:3,1:3,en)) phase_mechanical_S(ph)%data(1:3,1:3,en))
plasticType: select case (phase_plasticity(ph)) plasticType: select case (mechanical_plasticity_type(ph))
case (PLASTIC_KINEHARDENING_ID) plasticType case (MECHANICAL_PLASTICITY_KINEHARDENING) plasticType
call plastic_kinehardening_deltaState(Mp,ph,en) call plastic_kinehardening_deltaState(Mp,ph,en)
case (PLASTIC_NONLOCAL_ID) plasticType case (MECHANICAL_PLASTICITY_NONLOCAL) plasticType
call plastic_nonlocal_deltaState(Mp,ph,en) call plastic_nonlocal_deltaState(Mp,ph,en)
end select plasticType end select plasticType

4
src/phase_mechanical_plastic_nonlocal.f90
View File

@ -1252,7 +1252,7 @@ function rhoDotFlux(timestep,ph,en)
!* The entering flux from my neighbor will be distributed on my slip systems according to the !* The entering flux from my neighbor will be distributed on my slip systems according to the
!* compatibility !* compatibility
if (neighbor_n > 0) then if (neighbor_n > 0) then
if (phase_plasticity(np) == PLASTIC_NONLOCAL_ID .and. & if (mechanical_plasticity_type(np) == MECHANICAL_PLASTICITY_NONLOCAL .and. &
any(dependentState(ph)%compatibility(:,:,:,n,en) > 0.0_pREAL)) then any(dependentState(ph)%compatibility(:,:,:,n,en) > 0.0_pREAL)) then
forall (s = 1:ns, t = 1:4) forall (s = 1:ns, t = 1:4)
@ -1298,7 +1298,7 @@ function rhoDotFlux(timestep,ph,en)
!* In case of reduced transmissivity, part of the leaving flux is stored as dead dislocation density. !* In case of reduced transmissivity, part of the leaving flux is stored as dead dislocation density.
!* That means for an interface of zero transmissivity the leaving flux is fully converted to dead dislocations. !* That means for an interface of zero transmissivity the leaving flux is fully converted to dead dislocations.
if (opposite_n > 0) then if (opposite_n > 0) then
if (phase_plasticity(np) == PLASTIC_NONLOCAL_ID) then if (mechanical_plasticity_type(np) == MECHANICAL_PLASTICITY_NONLOCAL) then
normal_me2neighbor_defConf = math_det33(Favg) & normal_me2neighbor_defConf = math_det33(Favg) &
* matmul(math_inv33(transpose(Favg)),geom(ph)%IPareaNormal(1:3,n,en)) ! normal of the interface in (average) deformed configuration (pointing en => neighbor) * matmul(math_inv33(transpose(Favg)),geom(ph)%IPareaNormal(1:3,n,en)) ! normal of the interface in (average) deformed configuration (pointing en => neighbor)

101
src/phase_thermal.f90
View File

@ -15,19 +15,11 @@ submodule(phase) thermal
type(tSourceState), allocatable, dimension(:) :: & type(tSourceState), allocatable, dimension(:) :: &
thermalState thermalState
enum, bind(c); enumerator :: & type :: tFieldQuantities
THERMAL_UNDEFINED_ID ,&
THERMAL_DISSIPATION_ID, &
THERMAL_EXTERNALHEAT_ID
end enum
type :: tDataContainer ! ?? not very telling name. Better: "fieldQuantities" ??
real(pREAL), dimension(:), allocatable :: T, dot_T real(pREAL), dimension(:), allocatable :: T, dot_T
end type tDataContainer end type tFieldQuantities
integer(kind(THERMAL_UNDEFINED_ID)), dimension(:,:), allocatable :: &
thermal_source
type(tDataContainer), dimension(:), allocatable :: current ! ?? not very telling name. Better: "field" ?? MD: current(ho)%T(en) reads quite good type(tFieldQuantities), dimension(:), allocatable :: current
type(tThermalParameters), dimension(:), allocatable :: param type(tThermalParameters), dimension(:), allocatable :: param
@ -36,36 +28,36 @@ submodule(phase) thermal
interface interface
module function dissipation_init(source_length) result(mySources) module function source_dissipation_init(maxNsources) result(isMySource)
integer, intent(in) :: source_length integer, intent(in) :: maxNsources
logical, dimension(:,:), allocatable :: mySources logical, dimension(:,:), allocatable :: isMySource
end function dissipation_init end function source_dissipation_init
module function externalheat_init(source_length) result(mySources) module function source_externalheat_init(maxNsources) result(isMySource)
integer, intent(in) :: source_length integer, intent(in) :: maxNsources
logical, dimension(:,:), allocatable :: mySources logical, dimension(:,:), allocatable :: isMySource
end function externalheat_init end function source_externalheat_init
module subroutine externalheat_dotState(ph, en) module subroutine source_externalheat_dotState(ph, en)
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
en en
end subroutine externalheat_dotState end subroutine source_externalheat_dotState
module function dissipation_f_T(ph,en) result(f_T) module function source_dissipation_f_T(ph,en) result(f_T)
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
en en
real(pREAL) :: f_T real(pREAL) :: f_T
end function dissipation_f_T end function source_dissipation_f_T
module function externalheat_f_T(ph,en) result(f_T) module function source_externalheat_f_T(ph,en) result(f_T)
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
en en
real(pREAL) :: f_T real(pREAL) :: f_T
end function externalheat_f_T end function source_externalheat_f_T
end interface end interface
@ -129,11 +121,11 @@ module subroutine thermal_init(phases)
end do end do
allocate(thermal_source(maxval(thermal_Nsources),phases%length), source = THERMAL_UNDEFINED_ID) allocate(thermal_source_type(maxval(thermal_Nsources),phases%length), source = UNDEFINED)
if (maxval(thermal_Nsources) /= 0) then if (maxval(thermal_Nsources) /= 0) then
where(dissipation_init (maxval(thermal_Nsources))) thermal_source = THERMAL_DISSIPATION_ID where(source_dissipation_init (maxval(thermal_Nsources))) thermal_source_type = THERMAL_SOURCE_DISSIPATION
where(externalheat_init(maxval(thermal_Nsources))) thermal_source = THERMAL_EXTERNALHEAT_ID where(source_externalheat_init(maxval(thermal_Nsources))) thermal_source_type = THERMAL_SOURCE_EXTERNALHEAT
end if end if
thermal_source_maxSizeDotState = 0 thermal_source_maxSizeDotState = 0
@ -151,7 +143,7 @@ end subroutine thermal_init
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief Calculate thermal source. !< @brief Calculate thermal source (forcing term).
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module function phase_f_T(ph,en) result(f) module function phase_f_T(ph,en) result(f)
@ -165,13 +157,13 @@ module function phase_f_T(ph,en) result(f)
f = 0.0_pREAL f = 0.0_pREAL
do so = 1, thermal_Nsources(ph) do so = 1, thermal_Nsources(ph)
select case(thermal_source(so,ph)) select case(thermal_source_type(so,ph))
case (THERMAL_DISSIPATION_ID) case (THERMAL_SOURCE_DISSIPATION)
f = f + dissipation_f_T(ph,en) f = f + source_dissipation_f_T(ph,en)
case (THERMAL_EXTERNALHEAT_ID) case (THERMAL_SOURCE_EXTERNALHEAT)
f = f + externalheat_f_T(ph,en) f = f + source_externalheat_f_T(ph,en)
end select end select
@ -183,22 +175,22 @@ end function phase_f_T
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief tbd. !> @brief tbd.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function phase_thermal_collectDotState(ph,en) result(broken) function phase_thermal_collectDotState(ph,en) result(ok)
integer, intent(in) :: ph, en integer, intent(in) :: ph, en
logical :: broken logical :: ok
integer :: i integer :: i
broken = .false. ok = .true.
SourceLoop: do i = 1, thermal_Nsources(ph) SourceLoop: do i = 1, thermal_Nsources(ph)
if (thermal_source(i,ph) == THERMAL_EXTERNALHEAT_ID) & if (thermal_source_type(i,ph) == THERMAL_SOURCE_EXTERNALHEAT) &
call externalheat_dotState(ph,en) call source_externalheat_dotState(ph,en)
broken = broken .or. any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,en))) ok = ok .and. .not. any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,en)))
end do SourceLoop end do SourceLoop
@ -241,34 +233,35 @@ module function phase_thermal_constitutive(Delta_t,ph,en) result(converged_)
logical :: converged_ logical :: converged_
converged_ = .not. integrateThermalState(Delta_t,ph,en) converged_ = integrateThermalState(Delta_t,ph,en)
end function phase_thermal_constitutive end function phase_thermal_constitutive
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief integrate state with 1st order explicit Euler method !> @brief Integrate state with 1st order explicit Euler method.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function integrateThermalState(Delta_t, ph,en) result(broken) function integrateThermalState(Delta_t, ph,en) result(converged)
real(pREAL), intent(in) :: Delta_t real(pREAL), intent(in) :: Delta_t
integer, intent(in) :: ph, en integer, intent(in) :: ph, en
logical :: & logical :: converged
broken
integer :: & integer :: &
so, & so, &
sizeDotState sizeDotState
broken = phase_thermal_collectDotState(ph,en) converged = phase_thermal_collectDotState(ph,en)
if (broken) return if (converged) then
do so = 1, thermal_Nsources(ph) do so = 1, thermal_Nsources(ph)
sizeDotState = thermalState(ph)%p(so)%sizeDotState sizeDotState = thermalState(ph)%p(so)%sizeDotState
thermalState(ph)%p(so)%state(1:sizeDotState,en) = thermalState(ph)%p(so)%state0(1:sizeDotState,en) & thermalState(ph)%p(so)%state(1:sizeDotState,en) = thermalState(ph)%p(so)%state0(1:sizeDotState,en) &
+ thermalState(ph)%p(so)%dotState(1:sizeDotState,en) * Delta_t + thermalState(ph)%p(so)%dotState(1:sizeDotState,en) * Delta_t
end do end do
end if
end function integrateThermalState end function integrateThermalState
@ -318,7 +311,7 @@ end subroutine thermal_forward
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief Get temperature (for use by non-thermal physics) !< @brief Get temperature (for use by non-thermal physics).
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
pure module function thermal_T(ph,en) result(T) pure module function thermal_T(ph,en) result(T)
@ -332,7 +325,7 @@ end function thermal_T
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief Get rate of temperature (for use by non-thermal physics) !< @brief Get rate of temperature (for use by non-thermal physics).
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module function thermal_dot_T(ph,en) result(dot_T) module function thermal_dot_T(ph,en) result(dot_T)

36
src/phase_thermal_dissipation.f90 → src/phase_thermal_source_dissipation.f90
View File

@ -5,7 +5,7 @@
!> @brief material subroutine for thermal source due to plastic dissipation !> @brief material subroutine for thermal source due to plastic dissipation
!> @details to be done !> @details to be done
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
submodule(phase:thermal) dissipation submodule(phase:thermal) source_dissipation
type :: tParameters !< container type for internal constitutive parameters type :: tParameters !< container type for internal constitutive parameters
real(pREAL) :: & real(pREAL) :: &
@ -22,10 +22,10 @@ contains
!> @brief module initialization !> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function dissipation_init(source_length) result(mySources) module function source_dissipation_init(maxNsources) result(isMySource)
integer, intent(in) :: source_length integer, intent(in) :: maxNsources
logical, dimension(:,:), allocatable :: mySources logical, dimension(:,:), allocatable :: isMySource
type(tDict), pointer :: & type(tDict), pointer :: &
phases, & phases, &
@ -35,26 +35,29 @@ module function dissipation_init(source_length) result(mySources)
class(tList), pointer :: & class(tList), pointer :: &
sources sources
character(len=:), allocatable :: refs character(len=:), allocatable :: refs
integer :: so,Nmembers,ph integer :: ph,Nmembers,so,Nsources
mySources = thermal_active('dissipation',source_length) isMySource = thermal_active('dissipation',maxNsources)
if (count(mySources) == 0) return if (count(isMySource) == 0) return
print'(/,1x,a)', '<<<+- phase:thermal:dissipation init -+>>>' print'(/,1x,a)', '<<<+- phase:thermal:source_dissipation init -+>>>'
print'(/,a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT) print'(/,a,i2)', ' # phases: ',count(isMySource); flush(IO_STDOUT)
phases => config_material%get_dict('phase') phases => config_material%get_dict('phase')
allocate(param(phases%length)) allocate(param(phases%length))
do ph = 1, phases%length do ph = 1, phases%length
Nsources = count(isMySource(:,ph))
if (Nsources == 0) cycle
if (Nsources > 1) call IO_error(600,ext_msg='dissipation')
Nmembers = count(material_ID_phase == ph)
phase => phases%get_dict(ph) phase => phases%get_dict(ph)
if (count(mySources(:,ph)) == 0) cycle !ToDo: error if > 1
thermal => phase%get_dict('thermal') thermal => phase%get_dict('thermal')
sources => thermal%get_list('source') sources => thermal%get_list('source')
do so = 1, sources%length do so = 1, sources%length
if (mySources(so,ph)) then if (isMySource(so,ph)) then
associate(prm => param(ph)) associate(prm => param(ph))
src => sources%get_dict(so) src => sources%get_dict(so)
print'(1x,a,i0,a,i0)', 'phase ',ph,' source ',so print'(1x,a,i0,a,i0)', 'phase ',ph,' source ',so
@ -62,22 +65,21 @@ module function dissipation_init(source_length) result(mySources)
if (len(refs) > 0) print'(/,1x,a)', refs if (len(refs) > 0) print'(/,1x,a)', refs
prm%kappa = src%get_asReal('kappa') prm%kappa = src%get_asReal('kappa')
Nmembers = count(material_ID_phase == ph)
call phase_allocateState(thermalState(ph)%p(so),Nmembers,0,0,0) call phase_allocateState(thermalState(ph)%p(so),Nmembers,0,0,0)
end associate end associate
exit
end if end if
end do end do
end do end do
end function dissipation_init end function source_dissipation_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Ninstancess dissipation rate !> @brief Ninstancess dissipation rate
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function dissipation_f_T(ph,en) result(f_T) module function source_dissipation_f_T(ph,en) result(f_T)
integer, intent(in) :: ph, en integer, intent(in) :: ph, en
real(pREAL) :: & real(pREAL) :: &
@ -91,6 +93,6 @@ module function dissipation_f_T(ph,en) result(f_T)
f_T = prm%kappa*sum(abs(Mp*mechanical_L_p(ph,en))) f_T = prm%kappa*sum(abs(Mp*mechanical_L_p(ph,en)))
end associate end associate
end function dissipation_f_T end function source_dissipation_f_T
end submodule dissipation end submodule source_dissipation

61
src/phase_thermal_externalheat.f90 → src/phase_thermal_source_externalheat.f90
View File

@ -4,14 +4,14 @@
!> @author Philip Eisenlohr, Michigan State University !> @author Philip Eisenlohr, Michigan State University
!> @brief material subroutine for variable heat source !> @brief material subroutine for variable heat source
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
submodule(phase:thermal) externalheat submodule(phase:thermal) source_externalheat
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
source_thermal_externalheat_offset !< which source is my current thermal dissipation mechanism? source_ID !< index in phase source list corresponding to this source
type :: tParameters !< container type for internal constitutive parameters type :: tParameters !< container type for internal constitutive parameters
type(tTable) :: f type(tTable) :: f !< external heat power as (tabulated) function of time
end type tParameters end type tParameters
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances)
@ -24,10 +24,10 @@ contains
!> @brief module initialization !> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function externalheat_init(source_length) result(mySources) module function source_externalheat_init(maxNsources) result(isMySource)
integer, intent(in) :: source_length integer, intent(in) :: maxNsources
logical, dimension(:,:), allocatable :: mySources logical, dimension(:,:), allocatable :: isMySource
type(tDict), pointer :: & type(tDict), pointer :: &
phases, & phases, &
@ -37,28 +37,31 @@ module function externalheat_init(source_length) result(mySources)
type(tList), pointer :: & type(tList), pointer :: &
sources sources
character(len=:), allocatable :: refs character(len=:), allocatable :: refs
integer :: so,Nmembers,ph integer :: ph,Nmembers,so,Nsources
mySources = thermal_active('externalheat',source_length) isMySource = thermal_active('externalheat',maxNsources)
if (count(mySources) == 0) return if (count(isMySource) == 0) return
print'(/,1x,a)', '<<<+- phase:thermal:externalheat init -+>>>' print'(/,1x,a)', '<<<+- phase:thermal:source_externalheat init -+>>>'
print'(/,a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT) print'(/,a,i2)', ' # phases: ',count(isMySource); flush(IO_STDOUT)
phases => config_material%get_dict('phase') phases => config_material%get_dict('phase')
allocate(param(phases%length)) allocate(param(phases%length))
allocate(source_thermal_externalheat_offset (phases%length), source=0) allocate(source_ID(phases%length), source=0)
do ph = 1, phases%length do ph = 1, phases%length
Nsources = count(isMySource(:,ph))
if (Nsources == 0) cycle
if (Nsources > 1) call IO_error(600,ext_msg='externalheat')
Nmembers = count(material_ID_phase == ph)
phase => phases%get_dict(ph) phase => phases%get_dict(ph)
if (count(mySources(:,ph)) == 0) cycle
thermal => phase%get_dict('thermal') thermal => phase%get_dict('thermal')
sources => thermal%get_list('source') sources => thermal%get_list('source')
do so = 1, sources%length do so = 1, sources%length
if (mySources(so,ph)) then if (isMySource(so,ph)) then
source_thermal_externalheat_offset(ph) = so source_ID(ph) = so
associate(prm => param(ph)) associate(prm => param(ph))
src => sources%get_dict(so) src => sources%get_dict(so)
print'(1x,a,i0,a,i0)', 'phase ',ph,' source ',so print'(1x,a,i0,a,i0)', 'phase ',ph,' source ',so
@ -66,41 +69,36 @@ module function externalheat_init(source_length) result(mySources)
if (len(refs) > 0) print'(/,1x,a)', refs if (len(refs) > 0) print'(/,1x,a)', refs
prm%f = table(src,'t','f') prm%f = table(src,'t','f')
Nmembers = count(material_ID_phase == ph)
call phase_allocateState(thermalState(ph)%p(so),Nmembers,1,1,0) call phase_allocateState(thermalState(ph)%p(so),Nmembers,1,1,0)
end associate end associate
exit
end if end if
end do end do
end do end do
end function externalheat_init end function source_externalheat_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief rate of change of state !> @brief rate of change of state
!> @details state only contains current time to linearly interpolate given heat powers !> @details state only contains current time to linearly interpolate given heat powers
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine externalheat_dotState(ph, en) module subroutine source_externalheat_dotState(ph, en)
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
en en
integer :: &
so
so = source_thermal_externalheat_offset(ph) thermalState(ph)%p(source_ID(ph))%dotState(1,en) = 1.0_pREAL ! state is current time
thermalState(ph)%p(so)%dotState(1,en) = 1.0_pREAL ! state is current time end subroutine source_externalheat_dotState
end subroutine externalheat_dotState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns local heat generation rate !> @brief returns local heat generation rate
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function externalheat_f_T(ph,en) result(f_T) module function source_externalheat_f_T(ph,en) result(f_T)
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
@ -108,16 +106,11 @@ module function externalheat_f_T(ph,en) result(f_T)
real(pREAL) :: & real(pREAL) :: &
f_T f_T
integer :: &
so
so = source_thermal_externalheat_offset(ph)
associate(prm => param(ph)) associate(prm => param(ph))
f_T = prm%f%at(thermalState(ph)%p(so)%state(1,en)) f_T = prm%f%at(thermalState(ph)%p(source_ID(ph))%state(1,en))
end associate end associate
end function externalheat_f_T end function source_externalheat_f_T
end submodule externalheat end submodule source_externalheat