diff --git a/PRIVATE b/PRIVATE index e74cf0062..2ed5cd4ba 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit e74cf00628285a587ced1e551cc9837c1011ca1c +Subproject commit 2ed5cd4ba97b44ad9c8b61ced94060aee57a2dd8 diff --git a/examples/SpectralMethod/Polycrystal/20grains.seeds b/examples/grid/20grains.seeds similarity index 100% rename from examples/SpectralMethod/Polycrystal/20grains.seeds rename to examples/grid/20grains.seeds diff --git a/examples/SpectralMethod/Polycrystal/20grains16x16x16.vtr b/examples/grid/20grains16x16x16.vtr similarity index 100% rename from examples/SpectralMethod/Polycrystal/20grains16x16x16.vtr rename to examples/grid/20grains16x16x16.vtr diff --git a/examples/SpectralMethod/Polycrystal/20grains32x32x32.vtr b/examples/grid/20grains32x32x32.vtr similarity index 100% rename from examples/SpectralMethod/Polycrystal/20grains32x32x32.vtr rename to examples/grid/20grains32x32x32.vtr diff --git a/examples/SpectralMethod/Polycrystal/20grains64x64x64.vtr b/examples/grid/20grains64x64x64.vtr similarity index 100% rename from examples/SpectralMethod/Polycrystal/20grains64x64x64.vtr rename to examples/grid/20grains64x64x64.vtr diff --git a/examples/SpectralMethod/Polycrystal/material.yaml b/examples/grid/material.yaml similarity index 100% rename from examples/SpectralMethod/Polycrystal/material.yaml rename to examples/grid/material.yaml diff --git a/examples/SpectralMethod/Polycrystal/numerics.yaml b/examples/grid/numerics.yaml similarity index 100% rename from examples/SpectralMethod/Polycrystal/numerics.yaml rename to examples/grid/numerics.yaml diff --git a/examples/SpectralMethod/Polycrystal/shearXY.yaml b/examples/grid/shearXY.yaml similarity index 100% rename from examples/SpectralMethod/Polycrystal/shearXY.yaml rename to examples/grid/shearXY.yaml diff --git a/examples/SpectralMethod/Polycrystal/shearZX.yaml b/examples/grid/shearZX.yaml similarity index 100% rename from examples/SpectralMethod/Polycrystal/shearZX.yaml rename to examples/grid/shearZX.yaml diff --git a/examples/SpectralMethod/Polycrystal/tensionX.yaml b/examples/grid/tensionX.yaml similarity index 100% rename from examples/SpectralMethod/Polycrystal/tensionX.yaml rename to examples/grid/tensionX.yaml diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index a191298a3..97e7520f9 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -41,8 +41,6 @@ #include "phase_damage_isoductile.f90" #include "phase_damage_anisobrittle.f90" #include "phase_damage_anisoductile.f90" -#include "damage_none.f90" -#include "damage_nonlocal.f90" #include "homogenization.f90" #include "homogenization_mechanical.f90" #include "homogenization_mechanical_pass.f90" diff --git a/src/config.f90 b/src/config.f90 index b10edf013..02b16f2a2 100644 --- a/src/config.f90 +++ b/src/config.f90 @@ -5,16 +5,10 @@ !! precedence over material.yaml. !-------------------------------------------------------------------------------------------------- module config - use prec - use DAMASK_interface use IO use YAML_parse use YAML_types -#ifdef PETSc -#include - use petscsys -#endif implicit none private @@ -50,17 +44,12 @@ end subroutine config_init subroutine parse_material logical :: fileExists - character(len=:), allocatable :: fname - fname = getSolverJobName()//'.yaml' - inquire(file=fname,exist=fileExists) - if(.not. fileExists) then - fname = 'material.yaml' - inquire(file=fname,exist=fileExists) - if(.not. fileExists) call IO_error(100,ext_msg=fname) - endif - print*, 'reading '//fname; flush(IO_STDOUT) - config_material => YAML_parse_file(fname) + + inquire(file='material.yaml',exist=fileExists) + if(.not. fileExists) call IO_error(100,ext_msg='material.yaml') + print*, 'reading material.yaml'; flush(IO_STDOUT) + config_material => YAML_parse_file('material.yaml') end subroutine parse_material @@ -72,6 +61,7 @@ subroutine parse_numerics logical :: fexist + config_numerics => emptyDict inquire(file='numerics.yaml', exist=fexist) if (fexist) then @@ -89,6 +79,7 @@ subroutine parse_debug logical :: fexist + config_debug => emptyDict inquire(file='debug.yaml', exist=fexist) fileExists: if (fexist) then diff --git a/src/damage_none.f90 b/src/damage_none.f90 deleted file mode 100644 index 680110d55..000000000 --- a/src/damage_none.f90 +++ /dev/null @@ -1,38 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine for constant damage field -!-------------------------------------------------------------------------------------------------- -module damage_none - use prec - use config - use material - - implicit none - public - -contains - -!-------------------------------------------------------------------------------------------------- -!> @brief allocates all neccessary fields, reads information from material configuration file -!-------------------------------------------------------------------------------------------------- -subroutine damage_none_init - - integer :: h,Nmaterialpoints - - print'(/,a)', ' <<<+- damage_none init -+>>>'; flush(6) - - do h = 1, size(material_name_homogenization) - if (damage_type(h) /= DAMAGE_NONE_ID) cycle - - Nmaterialpoints = count(material_homogenizationAt == h) - damageState_h(h)%sizeState = 0 - allocate(damageState_h(h)%state0 (0,Nmaterialpoints)) - allocate(damageState_h(h)%state (0,Nmaterialpoints)) - - allocate (damage(h)%p(Nmaterialpoints), source=1.0_pReal) - - enddo - -end subroutine damage_none_init - -end module damage_none diff --git a/src/damage_nonlocal.f90 b/src/damage_nonlocal.f90 deleted file mode 100644 index 807231889..000000000 --- a/src/damage_nonlocal.f90 +++ /dev/null @@ -1,94 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine for non-locally evolving damage field -!-------------------------------------------------------------------------------------------------- -module damage_nonlocal - use prec - use material - use config - use YAML_types - use lattice - use phase - use results - - implicit none - private - - type, private :: tNumerics - real(pReal) :: & - charLength !< characteristic length scale for gradient problems - end type tNumerics - - type(tNumerics), private :: & - num - - public :: & - damage_nonlocal_init, & - damage_nonlocal_getDiffusion - -contains - -!-------------------------------------------------------------------------------------------------- -!> @brief module initialization -!> @details reads in material parameters, allocates arrays, and does sanity checks -!-------------------------------------------------------------------------------------------------- -subroutine damage_nonlocal_init - - integer :: Ninstances,Nmaterialpoints,h - class(tNode), pointer :: & - num_generic, & - material_homogenization - - print'(/,a)', ' <<<+- damage_nonlocal init -+>>>'; flush(6) - -!------------------------------------------------------------------------------------ -! read numerics parameter - num_generic => config_numerics%get('generic',defaultVal= emptyDict) - num%charLength = num_generic%get_asFloat('charLength',defaultVal=1.0_pReal) - - Ninstances = count(damage_type == DAMAGE_nonlocal_ID) - - material_homogenization => config_material%get('homogenization') - do h = 1, material_homogenization%length - if (damage_type(h) /= DAMAGE_NONLOCAL_ID) cycle - - Nmaterialpoints = count(material_homogenizationAt == h) - damageState_h(h)%sizeState = 1 - allocate(damageState_h(h)%state0 (1,Nmaterialpoints), source=1.0_pReal) - allocate(damageState_h(h)%state (1,Nmaterialpoints), source=1.0_pReal) - - damage(h)%p => damageState_h(h)%state(1,:) - - enddo - -end subroutine damage_nonlocal_init - - -!-------------------------------------------------------------------------------------------------- -!> @brief returns homogenized non local damage diffusion tensor in reference configuration -!-------------------------------------------------------------------------------------------------- -function damage_nonlocal_getDiffusion(ip,el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), dimension(3,3) :: & - damage_nonlocal_getDiffusion - integer :: & - homog, & - grain - - homog = material_homogenizationAt(el) - damage_nonlocal_getDiffusion = 0.0_pReal - do grain = 1, homogenization_Nconstituents(homog) - damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + & - crystallite_push33ToRef(grain,ip,el,lattice_D(1:3,1:3,material_phaseAt(grain,el))) - enddo - - damage_nonlocal_getDiffusion = & - num%charLength**2*damage_nonlocal_getDiffusion/real(homogenization_Nconstituents(homog),pReal) - -end function damage_nonlocal_getDiffusion - - -end module damage_nonlocal diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index dc78d9e44..afbae027f 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -15,7 +15,6 @@ module grid_damage_spectral use IO use spectral_utilities use discretization_grid - use damage_nonlocal use YAML_types use homogenization use config diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 62ac52f06..3e1b939b5 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -12,14 +12,53 @@ module homogenization use material use phase use discretization - use damage_none - use damage_nonlocal use HDF5_utilities use results + use lattice implicit none private + + enum, bind(c); enumerator :: & + THERMAL_ISOTHERMAL_ID, & + THERMAL_CONDUCTION_ID, & + DAMAGE_NONE_ID, & + DAMAGE_NONLOCAL_ID, & + HOMOGENIZATION_UNDEFINED_ID, & + HOMOGENIZATION_NONE_ID, & + HOMOGENIZATION_ISOSTRAIN_ID, & + HOMOGENIZATION_RGC_ID + end enum + + type(tState), allocatable, dimension(:), public :: & + homogState, & + damageState_h + + integer, dimension(:), allocatable, public, protected :: & + homogenization_typeInstance, & !< instance of particular type of each homogenization + thermal_typeInstance, & !< instance of particular type of each thermal transport + damage_typeInstance !< instance of particular type of each nonlocal damage + + real(pReal), dimension(:), allocatable, public, protected :: & + thermal_initialT + + integer(kind(THERMAL_isothermal_ID)), dimension(:), allocatable, public, protected :: & + thermal_type !< thermal transport model + integer(kind(DAMAGE_none_ID)), dimension(:), allocatable, public, protected :: & + damage_type !< nonlocal damage model + integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:), allocatable, public, protected :: & + homogenization_type !< type of each homogenization + + type, private :: tNumerics_damage + real(pReal) :: & + charLength !< characteristic length scale for gradient problems + end type tNumerics_damage + + type(tNumerics_damage), private :: & + num_damage + + logical, public :: & terminallyIll = .false. !< at least one material point is terminally ill @@ -193,7 +232,13 @@ module homogenization homogenization_forward, & homogenization_results, & homogenization_restartRead, & - homogenization_restartWrite + homogenization_restartWrite, & + THERMAL_CONDUCTION_ID, & + DAMAGE_NONLOCAL_ID + + public :: & + damage_nonlocal_init, & + damage_nonlocal_getDiffusion contains @@ -209,6 +254,12 @@ subroutine homogenization_init() print'(/,a)', ' <<<+- homogenization init -+>>>'; flush(IO_STDOUT) + + allocate(homogState (size(material_name_homogenization))) + allocate(damageState_h (size(material_name_homogenization))) + call material_parseHomogenization + + num_homog => config_numerics%get('homogenization',defaultVal=emptyDict) num_homogGeneric => num_homog%get('generic',defaultVal=emptyDict) @@ -220,8 +271,7 @@ subroutine homogenization_init() call thermal_init() call damage_init() - if (any(damage_type == DAMAGE_none_ID)) call damage_none_init - if (any(damage_type == DAMAGE_nonlocal_ID)) call damage_nonlocal_init + call damage_nonlocal_init end subroutine homogenization_init @@ -257,6 +307,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE if(homogState(ho)%sizeState > 0) homogState(ho)%State(:,me) = homogState(ho)%State0(:,me) if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%State(:,me) = damageState_h(ho)%State0(:,me) + call damage_partition(ce) doneAndHappy = [.false.,.true.] @@ -311,26 +362,6 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE enddo !$OMP END DO -! !$OMP DO PRIVATE(ho,ph,ce) -! do el = FEsolving_execElem(1),FEsolving_execElem(2) -! if (terminallyIll) continue -! ho = material_homogenizationAt(el) -! do ip = FEsolving_execIP(1),FEsolving_execIP(2) -! ce = (el-1)*discretization_nIPs + ip -! call damage_partition(ce) -! do co = 1, homogenization_Nconstituents(ho) -! ph = material_phaseAt(co,el) -! if (.not. thermal_stress(dt,ph,material_phaseMemberAt(co,ip,el))) then -! if (.not. terminallyIll) & ! so first signals terminally ill... -! print*, ' Integration point ', ip,' at element ', el, ' terminally ill' -! terminallyIll = .true. ! ...and kills all others -! endif -! call thermal_homogenize(ip,el) -! enddo -! enddo -! enddo -! !$OMP END DO - !$OMP DO PRIVATE(ho) elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2) ho = material_homogenizationAt(el) @@ -397,7 +428,8 @@ subroutine homogenization_forward do ho = 1, size(material_name_homogenization) homogState (ho)%state0 = homogState (ho)%state - damageState_h(ho)%state0 = damageState_h(ho)%state + if(damageState_h(ho)%sizeState > 0) & + damageState_h(ho)%state0 = damageState_h(ho)%state enddo end subroutine homogenization_forward @@ -457,4 +489,143 @@ subroutine homogenization_restartRead(fileHandle) end subroutine homogenization_restartRead + +!-------------------------------------------------------------------------------------------------- +!> @brief module initialization +!> @details reads in material parameters, allocates arrays, and does sanity checks +!-------------------------------------------------------------------------------------------------- +subroutine damage_nonlocal_init + + integer :: Ninstances,Nmaterialpoints,h + class(tNode), pointer :: & + num_generic, & + material_homogenization + + print'(/,a)', ' <<<+- damage_nonlocal init -+>>>'; flush(6) + +!------------------------------------------------------------------------------------ +! read numerics parameter + num_generic => config_numerics%get('generic',defaultVal= emptyDict) + num_damage%charLength = num_generic%get_asFloat('charLength',defaultVal=1.0_pReal) + + Ninstances = count(damage_type == DAMAGE_nonlocal_ID) + + material_homogenization => config_material%get('homogenization') + do h = 1, material_homogenization%length + if (damage_type(h) /= DAMAGE_NONLOCAL_ID) cycle + + Nmaterialpoints = count(material_homogenizationAt == h) + damageState_h(h)%sizeState = 1 + allocate(damageState_h(h)%state0 (1,Nmaterialpoints), source=1.0_pReal) + allocate(damageState_h(h)%state (1,Nmaterialpoints), source=1.0_pReal) + + enddo + +end subroutine damage_nonlocal_init + + +!-------------------------------------------------------------------------------------------------- +!> @brief returns homogenized non local damage diffusion tensor in reference configuration +!-------------------------------------------------------------------------------------------------- +function damage_nonlocal_getDiffusion(ip,el) + + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), dimension(3,3) :: & + damage_nonlocal_getDiffusion + integer :: & + homog, & + grain + + homog = material_homogenizationAt(el) + damage_nonlocal_getDiffusion = 0.0_pReal + do grain = 1, homogenization_Nconstituents(homog) + damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + & + crystallite_push33ToRef(grain,ip,el,lattice_D(1:3,1:3,material_phaseAt(grain,el))) + enddo + + damage_nonlocal_getDiffusion = & + num_damage%charLength**2*damage_nonlocal_getDiffusion/real(homogenization_Nconstituents(homog),pReal) + +end function damage_nonlocal_getDiffusion + + +!-------------------------------------------------------------------------------------------------- +!> @brief parses the homogenization part from the material configuration +! ToDo: This should be done in homogenization +!-------------------------------------------------------------------------------------------------- +subroutine material_parseHomogenization + + class(tNode), pointer :: & + material_homogenization, & + homog, & + homogMech, & + homogThermal, & + homogDamage + + integer :: h + + material_homogenization => config_material%get('homogenization') + + allocate(homogenization_type(size(material_name_homogenization)), source=HOMOGENIZATION_undefined_ID) + allocate(thermal_type(size(material_name_homogenization)), source=THERMAL_isothermal_ID) + allocate(damage_type (size(material_name_homogenization)), source=DAMAGE_none_ID) + allocate(homogenization_typeInstance(size(material_name_homogenization)), source=0) + allocate(thermal_typeInstance(size(material_name_homogenization)), source=0) + allocate(damage_typeInstance(size(material_name_homogenization)), source=0) + allocate(thermal_initialT(size(material_name_homogenization)), source=300.0_pReal) + + do h=1, size(material_name_homogenization) + homog => material_homogenization%get(h) + homogMech => homog%get('mechanics') + select case (homogMech%get_asString('type')) + case('pass') + homogenization_type(h) = HOMOGENIZATION_NONE_ID + case('isostrain') + homogenization_type(h) = HOMOGENIZATION_ISOSTRAIN_ID + case('RGC') + homogenization_type(h) = HOMOGENIZATION_RGC_ID + case default + call IO_error(500,ext_msg=homogMech%get_asString('type')) + end select + + homogenization_typeInstance(h) = count(homogenization_type==homogenization_type(h)) + + if(homog%contains('thermal')) then + homogThermal => homog%get('thermal') + thermal_initialT(h) = homogThermal%get_asFloat('T_0',defaultVal=300.0_pReal) + + select case (homogThermal%get_asString('type')) + case('isothermal') + thermal_type(h) = THERMAL_isothermal_ID + case('conduction') + thermal_type(h) = THERMAL_conduction_ID + case default + call IO_error(500,ext_msg=homogThermal%get_asString('type')) + end select + endif + + if(homog%contains('damage')) then + homogDamage => homog%get('damage') + select case (homogDamage%get_asString('type')) + case('none') + damage_type(h) = DAMAGE_none_ID + case('nonlocal') + damage_type(h) = DAMAGE_nonlocal_ID + case default + call IO_error(500,ext_msg=homogDamage%get_asString('type')) + end select + endif + enddo + + do h=1, size(material_name_homogenization) + homogenization_typeInstance(h) = count(homogenization_type(1:h) == homogenization_type(h)) + thermal_typeInstance(h) = count(thermal_type (1:h) == thermal_type (h)) + damage_typeInstance(h) = count(damage_type (1:h) == damage_type (h)) + enddo + +end subroutine material_parseHomogenization + + end module homogenization diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90 index 8b73ed54a..502174309 100644 --- a/src/homogenization_damage.f90 +++ b/src/homogenization_damage.f90 @@ -72,7 +72,8 @@ module subroutine damage_partition(ce) integer :: co - phi = current(material_homogenizationAt2(ce))%phi(material_homogenizationMemberAt2(ce)) + if(damageState_h(material_homogenizationAt2(ce))%sizeState < 1) return + phi = damagestate_h(material_homogenizationAt2(ce))%state(1,material_homogenizationMemberAt2(ce)) do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) call phase_damage_set_phi(phi,co,ce) enddo @@ -143,7 +144,7 @@ module subroutine damage_nonlocal_putNonLocalDamage(phi,ip,el) homog = material_homogenizationAt(el) offset = material_homogenizationMemberAt(ip,el) - damage(homog)%p(offset) = phi + damagestate_h(homog)%state(1,offset) = phi end subroutine damage_nonlocal_putNonLocalDamage @@ -162,7 +163,7 @@ module subroutine damage_nonlocal_results(homog,group) outputsLoop: do o = 1,size(prm%output) select case(prm%output(o)) case ('phi') - call results_writeDataset(group,damage(homog)%p,prm%output(o),& + call results_writeDataset(group,damagestate_h(homog)%state(1,:),prm%output(o),& 'damage indicator','-') end select enddo outputsLoop diff --git a/src/homogenization_mechanical.f90 b/src/homogenization_mechanical.f90 index 91e1d3d84..a19a61f72 100644 --- a/src/homogenization_mechanical.f90 +++ b/src/homogenization_mechanical.f90 @@ -2,7 +2,7 @@ !> @author Martin Diehl, KU Leuven !> @brief Partition F and homogenize P/dPdF !-------------------------------------------------------------------------------------------------- -submodule(homogenization) mechanics +submodule(homogenization) mechanical interface @@ -86,7 +86,7 @@ module subroutine mechanical_init(num_homog) class(tNode), pointer :: & num_homogMech - print'(/,a)', ' <<<+- homogenization:mechanics init -+>>>' + print'(/,a)', ' <<<+- homogenization:mechanical init -+>>>' allocate(homogenization_dPdF(3,3,3,3,discretization_nIPs*discretization_Nelems), source=0.0_pReal) homogenization_F0 = spread(math_I3,3,discretization_nIPs*discretization_Nelems) ! initialize to identity @@ -225,8 +225,6 @@ end function mechanical_updateState !> @brief Write results to file. !-------------------------------------------------------------------------------------------------- module subroutine mechanical_results(group_base,h) - use material, only: & - material_homogenization_type => homogenization_type character(len=*), intent(in) :: group_base integer, intent(in) :: h @@ -236,7 +234,7 @@ module subroutine mechanical_results(group_base,h) group = trim(group_base)//'/mech' call results_closeGroup(results_addGroup(group)) - select case(material_homogenization_type(h)) + select case(homogenization_type(h)) case(HOMOGENIZATION_rgc_ID) call mechanical_RGC_results(homogenization_typeInstance(h),group) @@ -253,4 +251,4 @@ module subroutine mechanical_results(group_base,h) end subroutine mechanical_results -end submodule mechanics +end submodule mechanical diff --git a/src/homogenization_mechanical_RGC.f90 b/src/homogenization_mechanical_RGC.f90 index 546c2b65c..f1bcda380 100644 --- a/src/homogenization_mechanical_RGC.f90 +++ b/src/homogenization_mechanical_RGC.f90 @@ -6,7 +6,7 @@ !> @brief Relaxed grain cluster (RGC) homogenization scheme !> N_constituents is defined as p x q x r (cluster) !-------------------------------------------------------------------------------------------------- -submodule(homogenization:mechanics) RGC +submodule(homogenization:mechanical) RGC use rotations use lattice @@ -88,7 +88,7 @@ module subroutine mechanical_RGC_init(num_homogMech) homog, & homogMech - print'(/,a)', ' <<<+- homogenization:mechanics:RGC init -+>>>' + print'(/,a)', ' <<<+- homogenization:mechanical:RGC init -+>>>' Ninstances = count(homogenization_type == HOMOGENIZATION_RGC_ID) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) diff --git a/src/homogenization_mechanical_isostrain.f90 b/src/homogenization_mechanical_isostrain.f90 index f9616de42..b492499d8 100644 --- a/src/homogenization_mechanical_isostrain.f90 +++ b/src/homogenization_mechanical_isostrain.f90 @@ -4,7 +4,7 @@ !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @brief Isostrain (full constraint Taylor assuption) homogenization scheme !-------------------------------------------------------------------------------------------------- -submodule(homogenization:mechanics) isostrain +submodule(homogenization:mechanical) isostrain enum, bind(c); enumerator :: & parallel_ID, & @@ -37,7 +37,7 @@ module subroutine mechanical_isostrain_init homog, & homogMech - print'(/,a)', ' <<<+- homogenization:mechanics:isostrain init -+>>>' + print'(/,a)', ' <<<+- homogenization:mechanical:isostrain init -+>>>' Ninstances = count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) diff --git a/src/homogenization_mechanical_pass.f90 b/src/homogenization_mechanical_pass.f90 index 9e8f3e44c..6217e6836 100644 --- a/src/homogenization_mechanical_pass.f90 +++ b/src/homogenization_mechanical_pass.f90 @@ -4,7 +4,7 @@ !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief dummy homogenization homogenization scheme for 1 constituent per material point !-------------------------------------------------------------------------------------------------- -submodule(homogenization:mechanics) none +submodule(homogenization:mechanical) none contains @@ -18,7 +18,7 @@ module subroutine mechanical_pass_init h, & Nmaterialpoints - print'(/,a)', ' <<<+- homogenization:mechanics:pass init -+>>>' + print'(/,a)', ' <<<+- homogenization:mechanical:pass init -+>>>' Ninstances = count(homogenization_type == HOMOGENIZATION_NONE_ID) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) diff --git a/src/lattice.f90 b/src/lattice.f90 index c9b6d99ef..e5f5453d4 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -459,7 +459,8 @@ subroutine lattice_init phase, & mech, & elasticity, & - thermal + thermal, & + damage print'(/,a)', ' <<<+- lattice init -+>>>'; flush(IO_STDOUT) @@ -535,13 +536,17 @@ subroutine lattice_init endif - lattice_D(1,1,ph) = phase%get_asFloat('D_11',defaultVal=0.0_pReal) - lattice_D(2,2,ph) = phase%get_asFloat('D_22',defaultVal=0.0_pReal) - lattice_D(3,3,ph) = phase%get_asFloat('D_33',defaultVal=0.0_pReal) - lattice_D(1:3,1:3,ph) = lattice_applyLatticeSymmetry33(lattice_D(1:3,1:3,ph), & + if (phase%contains('damage')) then + damage => phase%get('damage') + damage => damage%get(1) + lattice_D(1,1,ph) = damage%get_asFloat('D_11',defaultVal=0.0_pReal) + lattice_D(2,2,ph) = damage%get_asFloat('D_22',defaultVal=0.0_pReal) + lattice_D(3,3,ph) = damage%get_asFloat('D_33',defaultVal=0.0_pReal) + lattice_D(1:3,1:3,ph) = lattice_applyLatticeSymmetry33(lattice_D(1:3,1:3,ph), & phase%get_asString('lattice')) - lattice_M(ph) = phase%get_asFloat('M',defaultVal=0.0_pReal) + lattice_M(ph) = damage%get_asFloat('M',defaultVal=0.0_pReal) + endif ! SHOULD NOT BE PART OF LATTICE END call selfTest diff --git a/src/material.f90 b/src/material.f90 index ccc6898be..3656f0f0a 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -6,50 +6,26 @@ !-------------------------------------------------------------------------------------------------- module material use prec - use math use config use results use IO use rotations use discretization + use YAML_types implicit none private - enum, bind(c); enumerator :: & - THERMAL_ISOTHERMAL_ID, & - THERMAL_CONDUCTION_ID, & - DAMAGE_NONE_ID, & - DAMAGE_NONLOCAL_ID, & - HOMOGENIZATION_UNDEFINED_ID, & - HOMOGENIZATION_NONE_ID, & - HOMOGENIZATION_ISOSTRAIN_ID, & - HOMOGENIZATION_RGC_ID - end enum + integer, dimension(:), allocatable, public, protected :: & + homogenization_Nconstituents !< number of grains in each homogenization character(len=:), public, protected, allocatable, dimension(:) :: & material_name_phase, & !< name of each phase material_name_homogenization !< name of each homogenization - integer(kind(THERMAL_isothermal_ID)), dimension(:), allocatable, public, protected :: & - thermal_type !< thermal transport model - integer(kind(DAMAGE_none_ID)), dimension(:), allocatable, public, protected :: & - damage_type !< nonlocal damage model - integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:), allocatable, public, protected :: & - homogenization_type !< type of each homogenization - integer, public, protected :: & homogenization_maxNconstituents !< max number of grains in any USED homogenization - integer, dimension(:), allocatable, public, protected :: & - homogenization_Nconstituents, & !< number of grains in each homogenization - homogenization_typeInstance, & !< instance of particular type of each homogenization - thermal_typeInstance, & !< instance of particular type of each thermal transport - damage_typeInstance !< instance of particular type of each nonlocal damage - - real(pReal), dimension(:), allocatable, public, protected :: & - thermal_initialT !< initial temperature per each homogenization - integer, dimension(:), allocatable, public, protected :: & ! (elem) material_homogenizationAt, & !< homogenization ID of each element material_homogenizationAt2, & !< per cell @@ -63,50 +39,28 @@ module material integer, dimension(:,:,:), allocatable, public, protected :: & ! (constituent,IP,elem) material_phaseMemberAt !< position of the element within its phase instance - type(tState), allocatable, dimension(:), public :: & - homogState, & - damageState_h - type(Rotation), dimension(:,:,:), allocatable, public, protected :: & material_orientation0 !< initial orientation of each grain,IP,element - type(group_float), allocatable, dimension(:), public :: & - damage !< damage field - public :: & - material_init, & - THERMAL_ISOTHERMAL_ID, & - THERMAL_CONDUCTION_ID, & - DAMAGE_NONE_ID, & - DAMAGE_NONLOCAL_ID, & - HOMOGENIZATION_NONE_ID, & - HOMOGENIZATION_ISOSTRAIN_ID, & - HOMOGENIZATION_RGC_ID + material_init contains !-------------------------------------------------------------------------------------------------- -!> @brief parses material configuration file +!> @brief Parse material configuration file (material.yaml). !-------------------------------------------------------------------------------------------------- subroutine material_init(restart) logical, intent(in) :: restart + print'(/,a)', ' <<<+- material init -+>>>'; flush(IO_STDOUT) call material_parseMaterial print*, 'Material parsed' - call material_parseHomogenization - print*, 'Homogenization parsed' - - - allocate(homogState (size(material_name_homogenization))) - allocate(damageState_h (size(material_name_homogenization))) - - allocate(damage (size(material_name_homogenization))) - if (.not. restart) then call results_openJobFile @@ -118,82 +72,6 @@ subroutine material_init(restart) end subroutine material_init -!-------------------------------------------------------------------------------------------------- -!> @brief parses the homogenization part from the material configuration -! ToDo: This should be done in homogenization -!-------------------------------------------------------------------------------------------------- -subroutine material_parseHomogenization - - class(tNode), pointer :: & - material_homogenization, & - homog, & - homogMech, & - homogThermal, & - homogDamage - - integer :: h - - material_homogenization => config_material%get('homogenization') - - allocate(homogenization_type(size(material_name_homogenization)), source=HOMOGENIZATION_undefined_ID) - allocate(thermal_type(size(material_name_homogenization)), source=THERMAL_isothermal_ID) - allocate(damage_type (size(material_name_homogenization)), source=DAMAGE_none_ID) - allocate(homogenization_typeInstance(size(material_name_homogenization)), source=0) - allocate(thermal_typeInstance(size(material_name_homogenization)), source=0) - allocate(damage_typeInstance(size(material_name_homogenization)), source=0) - allocate(thermal_initialT(size(material_name_homogenization)), source=300.0_pReal) - - do h=1, size(material_name_homogenization) - homog => material_homogenization%get(h) - homogMech => homog%get('mechanics') - select case (homogMech%get_asString('type')) - case('pass') - homogenization_type(h) = HOMOGENIZATION_NONE_ID - case('isostrain') - homogenization_type(h) = HOMOGENIZATION_ISOSTRAIN_ID - case('RGC') - homogenization_type(h) = HOMOGENIZATION_RGC_ID - case default - call IO_error(500,ext_msg=homogMech%get_asString('type')) - end select - - homogenization_typeInstance(h) = count(homogenization_type==homogenization_type(h)) - - if(homog%contains('thermal')) then - homogThermal => homog%get('thermal') - thermal_initialT(h) = homogThermal%get_asFloat('T_0',defaultVal=300.0_pReal) - - select case (homogThermal%get_asString('type')) - case('isothermal') - thermal_type(h) = THERMAL_isothermal_ID - case('conduction') - thermal_type(h) = THERMAL_conduction_ID - case default - call IO_error(500,ext_msg=homogThermal%get_asString('type')) - end select - endif - - if(homog%contains('damage')) then - homogDamage => homog%get('damage') - select case (homogDamage%get_asString('type')) - case('none') - damage_type(h) = DAMAGE_none_ID - case('nonlocal') - damage_type(h) = DAMAGE_nonlocal_ID - case default - call IO_error(500,ext_msg=homogDamage%get_asString('type')) - end select - endif - enddo - - do h=1, size(material_name_homogenization) - homogenization_typeInstance(h) = count(homogenization_type(1:h) == homogenization_type(h)) - thermal_typeInstance(h) = count(thermal_type (1:h) == thermal_type (h)) - damage_typeInstance(h) = count(damage_type (1:h) == damage_type (h)) - enddo - -end subroutine material_parseHomogenization - !-------------------------------------------------------------------------------------------------- !> @brief parses the material part in the material configuration file diff --git a/src/phase.f90 b/src/phase.f90 index 628bfc443..f129ba92f 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -58,20 +58,17 @@ module phase type(tDebugOptions) :: debugCrystallite integer, dimension(:), allocatable, public :: & !< ToDo: should be protected (bug in Intel compiler) - thermal_Nsources, & - phase_Nsources, & !< number of source mechanisms active in each phase - phase_Nkinematics, & !< number of kinematic mechanisms active in each phase + phase_elasticityInstance, & phase_NstiffnessDegradations, & !< number of stiffness degradation mechanisms active in each phase - phase_plasticInstance, & !< instance of particular plasticity of each phase - phase_elasticityInstance !< instance of particular elasticity of each phase + phase_plasticInstance logical, dimension(:), allocatable, public :: & ! ToDo: should be protected (bug in Intel Compiler) phase_localPlasticity !< flags phases with local constitutive law type(tPlasticState), allocatable, dimension(:), public :: & plasticState - type(tSourceState), allocatable, dimension(:), public :: & - damageState, thermalState + type(tState), allocatable, dimension(:), public :: & + damageState integer, public, protected :: & @@ -181,6 +178,11 @@ module phase real(pReal) :: dot_T end function thermal_dot_T + module function damage_phi(ph,me) result(phi) + integer, intent(in) :: ph,me + real(pReal) :: phi + end function damage_phi + module subroutine phase_mechanical_setF(F,co,ip,el) real(pReal), dimension(3,3), intent(in) :: F @@ -261,6 +263,27 @@ module phase el !< element end subroutine plastic_dependentState + + module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) + integer, intent(in) :: ph, me + real(pReal), intent(in), dimension(3,3) :: & + S + real(pReal), intent(out), dimension(3,3) :: & + Ld !< damage velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) + end subroutine kinematics_cleavage_opening_LiAndItsTangent + + module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) + integer, intent(in) :: ph, me + real(pReal), intent(in), dimension(3,3) :: & + S + real(pReal), intent(out), dimension(3,3) :: & + Ld !< damage velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) + end subroutine kinematics_slipplane_opening_LiAndItsTangent + end interface @@ -345,14 +368,12 @@ subroutine phase_init PhaseLoop2:do ph = 1,phases%length !-------------------------------------------------------------------------------------------------- ! partition and initialize state - plasticState(ph)%state = plasticState(ph)%state0 - forall(so = 1:phase_Nsources(ph)) - damageState(ph)%p(so)%state = damageState(ph)%p(so)%state0 - end forall - - phase_source_maxSizeDotState = max(phase_source_maxSizeDotState, & - maxval(damageState(ph)%p%sizeDotState)) + plasticState(ph)%state = plasticState(ph)%state0 + if(damageState(ph)%sizeState > 0) & + damageState(ph)%state = damageState(ph)%state0 enddo PhaseLoop2 + + phase_source_maxSizeDotState = maxval(damageState%sizeDotState) phase_plasticity_maxSizeDotState = maxval(plasticState%sizeDotState) end subroutine phase_init @@ -362,7 +383,7 @@ end subroutine phase_init !> @brief Allocate the components of the state structure for a given phase !-------------------------------------------------------------------------------------------------- subroutine phase_allocateState(state, & - Nconstituents,sizeState,sizeDotState,sizeDeltaState) + Nconstituents,sizeState,sizeDotState,sizeDeltaState) class(tState), intent(out) :: & state @@ -399,15 +420,13 @@ subroutine phase_restore(ce,includeL) integer, intent(in) :: ce integer :: & - co, & !< constituent number - so + co do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce)) - do so = 1, phase_Nsources(material_phaseAt2(co,ce)) - damageState(material_phaseAt2(co,ce))%p(so)%state( :,material_phasememberAt2(co,ce)) = & - damageState(material_phaseAt2(co,ce))%p(so)%state0(:,material_phasememberAt2(co,ce)) - enddo + if (damageState(material_phaseAt2(co,ce))%sizeState > 0) & + damageState(material_phaseAt2(co,ce))%state( :,material_phasememberAt2(co,ce)) = & + damageState(material_phaseAt2(co,ce))%state0(:,material_phasememberAt2(co,ce)) enddo call mechanical_restore(ce,includeL) @@ -421,16 +440,16 @@ end subroutine phase_restore !-------------------------------------------------------------------------------------------------- subroutine phase_forward() - integer :: ph, so + integer :: ph call mechanical_forward() call thermal_forward() do ph = 1, size(damageState) - do so = 1,phase_Nsources(ph) - damageState(ph)%p(so)%state0 = damageState(ph)%p(so)%state - enddo; enddo + if (damageState(ph)%sizeState > 0) & + damageState(ph)%state0 = damageState(ph)%state + enddo end subroutine phase_forward @@ -526,9 +545,8 @@ subroutine crystallite_init() phases => config_material%get('phase') do ph = 1, phases%length - do so = 1, phase_Nsources(ph) - allocate(damageState(ph)%p(so)%subState0,source=damageState(ph)%p(so)%state0) ! ToDo: hack - enddo + if (damageState(ph)%sizeState > 0) & + allocate(damageState(ph)%subState0,source=damageState(ph)%state0) ! ToDo: hack enddo print'(a42,1x,i10)', ' # of elements: ', eMax @@ -574,9 +592,8 @@ subroutine phase_windForward(ip,el) call mechanical_windForward(ph,me) - do so = 1, phase_Nsources(material_phaseAt(co,el)) - damageState(ph)%p(so)%state0(:,me) = damageState(ph)%p(so)%state(:,me) - enddo + if(damageState(ph)%sizeState > 0) damageState(ph)%state0(:,me) = damageState(ph)%state(:,me) + enddo diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 57145c550..5e4e43017 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -15,110 +15,92 @@ submodule(phase) damagee real(pReal), dimension(:), allocatable :: phi, d_phi_d_dot_phi end type tDataContainer - integer(kind(DAMAGE_UNDEFINED_ID)), dimension(:,:), allocatable :: & + integer(kind(DAMAGE_UNDEFINED_ID)), dimension(:), allocatable :: & phase_source !< active sources mechanisms of each phase + integer, dimension(:), allocatable :: & + phase_Nsources + type(tDataContainer), dimension(:), allocatable :: current interface - module function anisobrittle_init(source_length) result(mySources) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources + module function anisobrittle_init() result(mySources) + logical, dimension(:), allocatable :: mySources end function anisobrittle_init - module function anisoductile_init(source_length) result(mySources) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources + module function anisoductile_init() result(mySources) + logical, dimension(:), allocatable :: mySources end function anisoductile_init - module function isobrittle_init(source_length) result(mySources) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources + module function isobrittle_init() result(mySources) + logical, dimension(:), allocatable :: mySources end function isobrittle_init - module function isoductile_init(source_length) result(mySources) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources + module function isoductile_init() result(mySources) + logical, dimension(:), allocatable :: mySources end function isoductile_init - module subroutine source_damage_isoBrittle_deltaState(C, Fe, ph, me) + module subroutine isobrittle_deltaState(C, Fe, ph, me) integer, intent(in) :: ph,me real(pReal), intent(in), dimension(3,3) :: & Fe real(pReal), intent(in), dimension(6,6) :: & C - end subroutine source_damage_isoBrittle_deltaState + end subroutine isobrittle_deltaState - module subroutine anisobrittle_dotState(S, co, ip, el) - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element + module subroutine anisobrittle_dotState(S, ph, me) + integer, intent(in) :: ph,me real(pReal), intent(in), dimension(3,3) :: & S end subroutine anisobrittle_dotState - module subroutine anisoductile_dotState(co, ip, el) - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element + module subroutine anisoductile_dotState(ph,me) + integer, intent(in) :: ph,me end subroutine anisoductile_dotState - module subroutine isoductile_dotState(co, ip, el) - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element + module subroutine isoductile_dotState(ph,me) + integer, intent(in) :: ph,me end subroutine isoductile_dotState - module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - integer, intent(in) :: & - phase, & !< phase ID of element - constituent !< position of element within its phase instance + module subroutine anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) + integer, intent(in) :: ph,me real(pReal), intent(in) :: & phi !< damage parameter real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi - end subroutine source_damage_anisoBrittle_getRateAndItsTangent + end subroutine anisobrittle_getRateAndItsTangent - module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - integer, intent(in) :: & - phase, & !< phase ID of element - constituent !< position of element within its phase instance + module subroutine anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me) + integer, intent(in) :: ph,me real(pReal), intent(in) :: & phi !< damage parameter real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi - end subroutine source_damage_anisoDuctile_getRateAndItsTangent + end subroutine anisoductile_getRateAndItsTangent - module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - integer, intent(in) :: & - phase, & !< phase ID of element - constituent !< position of element within its phase instance + module subroutine isobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me) + integer, intent(in) :: ph,me real(pReal), intent(in) :: & phi !< damage parameter real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi - end subroutine source_damage_isoBrittle_getRateAndItsTangent + end subroutine isobrittle_getRateAndItsTangent - module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - integer, intent(in) :: & - phase, & !< phase ID of element - constituent !< position of element within its phase instance + module subroutine isoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me) + integer, intent(in) :: ph,me real(pReal), intent(in) :: & phi !< damage parameter real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi - end subroutine source_damage_isoDuctile_getRateAndItsTangent + end subroutine isoductile_getRateAndItsTangent module subroutine anisobrittle_results(phase,group) integer, intent(in) :: phase @@ -175,19 +157,20 @@ module subroutine damage_init allocate(current(ph)%d_phi_d_dot_phi(Nconstituents),source=0.0_pReal) phase => phases%get(ph) - sources => phase%get('source',defaultVal=emptyList) + sources => phase%get('damage',defaultVal=emptyList) + if (sources%length > 1) error stop phase_Nsources(ph) = sources%length - allocate(damageState(ph)%p(phase_Nsources(ph))) + enddo - allocate(phase_source(maxval(phase_Nsources),phases%length), source = DAMAGE_UNDEFINED_ID) + allocate(phase_source(phases%length), source = DAMAGE_UNDEFINED_ID) ! initialize source mechanisms if(maxval(phase_Nsources) /= 0) then - where(isobrittle_init (maxval(phase_Nsources))) phase_source = DAMAGE_ISOBRITTLE_ID - where(isoductile_init (maxval(phase_Nsources))) phase_source = DAMAGE_ISODUCTILE_ID - where(anisobrittle_init (maxval(phase_Nsources))) phase_source = DAMAGE_ANISOBRITTLE_ID - where(anisoductile_init (maxval(phase_Nsources))) phase_source = DAMAGE_ANISODUCTILE_ID + where(isobrittle_init() ) phase_source = DAMAGE_ISOBRITTLE_ID + where(isoductile_init() ) phase_source = DAMAGE_ISODUCTILE_ID + where(anisobrittle_init()) phase_source = DAMAGE_ANISOBRITTLE_ID + where(anisoductile_init()) phase_source = DAMAGE_ANISODUCTILE_ID endif end subroutine damage_init @@ -213,7 +196,6 @@ module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, integer :: & ph, & co, & - so, & me phiDot = 0.0_pReal @@ -222,19 +204,19 @@ module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) ph = material_phaseAt(co,el) me = material_phasememberAt(co,ip,el) - do so = 1, phase_Nsources(ph) - select case(phase_source(so,ph)) + + select case(phase_source(ph)) case (DAMAGE_ISOBRITTLE_ID) - call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, ph, me) + call isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, ph, me) case (DAMAGE_ISODUCTILE_ID) - call source_damage_isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, ph, me) + call isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, ph, me) case (DAMAGE_ANISOBRITTLE_ID) - call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) + call anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) case (DAMAGE_ANISODUCTILE_ID) - call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) + call anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) case default localphiDot = 0.0_pReal @@ -243,7 +225,6 @@ module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, end select phiDot = phiDot + localphiDot dPhiDot_dPhi = dPhiDot_dPhi + dLocalphiDot_dPhi - enddo enddo end subroutine phase_damage_getRateAndItsTangents @@ -267,57 +248,52 @@ module function integrateDamageState(dt,co,ip,el) result(broken) NiterationState, & !< number of iterations in state loop ph, & me, & - so - integer, dimension(maxval(phase_Nsources)) :: & size_so real(pReal) :: & zeta real(pReal), dimension(phase_source_maxSizeDotState) :: & r ! state residuum - real(pReal), dimension(phase_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState + real(pReal), dimension(phase_source_maxSizeDotState,2) :: source_dotState logical :: & converged_ - ph = material_phaseAt(co,el) me = material_phaseMemberAt(co,ip,el) + if (damageState(ph)%sizeState == 0) then + broken = .false. + return + endif + converged_ = .true. - broken = phase_damage_collectDotState(co,ip,el,ph,me) + broken = phase_damage_collectDotState(ph,me) if(broken) return - do so = 1, phase_Nsources(ph) - size_so(so) = damageState(ph)%p(so)%sizeDotState - damageState(ph)%p(so)%state(1:size_so(so),me) = damageState(ph)%p(so)%subState0(1:size_so(so),me) & - + damageState(ph)%p(so)%dotState (1:size_so(so),me) * dt - source_dotState(1:size_so(so),2,so) = 0.0_pReal - enddo + size_so = damageState(ph)%sizeDotState + damageState(ph)%state(1:size_so,me) = damageState(ph)%subState0(1:size_so,me) & + + damageState(ph)%dotState (1:size_so,me) * dt + source_dotState(1:size_so,2) = 0.0_pReal iteration: do NiterationState = 1, num%nState - do so = 1, phase_Nsources(ph) - if(nIterationState > 1) source_dotState(1:size_so(so),2,so) = source_dotState(1:size_so(so),1,so) - source_dotState(1:size_so(so),1,so) = damageState(ph)%p(so)%dotState(:,me) - enddo + if(nIterationState > 1) source_dotState(1:size_so,2) = source_dotState(1:size_so,1) + source_dotState(1:size_so,1) = damageState(ph)%dotState(:,me) - broken = phase_damage_collectDotState(co,ip,el,ph,me) + broken = phase_damage_collectDotState(ph,me) if(broken) exit iteration - do so = 1, phase_Nsources(ph) - zeta = damper(damageState(ph)%p(so)%dotState(:,me), & - source_dotState(1:size_so(so),1,so),& - source_dotState(1:size_so(so),2,so)) - damageState(ph)%p(so)%dotState(:,me) = damageState(ph)%p(so)%dotState(:,me) * zeta & - + source_dotState(1:size_so(so),1,so)* (1.0_pReal - zeta) - r(1:size_so(so)) = damageState(ph)%p(so)%state (1:size_so(so),me) & - - damageState(ph)%p(so)%subState0(1:size_so(so),me) & - - damageState(ph)%p(so)%dotState (1:size_so(so),me) * dt - damageState(ph)%p(so)%state(1:size_so(so),me) = damageState(ph)%p(so)%state(1:size_so(so),me) & - - r(1:size_so(so)) - converged_ = converged_ .and. converged(r(1:size_so(so)), & - damageState(ph)%p(so)%state(1:size_so(so),me), & - damageState(ph)%p(so)%atol(1:size_so(so))) - enddo + + zeta = damper(damageState(ph)%dotState(:,me),source_dotState(1:size_so,1),source_dotState(1:size_so,2)) + damageState(ph)%dotState(:,me) = damageState(ph)%dotState(:,me) * zeta & + + source_dotState(1:size_so,1)* (1.0_pReal - zeta) + r(1:size_so) = damageState(ph)%state (1:size_so,me) & + - damageState(ph)%subState0(1:size_so,me) & + - damageState(ph)%dotState (1:size_so,me) * dt + damageState(ph)%state(1:size_so,me) = damageState(ph)%state(1:size_so,me) - r(1:size_so) + converged_ = converged_ .and. converged(r(1:size_so), & + damageState(ph)%state(1:size_so,me), & + damageState(ph)%atol(1:size_so)) + if(converged_) then broken = phase_damage_deltaState(mechanical_F_e(ph,me),ph,me) @@ -366,10 +342,10 @@ module subroutine damage_results(group,ph) sourceLoop: do so = 1, phase_Nsources(ph) - if (phase_source(so,ph) /= DAMAGE_UNDEFINED_ID) & + if (phase_source(ph) /= DAMAGE_UNDEFINED_ID) & call results_closeGroup(results_addGroup(group//'sources/')) ! should be 'damage' - sourceType: select case (phase_source(so,ph)) + sourceType: select case (phase_source(ph)) case (DAMAGE_ISOBRITTLE_ID) sourceType call isobrittle_results(ph,group//'sources/') @@ -393,39 +369,34 @@ end subroutine damage_results !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -function phase_damage_collectDotState(co,ip,el,ph,me) result(broken) +function phase_damage_collectDotState(ph,me) result(broken) integer, intent(in) :: & - co, & !< component-ID me integration point - ip, & !< integration point - el, & !< element ph, & - me - integer :: & - so !< counter in source loop + me !< counter in source loop logical :: broken broken = .false. - SourceLoop: do so = 1, phase_Nsources(ph) + if (damageState(ph)%sizeState > 0) then - sourceType: select case (phase_source(so,ph)) + sourceType: select case (phase_source(ph)) case (DAMAGE_ISODUCTILE_ID) sourceType - call isoductile_dotState(co, ip, el) + call isoductile_dotState(ph,me) case (DAMAGE_ANISODUCTILE_ID) sourceType - call anisoductile_dotState(co, ip, el) + call anisoductile_dotState(ph,me) case (DAMAGE_ANISOBRITTLE_ID) sourceType - call anisobrittle_dotState(mechanical_S(ph,me),co, ip, el) ! correct stress? + call anisobrittle_dotState(mechanical_S(ph,me), ph,me) ! correct stress? end select sourceType - broken = broken .or. any(IEEE_is_NaN(damageState(ph)%p(so)%dotState(:,me))) + broken = broken .or. any(IEEE_is_NaN(damageState(ph)%dotState(:,me))) - enddo SourceLoop + endif end function phase_damage_collectDotState @@ -443,7 +414,6 @@ function phase_damage_deltaState(Fe, ph, me) result(broken) real(pReal), intent(in), dimension(3,3) :: & Fe !< elastic deformation gradient integer :: & - so, & myOffset, & mySize logical :: & @@ -452,23 +422,22 @@ function phase_damage_deltaState(Fe, ph, me) result(broken) broken = .false. - sourceLoop: do so = 1, phase_Nsources(ph) + if (damageState(ph)%sizeState == 0) return - sourceType: select case (phase_source(so,ph)) + sourceType: select case (phase_source(ph)) case (DAMAGE_ISOBRITTLE_ID) sourceType - call source_damage_isoBrittle_deltaState(phase_homogenizedC(ph,me), Fe, ph,me) - broken = any(IEEE_is_NaN(damageState(ph)%p(so)%deltaState(:,me))) + call isobrittle_deltaState(phase_homogenizedC(ph,me), Fe, ph,me) + broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,me))) if(.not. broken) then - myOffset = damageState(ph)%p(so)%offsetDeltaState - mySize = damageState(ph)%p(so)%sizeDeltaState - damageState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,me) = & - damageState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,me) + damageState(ph)%p(so)%deltaState(1:mySize,me) + myOffset = damageState(ph)%offsetDeltaState + mySize = damageState(ph)%sizeDeltaState + damageState(ph)%state(myOffset + 1: myOffset + mySize,me) = & + damageState(ph)%state(myOffset + 1: myOffset + mySize,me) + damageState(ph)%deltaState(1:mySize,me) endif end select sourceType - enddo SourceLoop end function phase_damage_deltaState @@ -476,28 +445,25 @@ end function phase_damage_deltaState !-------------------------------------------------------------------------------------------------- !> @brief checks if a source mechanism is active or not !-------------------------------------------------------------------------------------------------- -function source_active(source_label,src_length) result(active_source) +function source_active(source_label) result(active_source) character(len=*), intent(in) :: source_label !< name of source mechanism - integer, intent(in) :: src_length !< max. number of sources in system - logical, dimension(:,:), allocatable :: active_source + logical, dimension(:), allocatable :: active_source class(tNode), pointer :: & phases, & phase, & sources, & src - integer :: p,s + integer :: ph phases => config_material%get('phase') - allocate(active_source(src_length,phases%length), source = .false. ) - do p = 1, phases%length - phase => phases%get(p) - sources => phase%get('source',defaultVal=emptyList) - do s = 1, sources%length - src => sources%get(s) - if(src%get_asString('type') == source_label) active_source(s,p) = .true. - enddo + allocate(active_source(phases%length)) + do ph = 1, phases%length + phase => phases%get(ph) + sources => phase%get('damage',defaultVal=emptyList) + src => sources%get(1) + active_source(ph) = src%get_asString('type',defaultVal = 'x') == source_label enddo @@ -528,4 +494,15 @@ module function phase_damage_get_phi(co,ip,el) result(phi) end function phase_damage_get_phi +module function damage_phi(ph,me) result(phi) + + integer, intent(in) :: ph, me + real(pReal) :: phi + + + phi = current(ph)%phi(me) + +end function damage_phi + + end submodule damagee diff --git a/src/phase_damage_anisobrittle.f90 b/src/phase_damage_anisobrittle.f90 index 51ca7786a..fd82b74c4 100644 --- a/src/phase_damage_anisobrittle.f90 +++ b/src/phase_damage_anisobrittle.f90 @@ -6,10 +6,6 @@ !-------------------------------------------------------------------------------------------------- submodule (phase:damagee) anisobrittle - integer, dimension(:), allocatable :: & - source_damage_anisoBrittle_offset, & !< which source is my current source mechanism? - source_damage_anisoBrittle_instance !< instance of source mechanism - type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & dot_o, & !< opening rate of cleavage planes @@ -35,42 +31,38 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function anisobrittle_init(source_length) result(mySources) +module function anisobrittle_init() result(mySources) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources + logical, dimension(:), allocatable :: mySources class(tNode), pointer :: & phases, & phase, & sources, & src - integer :: Ninstances,sourceOffset,Nconstituents,p + integer :: Nconstituents,p integer, dimension(:), allocatable :: N_cl character(len=pStringLen) :: extmsg = '' - print'(/,a)', ' <<<+- phase:damage:anisobrittle init -+>>>' - mySources = source_active('damage_anisoBrittle',source_length) - Ninstances = count(mySources) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) - if(Ninstances == 0) return + mySources = source_active('anisobrittle') + if(count(mySources) == 0) return + + print'(/,a)', ' <<<+- phase:damage:anisobrittle init -+>>>' + print'(a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT) + phases => config_material%get('phase') - allocate(param(Ninstances)) - allocate(source_damage_anisoBrittle_offset (phases%length), source=0) - allocate(source_damage_anisoBrittle_instance(phases%length), source=0) + allocate(param(phases%length)) + do p = 1, phases%length + if(mySources(p)) then phase => phases%get(p) - if(any(mySources(:,p))) source_damage_anisoBrittle_instance(p) = count(mySources(:,1:p)) - if(count(mySources(:,p)) == 0) cycle - sources => phase%get('source') - do sourceOffset = 1, sources%length - if(mySources(sourceOffset,p)) then - source_damage_anisoBrittle_offset(p) = sourceOffset - associate(prm => param(source_damage_anisoBrittle_instance(p))) - src => sources%get(sourceOffset) + sources => phase%get('damage') + + associate(prm => param(p)) + src => sources%get(1) N_cl = src%get_asInts('N_cl',defaultVal=emptyIntArray) prm%sum_N_cl = sum(abs(N_cl)) @@ -101,9 +93,9 @@ module function anisobrittle_init(source_length) result(mySources) if (any(prm%s_crit < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit' Nconstituents = count(material_phaseAt==p) * discretization_nIPs - call phase_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0) - damageState(p)%p(sourceOffset)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal) - if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol' + call phase_allocateState(damageState(p),Nconstituents,1,1,0) + damageState(p)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol' end associate @@ -111,7 +103,7 @@ module function anisobrittle_init(source_length) result(mySources) ! exit if any parameter is out of range if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_anisoBrittle)') endif - enddo + enddo end function anisobrittle_init @@ -120,18 +112,14 @@ end function anisobrittle_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -module subroutine anisobrittle_dotState(S, co, ip, el) +module subroutine anisobrittle_dotState(S, ph,me) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element + ph,me real(pReal), intent(in), dimension(3,3) :: & S integer :: & - ph, & - me, & sourceOffset, & damageOffset, & homog, & @@ -139,28 +127,22 @@ module subroutine anisobrittle_dotState(S, co, ip, el) real(pReal) :: & traction_d, traction_t, traction_n, traction_crit - ph = material_phaseAt(co,el) - me = material_phasememberAt(co,ip,el) - sourceOffset = source_damage_anisoBrittle_offset(ph) - homog = material_homogenizationAt(el) - damageOffset = material_homogenizationMemberAt(ip,el) - associate(prm => param(source_damage_anisoBrittle_instance(ph))) - damageState(ph)%p(sourceOffset)%dotState(1,me) = 0.0_pReal - do i = 1, prm%sum_N_cl - traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i)) - traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i)) - traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i)) + associate(prm => param(ph)) + damageState(ph)%dotState(1,me) = 0.0_pReal + do i = 1, prm%sum_N_cl + traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i)) + traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i)) + traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i)) - traction_crit = prm%g_crit(i)*damage(homog)%p(damageOffset)**2.0_pReal + traction_crit = prm%g_crit(i)*damage_phi(ph,me)**2.0_pReal - damageState(ph)%p(sourceOffset)%dotState(1,me) & - = damageState(ph)%p(sourceOffset)%dotState(1,me) & - + prm%dot_o / prm%s_crit(i) & - * ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%q + & - (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%q + & - (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**prm%q) - enddo + damageState(ph)%dotState(1,me) = damageState(ph)%dotState(1,me) & + + prm%dot_o / prm%s_crit(i) & + * ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%q + & + (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%q + & + (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**prm%q) + enddo end associate end subroutine anisobrittle_dotState @@ -169,28 +151,24 @@ end subroutine anisobrittle_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- -module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) +module subroutine anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) integer, intent(in) :: & - phase, & - constituent + ph, & + me real(pReal), intent(in) :: & phi real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi - integer :: & - sourceOffset - sourceOffset = source_damage_anisoBrittle_offset(phase) - - dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent) + dLocalphiDot_dPhi = -damageState(ph)%state(1,me) localphiDot = 1.0_pReal & + dLocalphiDot_dPhi*phi -end subroutine source_damage_anisoBrittle_getRateAndItsTangent +end subroutine anisobrittle_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- @@ -203,16 +181,78 @@ module subroutine anisobrittle_results(phase,group) integer :: o - associate(prm => param(source_damage_anisoBrittle_instance(phase)), & - stt => damageState(phase)%p(source_damage_anisoBrittle_offset(phase))%state) - outputsLoop: do o = 1,size(prm%output) - select case(trim(prm%output(o))) - case ('f_phi') - call results_writeDataset(group,stt,trim(prm%output(o)),'driving force','J/m³') - end select - enddo outputsLoop + + associate(prm => param(phase), stt => damageState(phase)%state) + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + case ('f_phi') + call results_writeDataset(group,stt,trim(prm%output(o)),'driving force','J/m³') + end select + enddo outputsLoop end associate end subroutine anisobrittle_results + +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the velocity gradient +!-------------------------------------------------------------------------------------------------- +module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) + + integer, intent(in) :: & + ph,me + real(pReal), intent(in), dimension(3,3) :: & + S + real(pReal), intent(out), dimension(3,3) :: & + Ld !< damage velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) + + integer :: & + i, k, l, m, n + real(pReal) :: & + traction_d, traction_t, traction_n, traction_crit, & + udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt + + + Ld = 0.0_pReal + dLd_dTstar = 0.0_pReal + associate(prm => param(ph)) + do i = 1,prm%sum_N_cl + traction_crit = prm%g_crit(i)*damage_phi(ph,me)**2.0_pReal + + traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i)) + if (abs(traction_d) > traction_crit + tol_math_check) then + udotd = sign(1.0_pReal,traction_d)* prm%dot_o * ((abs(traction_d) - traction_crit)/traction_crit)**prm%q + Ld = Ld + udotd*prm%cleavage_systems(1:3,1:3,1,i) + dudotd_dt = sign(1.0_pReal,traction_d)*udotd*prm%q / (abs(traction_d) - traction_crit) + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) & + + dudotd_dt*prm%cleavage_systems(k,l,1,i) * prm%cleavage_systems(m,n,1,i) + endif + + traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i)) + if (abs(traction_t) > traction_crit + tol_math_check) then + udott = sign(1.0_pReal,traction_t)* prm%dot_o * ((abs(traction_t) - traction_crit)/traction_crit)**prm%q + Ld = Ld + udott*prm%cleavage_systems(1:3,1:3,2,i) + dudott_dt = sign(1.0_pReal,traction_t)*udott*prm%q / (abs(traction_t) - traction_crit) + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) & + + dudott_dt*prm%cleavage_systems(k,l,2,i) * prm%cleavage_systems(m,n,2,i) + endif + + traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i)) + if (abs(traction_n) > traction_crit + tol_math_check) then + udotn = sign(1.0_pReal,traction_n)* prm%dot_o * ((abs(traction_n) - traction_crit)/traction_crit)**prm%q + Ld = Ld + udotn*prm%cleavage_systems(1:3,1:3,3,i) + dudotn_dt = sign(1.0_pReal,traction_n)*udotn*prm%q / (abs(traction_n) - traction_crit) + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) & + + dudotn_dt*prm%cleavage_systems(k,l,3,i) * prm%cleavage_systems(m,n,3,i) + endif + enddo + end associate + +end subroutine kinematics_cleavage_opening_LiAndItsTangent + end submodule anisobrittle diff --git a/src/phase_damage_anisoductile.f90 b/src/phase_damage_anisoductile.f90 index e54eff201..f50c05f07 100644 --- a/src/phase_damage_anisoductile.f90 +++ b/src/phase_damage_anisoductile.f90 @@ -6,10 +6,6 @@ !-------------------------------------------------------------------------------------------------- submodule(phase:damagee) anisoductile - integer, dimension(:), allocatable :: & - source_damage_anisoDuctile_offset, & !< which source is my current damage mechanism? - source_damage_anisoDuctile_instance !< instance of damage source mechanism - type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & q !< damage rate sensitivity @@ -19,7 +15,7 @@ submodule(phase:damagee) anisoductile output end type tParameters - type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances) + type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters contains @@ -28,10 +24,9 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function anisoductile_init(source_length) result(mySources) +module function anisoductile_init() result(mySources) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources + logical, dimension(:), allocatable :: mySources class(tNode), pointer :: & phases, & @@ -40,34 +35,31 @@ module function anisoductile_init(source_length) result(mySources) pl, & sources, & src - integer :: Ninstances,sourceOffset,Nconstituents,p + integer :: Ninstances,Nconstituents,p integer, dimension(:), allocatable :: N_sl character(len=pStringLen) :: extmsg = '' - print'(/,a)', ' <<<+- phase:damage:anisoductile init -+>>>' - mySources = source_active('damage_anisoDuctile',source_length) - Ninstances = count(mySources) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) - if(Ninstances == 0) return + mySources = source_active('anisoductile') + if(count(mySources) == 0) return + + print'(/,a)', ' <<<+- phase:damage:anisoductile init -+>>>' + print'(a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT) + phases => config_material%get('phase') - allocate(param(Ninstances)) - allocate(source_damage_anisoDuctile_offset (phases%length), source=0) - allocate(source_damage_anisoDuctile_instance(phases%length), source=0) + allocate(param(phases%length)) do p = 1, phases%length - phase => phases%get(p) - if(any(mySources(:,p))) source_damage_anisoDuctile_instance(p) = count(mySources(:,1:p)) - if(count(mySources(:,p)) == 0) cycle - mech => phase%get('mechanics') - pl => mech%get('plasticity') - sources => phase%get('source') - do sourceOffset = 1, sources%length - if(mySources(sourceOffset,p)) then - source_damage_anisoDuctile_offset(p) = sourceOffset - associate(prm => param(source_damage_anisoDuctile_instance(p))) - src => sources%get(sourceOffset) + if(mySources(p)) then + phase => phases%get(p) + mech => phase%get('mechanics') + pl => mech%get('plasticity') + sources => phase%get('damage') + + + associate(prm => param(p)) + src => sources%get(1) N_sl = pl%get_asInts('N_sl',defaultVal=emptyIntArray) prm%q = src%get_asFloat('q') @@ -86,10 +78,10 @@ module function anisoductile_init(source_length) result(mySources) if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit' - Nconstituents=count(material_phaseAt==p) * discretization_nIPs - call phase_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0) - damageState(p)%p(sourceOffset)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal) - if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol' + Nconstituents=count(material_phaseAt2==p) + call phase_allocateState(damageState(p),Nconstituents,1,1,0) + damageState(p)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol' end associate @@ -97,7 +89,7 @@ module function anisoductile_init(source_length) result(mySources) ! exit if any parameter is out of range if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_anisoDuctile)') endif - enddo + enddo @@ -107,29 +99,15 @@ end function anisoductile_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -module subroutine anisoductile_dotState(co, ip, el) +module subroutine anisoductile_dotState(ph,me) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element - - integer :: & ph, & - me, & - sourceOffset, & - damageOffset, & - homog + me - ph = material_phaseAt(co,el) - me = material_phasememberAt(co,ip,el) - sourceOffset = source_damage_anisoDuctile_offset(ph) - homog = material_homogenizationAt(el) - damageOffset = material_homogenizationMemberAt(ip,el) - associate(prm => param(source_damage_anisoDuctile_instance(ph))) - damageState(ph)%p(sourceOffset)%dotState(1,me) & - = sum(plasticState(ph)%slipRate(:,me)/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit) + associate(prm => param(ph)) + damageState(ph)%dotState(1,me) = sum(plasticState(ph)%slipRate(:,me)/(damage_phi(ph,me)**prm%q)/prm%gamma_crit) end associate end subroutine anisoductile_dotState @@ -138,28 +116,24 @@ end subroutine anisoductile_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- -module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) +module subroutine anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me) integer, intent(in) :: & - phase, & - constituent + ph, & + me real(pReal), intent(in) :: & phi real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi - integer :: & - sourceOffset - sourceOffset = source_damage_anisoDuctile_offset(phase) - - dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent) + dLocalphiDot_dPhi = -damageState(ph)%state(1,me) localphiDot = 1.0_pReal & + dLocalphiDot_dPhi*phi -end subroutine source_damage_anisoDuctile_getRateAndItsTangent +end subroutine anisoductile_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- @@ -172,14 +146,14 @@ module subroutine anisoductile_results(phase,group) integer :: o - associate(prm => param(source_damage_anisoDuctile_instance(phase)), & - stt => damageState(phase)%p(source_damage_anisoDuctile_offset(phase))%state) - outputsLoop: do o = 1,size(prm%output) - select case(trim(prm%output(o))) - case ('f_phi') - call results_writeDataset(group,stt,trim(prm%output(o)),'driving force','J/m³') - end select - enddo outputsLoop + + associate(prm => param(phase), stt => damageState(phase)%state) + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + case ('f_phi') + call results_writeDataset(group,stt,trim(prm%output(o)),'driving force','J/m³') + end select + enddo outputsLoop end associate end subroutine anisoductile_results diff --git a/src/phase_damage_isobrittle.f90 b/src/phase_damage_isobrittle.f90 index 529ecd404..1cc0d7900 100644 --- a/src/phase_damage_isobrittle.f90 +++ b/src/phase_damage_isobrittle.f90 @@ -6,10 +6,6 @@ !-------------------------------------------------------------------------------------------------- submodule(phase:damagee) isobrittle - integer, dimension(:), allocatable :: & - source_damage_isoBrittle_offset, & - source_damage_isoBrittle_instance - type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & W_crit !< critical elastic strain energy @@ -26,41 +22,36 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function isobrittle_init(source_length) result(mySources) +module function isobrittle_init() result(mySources) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources + logical, dimension(:), allocatable :: mySources class(tNode), pointer :: & phases, & phase, & sources, & src - integer :: Ninstances,sourceOffset,Nconstituents,p + integer :: Nconstituents,p character(len=pStringLen) :: extmsg = '' - print'(/,a)', ' <<<+- phase:damage:isobrittle init -+>>>' - mySources = source_active('damage_isoBrittle',source_length) - Ninstances = count(mySources) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) - if(Ninstances == 0) return + mySources = source_active('isobrittle') + if(count(mySources) == 0) return + + print'(/,a)', ' <<<+- phase:damage:isobrittle init -+>>>' + print'(a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT) + phases => config_material%get('phase') - allocate(param(Ninstances)) - allocate(source_damage_isoBrittle_offset (phases%length), source=0) - allocate(source_damage_isoBrittle_instance(phases%length), source=0) + allocate(param(phases%length)) do p = 1, phases%length + if(mySources(p)) then phase => phases%get(p) - if(any(mySources(:,p))) source_damage_isoBrittle_instance(p) = count(mySources(:,1:p)) - if(count(mySources(:,p)) == 0) cycle - sources => phase%get('source') - do sourceOffset = 1, sources%length - if(mySources(sourceOffset,p)) then - source_damage_isoBrittle_offset(p) = sourceOffset - associate(prm => param(source_damage_isoBrittle_instance(p))) - src => sources%get(sourceOffset) + sources => phase%get('damage') + + associate(prm => param(p)) + src => sources%get(1) prm%W_crit = src%get_asFloat('W_crit') @@ -73,10 +64,10 @@ module function isobrittle_init(source_length) result(mySources) ! sanity checks if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit' - Nconstituents = count(material_phaseAt==p) * discretization_nIPs - call phase_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,1) - damageState(p)%p(sourceOffset)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal) - if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' + Nconstituents = count(material_phaseAt2==p) + call phase_allocateState(damageState(p),Nconstituents,1,1,1) + damageState(p)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' end associate @@ -84,7 +75,7 @@ module function isobrittle_init(source_length) result(mySources) ! exit if any parameter is out of range if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isoBrittle)') endif - enddo + enddo @@ -94,7 +85,7 @@ end function isobrittle_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -module subroutine source_damage_isoBrittle_deltaState(C, Fe, ph,me) +module subroutine isobrittle_deltaState(C, Fe, ph,me) integer, intent(in) :: ph,me real(pReal), intent(in), dimension(3,3) :: & @@ -102,60 +93,47 @@ module subroutine source_damage_isoBrittle_deltaState(C, Fe, ph,me) real(pReal), intent(in), dimension(6,6) :: & C - integer :: & - sourceOffset real(pReal), dimension(6) :: & strain real(pReal) :: & strainenergy - sourceOffset = source_damage_isoBrittle_offset(ph) strain = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3) - associate(prm => param(source_damage_isoBrittle_instance(ph))) - strainenergy = 2.0_pReal*sum(strain*matmul(C,strain))/prm%W_crit - ! ToDo: check strainenergy = 2.0_pReal*dot_product(strain,matmul(C,strain))/prm%W_crit + associate(prm => param(ph)) + strainenergy = 2.0_pReal*sum(strain*matmul(C,strain))/prm%W_crit + ! ToDo: check strainenergy = 2.0_pReal*dot_product(strain,matmul(C,strain))/prm%W_crit - if (strainenergy > damageState(ph)%p(sourceOffset)%subState0(1,me)) then - damageState(ph)%p(sourceOffset)%deltaState(1,me) = & - strainenergy - damageState(ph)%p(sourceOffset)%state(1,me) - else - damageState(ph)%p(sourceOffset)%deltaState(1,me) = & - damageState(ph)%p(sourceOffset)%subState0(1,me) - & - damageState(ph)%p(sourceOffset)%state(1,me) - endif + damageState(ph)%deltaState(1,me) = merge(strainenergy - damageState(ph)%state(1,me), & + damageState(ph)%subState0(1,me) - damageState(ph)%state(1,me), & + strainenergy > damageState(ph)%subState0(1,me)) end associate -end subroutine source_damage_isoBrittle_deltaState +end subroutine isobrittle_deltaState !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- -module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) +module subroutine isobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) integer, intent(in) :: & - phase, & - constituent + ph, me real(pReal), intent(in) :: & phi real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi - integer :: & - sourceOffset - sourceOffset = source_damage_isoBrittle_offset(phase) - - associate(prm => param(source_damage_isoBrittle_instance(phase))) - localphiDot = 1.0_pReal & - - phi*damageState(phase)%p(sourceOffset)%state(1,constituent) - dLocalphiDot_dPhi = - damageState(phase)%p(sourceOffset)%state(1,constituent) + associate(prm => param(ph)) + localphiDot = 1.0_pReal & + - phi*damageState(ph)%state(1,me) + dLocalphiDot_dPhi = - damageState(ph)%state(1,me) end associate -end subroutine source_damage_isoBrittle_getRateAndItsTangent +end subroutine isobrittle_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- @@ -168,14 +146,15 @@ module subroutine isobrittle_results(phase,group) integer :: o - associate(prm => param(source_damage_isoBrittle_instance(phase)), & - stt => damageState(phase)%p(source_damage_isoBrittle_offset(phase))%state) - outputsLoop: do o = 1,size(prm%output) - select case(trim(prm%output(o))) - case ('f_phi') - call results_writeDataset(group,stt,trim(prm%output(o)),'driving force','J/m³') - end select - enddo outputsLoop + + associate(prm => param(phase), & + stt => damageState(phase)%state) + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + case ('f_phi') + call results_writeDataset(group,stt,trim(prm%output(o)),'driving force','J/m³') + end select + enddo outputsLoop end associate end subroutine isobrittle_results diff --git a/src/phase_damage_isoductile.f90 b/src/phase_damage_isoductile.f90 index 71f99078b..247cd715a 100644 --- a/src/phase_damage_isoductile.f90 +++ b/src/phase_damage_isoductile.f90 @@ -6,10 +6,6 @@ !-------------------------------------------------------------------------------------------------- submodule(phase:damagee) isoductile - integer, dimension(:), allocatable :: & - source_damage_isoDuctile_offset, & !< which source is my current damage mechanism? - source_damage_isoDuctile_instance !< instance of damage source mechanism - type:: tParameters !< container type for internal constitutive parameters real(pReal) :: & gamma_crit, & !< critical plastic strain @@ -28,41 +24,36 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function isoductile_init(source_length) result(mySources) +module function isoductile_init() result(mySources) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources + logical, dimension(:), allocatable :: mySources class(tNode), pointer :: & phases, & phase, & sources, & src - integer :: Ninstances,sourceOffset,Nconstituents,p + integer :: Ninstances,Nconstituents,p character(len=pStringLen) :: extmsg = '' - print'(/,a)', ' <<<+- phase:damage:isoductile init -+>>>' - mySources = source_active('damage_isoDuctile',source_length) - Ninstances = count(mySources) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) - if(Ninstances == 0) return + mySources = source_active('isoductile') + if(count(mySources) == 0) return + + print'(/,a)', ' <<<+- phase:damage:isoductile init -+>>>' + print'(a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT) + phases => config_material%get('phase') - allocate(param(Ninstances)) - allocate(source_damage_isoDuctile_offset (phases%length), source=0) - allocate(source_damage_isoDuctile_instance(phases%length), source=0) + allocate(param(phases%length)) do p = 1, phases%length - phase => phases%get(p) - if(count(mySources(:,p)) == 0) cycle - if(any(mySources(:,p))) source_damage_isoDuctile_instance(p) = count(mySources(:,1:p)) - sources => phase%get('source') - do sourceOffset = 1, sources%length - if(mySources(sourceOffset,p)) then - source_damage_isoDuctile_offset(p) = sourceOffset - associate(prm => param(source_damage_isoDuctile_instance(p))) - src => sources%get(sourceOffset) + if(mySources(p)) then + phase => phases%get(p) + sources => phase%get('damage') + + associate(prm => param(p)) + src => sources%get(1) prm%q = src%get_asFloat('q') prm%gamma_crit = src%get_asFloat('gamma_crit') @@ -77,10 +68,10 @@ module function isoductile_init(source_length) result(mySources) if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit' - Nconstituents=count(material_phaseAt==p) * discretization_nIPs - call phase_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0) - damageState(p)%p(sourceOffset)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal) - if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' + Nconstituents=count(material_phaseAt2==p) + call phase_allocateState(damageState(p),Nconstituents,1,1,0) + damageState(p)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' end associate @@ -88,7 +79,6 @@ module function isoductile_init(source_length) result(mySources) ! exit if any parameter is out of range if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isoDuctile)') endif - enddo enddo @@ -98,29 +88,16 @@ end function isoductile_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -module subroutine isoductile_dotState(co, ip, el) +module subroutine isoductile_dotState(ph, me) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element - - integer :: & ph, & - me, & - sourceOffset, & - damageOffset, & - homog + me - ph = material_phaseAt(co,el) - me = material_phasememberAt(co,ip,el) - sourceOffset = source_damage_isoDuctile_offset(ph) - homog = material_homogenizationAt(el) - damageOffset = material_homogenizationMemberAt(ip,el) - associate(prm => param(source_damage_isoDuctile_instance(ph))) - damageState(ph)%p(sourceOffset)%dotState(1,me) = & - sum(plasticState(ph)%slipRate(:,me))/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit + associate(prm => param(ph)) + damageState(ph)%dotState(1,me) = sum(plasticState(ph)%slipRate(:,me)) & + / (prm%gamma_crit*damage_phi(ph,me)**prm%q) end associate end subroutine isoductile_dotState @@ -129,28 +106,24 @@ end subroutine isoductile_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- -module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) +module subroutine isoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) integer, intent(in) :: & - phase, & - constituent + ph, & + me real(pReal), intent(in) :: & phi real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi - integer :: & - sourceOffset - sourceOffset = source_damage_isoDuctile_offset(phase) - - dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent) + dLocalphiDot_dPhi = -damageState(ph)%state(1,me) localphiDot = 1.0_pReal & + dLocalphiDot_dPhi*phi -end subroutine source_damage_isoDuctile_getRateAndItsTangent +end subroutine isoductile_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- @@ -163,14 +136,13 @@ module subroutine isoductile_results(phase,group) integer :: o - associate(prm => param(source_damage_isoDuctile_instance(phase)), & - stt => damageState(phase)%p(source_damage_isoDuctile_offset(phase))%state) - outputsLoop: do o = 1,size(prm%output) - select case(trim(prm%output(o))) - case ('f_phi') - call results_writeDataset(group,stt,trim(prm%output(o)),'driving force','J/m³') - end select - enddo outputsLoop + associate(prm => param(phase), stt => damageState(phase)%state) + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + case ('f_phi') + call results_writeDataset(group,stt,trim(prm%output(o)),'driving force','J/m³') + end select + enddo outputsLoop end associate end subroutine isoductile_results diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index efd4d51d2..c10a916da 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -1,7 +1,8 @@ !---------------------------------------------------------------------------------------------------- !> @brief internal microstructure state for all plasticity constitutive models !---------------------------------------------------------------------------------------------------- -submodule(phase) mechanics +submodule(phase) mechanical + enum, bind(c); enumerator :: & ELASTICITY_UNDEFINED_ID, & @@ -22,8 +23,6 @@ submodule(phase) mechanics KINEMATICS_THERMAL_EXPANSION_ID end enum - integer(kind(KINEMATICS_UNDEFINED_ID)), dimension(:,:), allocatable :: & - phase_kinematics integer(kind(ELASTICITY_UNDEFINED_ID)), dimension(:), allocatable :: & phase_elasticity !< elasticity of each phase integer(kind(STIFFNESS_DEGRADATION_UNDEFINED_ID)), dimension(:,:), allocatable :: & @@ -98,12 +97,9 @@ submodule(phase) mechanics end function plastic_deltaState module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & - S, Fi, co, ip, el) - + S, Fi, ph,me) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element + ph,me real(pReal), intent(in), dimension(3,3) :: & S !< 2nd Piola-Kirchhoff stress real(pReal), intent(in), dimension(3,3) :: & @@ -206,7 +202,7 @@ module subroutine mechanical_init(phases) elastic, & stiffDegradation - print'(/,a)', ' <<<+- phase:mechanics init -+>>>' + print'(/,a)', ' <<<+- phase:mechanical init -+>>>' !------------------------------------------------------------------------------------------------- ! initialize elasticity (hooke) !ToDO: Maybe move to elastic submodule along with function homogenizedC? @@ -381,7 +377,7 @@ subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & DegradationLoop: do d = 1, phase_NstiffnessDegradations(material_phaseAt(co,el)) degradationType: select case(phase_stiffnessDegradation(d,material_phaseAt(co,el))) case (STIFFNESS_DEGRADATION_damage_ID) degradationType - C = C * damage(ho)%p(material_homogenizationMemberAt(ip,el))**2 + C = C * phase_damage_get_phi(co,ip,el)**2 end select degradationType enddo DegradationLoop @@ -583,7 +579,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken) enddo LpLoop call phase_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, & - S, Fi_new, co, ip, el) + S, Fi_new, ph,me) !* update current residuum and check for convergence of loop atol_Li = max(num%rtol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error @@ -1141,7 +1137,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) real(pReal) :: & formerSubStep integer :: & - so, ph, me, sizeDotState + ph, me, sizeDotState logical :: todo real(pReal) :: subFrac,subStep real(pReal), dimension(3,3) :: & @@ -1162,10 +1158,9 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,me) subState0 = plasticState(ph)%State0(:,me) + if (damageState(ph)%sizeState > 0) & + damageState(ph)%subState0(:,me) = damageState(ph)%state0(:,me) - do so = 1, phase_Nsources(ph) - damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state0(:,me) - enddo subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,me) subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,me) subF0 = phase_mechanical_F0(ph)%data(1:3,1:3,me) @@ -1191,9 +1186,9 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) subFp0 = phase_mechanical_Fp(ph)%data(1:3,1:3,me) subFi0 = phase_mechanical_Fi(ph)%data(1:3,1:3,me) subState0 = plasticState(ph)%state(:,me) - do so = 1, phase_Nsources(ph) - damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state(:,me) - enddo + if (damageState(ph)%sizeState > 0) & + damageState(ph)%subState0(:,me) = damageState(ph)%state(:,me) + endif !-------------------------------------------------------------------------------------------------- ! cut back (reduced time and restore) @@ -1207,9 +1202,8 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) phase_mechanical_Li(ph)%data(1:3,1:3,me) = subLi0 endif plasticState(ph)%state(:,me) = subState0 - do so = 1, phase_Nsources(ph) - damageState(ph)%p(so)%state(:,me) = damageState(ph)%p(so)%subState0(:,me) - enddo + if (damageState(ph)%sizeState > 0) & + damageState(ph)%state(:,me) = damageState(ph)%subState0(:,me) todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair) endif @@ -1303,7 +1297,7 @@ module function phase_mechanical_dPdF(dt,co,ip,el) result(dPdF) call phase_LiAndItsTangents(devNull,dLidS,dLidFi, & phase_mechanical_S(ph)%data(1:3,1:3,me), & phase_mechanical_Fi(ph)%data(1:3,1:3,me), & - co,ip,el) + ph,me) invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,me)) invFi = math_inv33(phase_mechanical_Fi(ph)%data(1:3,1:3,me)) @@ -1505,4 +1499,4 @@ module subroutine phase_mechanical_setF(F,co,ip,el) end subroutine phase_mechanical_setF -end submodule mechanics +end submodule mechanical diff --git a/src/phase_mechanical_eigen.f90 b/src/phase_mechanical_eigen.f90 index 6df993565..bee7dc2d2 100644 --- a/src/phase_mechanical_eigen.f90 +++ b/src/phase_mechanical_eigen.f90 @@ -1,50 +1,29 @@ -submodule(phase:mechanics) eigendeformation +submodule(phase:mechanical) eigen + + integer, dimension(:), allocatable :: & + Nmodels + + integer(kind(KINEMATICS_UNDEFINED_ID)), dimension(:,:), allocatable :: & + model + integer(kind(KINEMATICS_UNDEFINED_ID)), dimension(:), allocatable :: & + model_damage interface - module function kinematics_cleavage_opening_init(kinematics_length) result(myKinematics) - integer, intent(in) :: kinematics_length - logical, dimension(:,:), allocatable :: myKinematics + module function kinematics_cleavage_opening_init() result(myKinematics) + logical, dimension(:), allocatable :: myKinematics end function kinematics_cleavage_opening_init - module function kinematics_slipplane_opening_init(kinematics_length) result(myKinematics) - integer, intent(in) :: kinematics_length - logical, dimension(:,:), allocatable :: myKinematics + module function kinematics_slipplane_opening_init() result(myKinematics) + logical, dimension(:), allocatable :: myKinematics end function kinematics_slipplane_opening_init - module function kinematics_thermal_expansion_init(kinematics_length) result(myKinematics) + module function thermalexpansion_init(kinematics_length) result(myKinematics) integer, intent(in) :: kinematics_length logical, dimension(:,:), allocatable :: myKinematics - end function kinematics_thermal_expansion_init - - module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el) - integer, intent(in) :: & - co, & !< grain number - ip, & !< integration point number - el !< element number - real(pReal), intent(in), dimension(3,3) :: & - S - real(pReal), intent(out), dimension(3,3) :: & - Ld !< damage velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) - end subroutine kinematics_cleavage_opening_LiAndItsTangent - - module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el) - integer, intent(in) :: & - co, & !< grain number - ip, & !< integration point number - el !< element number - real(pReal), intent(in), dimension(3,3) :: & - S - real(pReal), intent(out), dimension(3,3) :: & - Ld !< damage velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) - end subroutine kinematics_slipplane_opening_LiAndItsTangent + end function thermalexpansion_init module subroutine thermalexpansion_LiAndItsTangent(Li, dLi_dTstar, ph,me) integer, intent(in) :: ph, me - !< element number real(pReal), intent(out), dimension(3,3) :: & Li !< thermal velocity gradient real(pReal), intent(out), dimension(3,3,3,3) :: & @@ -66,29 +45,36 @@ module subroutine eigendeformation_init(phases) ph class(tNode), pointer :: & phase, & - kinematics + kinematics, & + damage, & + mechanics - print'(/,a)', ' <<<+- phase:mechanics:eigendeformation init -+>>>' + print'(/,a)', ' <<<+- phase:mechanical:eigen init -+>>>' !-------------------------------------------------------------------------------------------------- -! initialize kinematic mechanisms - allocate(phase_Nkinematics(phases%length),source = 0) +! explicit eigen mechanisms + allocate(Nmodels(phases%length),source = 0) + do ph = 1,phases%length phase => phases%get(ph) - kinematics => phase%get('kinematics',defaultVal=emptyList) - phase_Nkinematics(ph) = kinematics%length + mechanics => phase%get('mechanics') + kinematics => mechanics%get('eigen',defaultVal=emptyList) + Nmodels(ph) = kinematics%length enddo - allocate(phase_kinematics(maxval(phase_Nkinematics),phases%length), source = KINEMATICS_undefined_ID) + allocate(model(maxval(Nmodels),phases%length), source = KINEMATICS_undefined_ID) - if(maxval(phase_Nkinematics) /= 0) then - where(kinematics_cleavage_opening_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_cleavage_opening_ID - where(kinematics_slipplane_opening_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_slipplane_opening_ID - where(kinematics_thermal_expansion_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_thermal_expansion_ID + if(maxval(Nmodels) /= 0) then + where(thermalexpansion_init(maxval(Nmodels))) model = KINEMATICS_thermal_expansion_ID endif -end subroutine eigendeformation_init + allocate(model_damage(phases%length), source = KINEMATICS_UNDEFINED_ID) + where(kinematics_cleavage_opening_init()) model_damage = KINEMATICS_cleavage_opening_ID + where(kinematics_slipplane_opening_init()) model_damage = KINEMATICS_slipplane_opening_ID + + +end subroutine eigendeformation_init !-------------------------------------------------------------------------------------------------- @@ -104,17 +90,19 @@ function kinematics_active(kinematics_label,kinematics_length) result(active_ki phases, & phase, & kinematics, & - kinematics_type + kinematics_type, & + mechanics integer :: p,k phases => config_material%get('phase') allocate(active_kinematics(kinematics_length,phases%length), source = .false. ) do p = 1, phases%length phase => phases%get(p) - kinematics => phase%get('kinematics',defaultVal=emptyList) + mechanics => phase%get('mechanics') + kinematics => mechanics%get('eigen',defaultVal=emptyList) do k = 1, kinematics%length kinematics_type => kinematics%get(k) - if(kinematics_type%get_asString('type') == kinematics_label) active_kinematics(k,p) = .true. + active_kinematics(k,p) = kinematics_type%get_asString('type') == kinematics_label enddo enddo @@ -122,17 +110,46 @@ function kinematics_active(kinematics_label,kinematics_length) result(active_ki end function kinematics_active + +!-------------------------------------------------------------------------------------------------- +!> @brief checks if a kinematic mechanism is active or not +!-------------------------------------------------------------------------------------------------- +function kinematics_active2(kinematics_label) result(active_kinematics) + + character(len=*), intent(in) :: kinematics_label !< name of kinematic mechanism + logical, dimension(:), allocatable :: active_kinematics + + class(tNode), pointer :: & + phases, & + phase, & + kinematics, & + kinematics_type + integer :: p + + phases => config_material%get('phase') + allocate(active_kinematics(phases%length), source = .false. ) + do p = 1, phases%length + phase => phases%get(p) + kinematics => phase%get('damage',defaultVal=emptyList) + if(kinematics%length < 1) return + kinematics_type => kinematics%get(1) + if (.not. kinematics_type%contains('type')) continue + active_kinematics(p) = kinematics_type%get_asString('type',defaultVal='n/a') == kinematics_label + enddo + + +end function kinematics_active2 + + !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient ! ToDo: MD: S is Mi? !-------------------------------------------------------------------------------------------------- module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & - S, Fi, co, ip, el) + S, Fi, ph,me) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element + ph,me real(pReal), intent(in), dimension(3,3) :: & S !< 2nd Piola-Kirchhoff stress real(pReal), intent(in), dimension(3,3) :: & @@ -152,44 +169,49 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & real(pReal) :: & detFi integer :: & - k, i, j, & - instance, of, me, ph + k, i, j + logical :: active + active = .false. Li = 0.0_pReal dLi_dS = 0.0_pReal dLi_dFi = 0.0_pReal - plasticType: select case (phase_plasticity(material_phaseAt(co,el))) + + plasticType: select case (phase_plasticity(ph)) case (PLASTICITY_isotropic_ID) plasticType - of = material_phasememberAt(co,ip,el) - instance = phase_plasticInstance(material_phaseAt(co,el)) - call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,instance,of) - case default plasticType - my_Li = 0.0_pReal - my_dLi_dS = 0.0_pReal + call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,phase_plasticInstance(ph),me) + Li = Li + my_Li + dLi_dS = dLi_dS + my_dLi_dS + active = .true. end select plasticType - Li = Li + my_Li - dLi_dS = dLi_dS + my_dLi_dS - KinematicsLoop: do k = 1, phase_Nkinematics(material_phaseAt(co,el)) - kinematicsType: select case (phase_kinematics(k,material_phaseAt(co,el))) - case (KINEMATICS_cleavage_opening_ID) kinematicsType - call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) - case (KINEMATICS_slipplane_opening_ID) kinematicsType - call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) + KinematicsLoop: do k = 1, Nmodels(ph) + kinematicsType: select case (model(k,ph)) case (KINEMATICS_thermal_expansion_ID) kinematicsType - me = material_phaseMemberAt(co,ip,el) - ph = material_phaseAt(co,el) call thermalexpansion_LiAndItsTangent(my_Li, my_dLi_dS, ph,me) - case default kinematicsType - my_Li = 0.0_pReal - my_dLi_dS = 0.0_pReal + Li = Li + my_Li + dLi_dS = dLi_dS + my_dLi_dS + active = .true. end select kinematicsType - Li = Li + my_Li - dLi_dS = dLi_dS + my_dLi_dS enddo KinematicsLoop + select case (model_damage(ph)) + case (KINEMATICS_cleavage_opening_ID) + call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me) + Li = Li + my_Li + dLi_dS = dLi_dS + my_dLi_dS + active = .true. + case (KINEMATICS_slipplane_opening_ID) + call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me) + Li = Li + my_Li + dLi_dS = dLi_dS + my_dLi_dS + active = .true. + end select + + if(.not. active) return + FiInv = math_inv33(Fi) detFi = math_det33(Fi) Li = matmul(matmul(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration @@ -204,4 +226,4 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & end subroutine phase_LiAndItsTangents -end submodule eigendeformation +end submodule eigen diff --git a/src/phase_mechanical_eigen_cleavageopening.f90 b/src/phase_mechanical_eigen_cleavageopening.f90 index 59be7837a..0f48be1a4 100644 --- a/src/phase_mechanical_eigen_cleavageopening.f90 +++ b/src/phase_mechanical_eigen_cleavageopening.f90 @@ -4,24 +4,7 @@ !> @brief material subroutine incorporating kinematics resulting from opening of cleavage planes !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule(phase:eigendeformation) cleavageopening - - integer, dimension(:), allocatable :: kinematics_cleavage_opening_instance - - type :: tParameters !< container type for internal constitutive parameters - integer :: & - sum_N_cl !< total number of cleavage planes - real(pReal) :: & - dot_o, & !< opening rate of cleavage planes - q !< damage rate sensitivity - real(pReal), dimension(:), allocatable :: & - g_crit - real(pReal), dimension(:,:,:,:), allocatable :: & - cleavage_systems - end type tParameters - - type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances) - +submodule(phase:eigen) cleavageopening contains @@ -30,136 +13,20 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function kinematics_cleavage_opening_init(kinematics_length) result(myKinematics) +module function kinematics_cleavage_opening_init() result(myKinematics) - integer, intent(in) :: kinematics_length - logical, dimension(:,:), allocatable :: myKinematics + logical, dimension(:), allocatable :: myKinematics - integer :: Ninstances,p,k - integer, dimension(:), allocatable :: N_cl !< active number of cleavage systems per family - character(len=pStringLen) :: extmsg = '' - class(tNode), pointer :: & - phases, & - phase, & - kinematics, & - kinematic_type - print'(/,a)', ' <<<+- phase:mechanics:eigendeformation:cleavageopening init -+>>>' - - myKinematics = kinematics_active('cleavage_opening',kinematics_length) - Ninstances = count(myKinematics) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) - if(Ninstances == 0) return - - phases => config_material%get('phase') - allocate(param(Ninstances)) - allocate(kinematics_cleavage_opening_instance(phases%length), source=0) - - do p = 1, phases%length - if(any(myKinematics(:,p))) kinematics_cleavage_opening_instance(p) = count(myKinematics(:,1:p)) - phase => phases%get(p) - if(count(myKinematics(:,p)) == 0) cycle - kinematics => phase%get('kinematics') - do k = 1, kinematics%length - if(myKinematics(k,p)) then - associate(prm => param(kinematics_cleavage_opening_instance(p))) - kinematic_type => kinematics%get(k) - - N_cl = kinematic_type%get_asInts('N_cl') - prm%sum_N_cl = sum(abs(N_cl)) - - prm%q = kinematic_type%get_asFloat('q') - prm%dot_o = kinematic_type%get_asFloat('dot_o') - - prm%g_crit = kinematic_type%get_asFloats('g_crit',requiredSize=size(N_cl)) - - prm%cleavage_systems = lattice_SchmidMatrix_cleavage(N_cl,phase%get_asString('lattice'),& - phase%get_asFloat('c/a',defaultVal=0.0_pReal)) - - ! expand: family => system - prm%g_crit = math_expand(prm%g_crit,N_cl) - - ! sanity checks - if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' - if (prm%dot_o <= 0.0_pReal) extmsg = trim(extmsg)//' dot_o' - if (any(prm%g_crit < 0.0_pReal)) extmsg = trim(extmsg)//' g_crit' - -!-------------------------------------------------------------------------------------------------- -! exit if any parameter is out of range - if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(cleavage_opening)') - end associate - endif - enddo - enddo + myKinematics = kinematics_active2('anisobrittle') + if(count(myKinematics) == 0) return + print'(/,a)', ' <<<+- phase:mechanical:eigen:cleavageopening init -+>>>' + print'(a,i2)', ' # phases: ',count(myKinematics); flush(IO_STDOUT) end function kinematics_cleavage_opening_init -!-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the velocity gradient -!-------------------------------------------------------------------------------------------------- -module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el) - integer, intent(in) :: & - co, & !< grain number - ip, & !< integration point number - el !< element number - real(pReal), intent(in), dimension(3,3) :: & - S - real(pReal), intent(out), dimension(3,3) :: & - Ld !< damage velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) - - integer :: & - homog, damageOffset, & - i, k, l, m, n - real(pReal) :: & - traction_d, traction_t, traction_n, traction_crit, & - udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt - - homog = material_homogenizationAt(el) - damageOffset = material_homogenizationMemberAt(ip,el) - - Ld = 0.0_pReal - dLd_dTstar = 0.0_pReal - associate(prm => param(kinematics_cleavage_opening_instance(material_phaseAt(co,el)))) - do i = 1,prm%sum_N_cl - traction_crit = prm%g_crit(i)* damage(homog)%p(damageOffset)**2.0_pReal - - traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i)) - if (abs(traction_d) > traction_crit + tol_math_check) then - udotd = sign(1.0_pReal,traction_d)* prm%dot_o * ((abs(traction_d) - traction_crit)/traction_crit)**prm%q - Ld = Ld + udotd*prm%cleavage_systems(1:3,1:3,1,i) - dudotd_dt = sign(1.0_pReal,traction_d)*udotd*prm%q / (abs(traction_d) - traction_crit) - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) & - + dudotd_dt*prm%cleavage_systems(k,l,1,i) * prm%cleavage_systems(m,n,1,i) - endif - - traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i)) - if (abs(traction_t) > traction_crit + tol_math_check) then - udott = sign(1.0_pReal,traction_t)* prm%dot_o * ((abs(traction_t) - traction_crit)/traction_crit)**prm%q - Ld = Ld + udott*prm%cleavage_systems(1:3,1:3,2,i) - dudott_dt = sign(1.0_pReal,traction_t)*udott*prm%q / (abs(traction_t) - traction_crit) - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) & - + dudott_dt*prm%cleavage_systems(k,l,2,i) * prm%cleavage_systems(m,n,2,i) - endif - - traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i)) - if (abs(traction_n) > traction_crit + tol_math_check) then - udotn = sign(1.0_pReal,traction_n)* prm%dot_o * ((abs(traction_n) - traction_crit)/traction_crit)**prm%q - Ld = Ld + udotn*prm%cleavage_systems(1:3,1:3,3,i) - dudotn_dt = sign(1.0_pReal,traction_n)*udotn*prm%q / (abs(traction_n) - traction_crit) - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) & - + dudotn_dt*prm%cleavage_systems(k,l,3,i) * prm%cleavage_systems(m,n,3,i) - endif - enddo - end associate - -end subroutine kinematics_cleavage_opening_LiAndItsTangent end submodule cleavageopening diff --git a/src/phase_mechanical_eigen_slipplaneopening.f90 b/src/phase_mechanical_eigen_slipplaneopening.f90 index a144d39a6..e8a7d65b9 100644 --- a/src/phase_mechanical_eigen_slipplaneopening.f90 +++ b/src/phase_mechanical_eigen_slipplaneopening.f90 @@ -4,7 +4,7 @@ !> @brief material subroutine incorporating kinematics resulting from opening of slip planes !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule(phase:eigendeformation) slipplaneopening +submodule(phase:eigen) slipplaneopening integer, dimension(:), allocatable :: kinematics_slipplane_opening_instance @@ -32,12 +32,11 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function kinematics_slipplane_opening_init(kinematics_length) result(myKinematics) +module function kinematics_slipplane_opening_init() result(myKinematics) - integer, intent(in) :: kinematics_length - logical, dimension(:,:), allocatable :: myKinematics + logical, dimension(:), allocatable :: myKinematics - integer :: Ninstances,p,i,k + integer :: p,i character(len=pStringLen) :: extmsg = '' integer, dimension(:), allocatable :: N_sl real(pReal), dimension(:,:), allocatable :: d,n,t @@ -49,28 +48,26 @@ module function kinematics_slipplane_opening_init(kinematics_length) result(myKi kinematics, & kinematic_type - print'(/,a)', ' <<<+- phase:mechanics:eigendeformation:slipplaneopening init -+>>>' - myKinematics = kinematics_active('slipplane_opening',kinematics_length) - Ninstances = count(myKinematics) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) - if(Ninstances == 0) return + myKinematics = kinematics_active2('isoductile') + if(count(myKinematics) == 0) return + print'(/,a)', ' <<<+- phase:mechanical:eigen:slipplaneopening init -+>>>' + print'(a,i2)', ' # phases: ',count(myKinematics); flush(IO_STDOUT) + phases => config_material%get('phase') - allocate(kinematics_slipplane_opening_instance(phases%length), source=0) - allocate(param(Ninstances)) + allocate(param(phases%length)) do p = 1, phases%length - if(any(myKinematics(:,p))) kinematics_slipplane_opening_instance(p) = count(myKinematics(:,1:p)) - phase => phases%get(p) - mech => phase%get('mechanics') - pl => mech%get('plasticity') - if(count(myKinematics(:,p)) == 0) cycle - kinematics => phase%get('kinematics') - do k = 1, kinematics%length - if(myKinematics(k,p)) then - associate(prm => param(kinematics_slipplane_opening_instance(p))) - kinematic_type => kinematics%get(k) + if(myKinematics(p)) then + phase => phases%get(p) + mech => phase%get('mechanics') + pl => mech%get('plasticity') + + kinematics => phase%get('damage') + + associate(prm => param(p)) + kinematic_type => kinematics%get(1) prm%dot_o = kinematic_type%get_asFloat('dot_o') prm%q = kinematic_type%get_asFloat('q') @@ -105,9 +102,8 @@ module function kinematics_slipplane_opening_init(kinematics_length) result(myKi ! exit if any parameter is out of range if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(slipplane_opening)') - end associate - endif - enddo + end associate + endif enddo @@ -117,12 +113,10 @@ end function kinematics_slipplane_opening_init !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el) +module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) integer, intent(in) :: & - co, & !< grain number - ip, & !< integration point number - el !< element number + ph, me real(pReal), intent(in), dimension(3,3) :: & S real(pReal), intent(out), dimension(3,3) :: & @@ -131,19 +125,13 @@ module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) integer :: & - instance, phase, & - homog, damageOffset, & i, k, l, m, n real(pReal) :: & traction_d, traction_t, traction_n, traction_crit, & udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt - phase = material_phaseAt(co,el) - instance = kinematics_slipplane_opening_instance(phase) - homog = material_homogenizationAt(el) - damageOffset = material_homogenizationMemberAt(ip,el) - associate(prm => param(instance)) + associate(prm => param(ph)) Ld = 0.0_pReal dLd_dTstar = 0.0_pReal do i = 1, prm%sum_N_sl @@ -152,7 +140,7 @@ module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S traction_t = math_tensordot(S,prm%P_t(1:3,1:3,i)) traction_n = math_tensordot(S,prm%P_n(1:3,1:3,i)) - traction_crit = prm%g_crit(i)* damage(homog)%p(damageOffset) ! degrading critical load carrying capacity by damage + traction_crit = prm%g_crit(i)* damage_phi(ph,me) udotd = sign(1.0_pReal,traction_d)* prm%dot_o* ( abs(traction_d)/traction_crit & - abs(traction_d)/prm%g_crit(i))**prm%q diff --git a/src/phase_mechanical_eigen_thermalexpansion.f90 b/src/phase_mechanical_eigen_thermalexpansion.f90 index 6630f0eb7..86e7fa907 100644 --- a/src/phase_mechanical_eigen_thermalexpansion.f90 +++ b/src/phase_mechanical_eigen_thermalexpansion.f90 @@ -3,7 +3,7 @@ !> @brief material subroutine incorporating kinematics resulting from thermal expansion !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule(phase:eigendeformation) thermalexpansion +submodule(phase:eigen) thermalexpansion integer, dimension(:), allocatable :: kinematics_thermal_expansion_instance @@ -23,7 +23,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function kinematics_thermal_expansion_init(kinematics_length) result(myKinematics) +module function thermalexpansion_init(kinematics_length) result(myKinematics) integer, intent(in) :: kinematics_length logical, dimension(:,:), allocatable :: myKinematics @@ -36,7 +36,7 @@ module function kinematics_thermal_expansion_init(kinematics_length) result(myKi kinematics, & kinematic_type - print'(/,a)', ' <<<+- phase:mechanics:eigendeformation:thermalexpansion init -+>>>' + print'(/,a)', ' <<<+- phase:mechanical:eigen:thermalexpansion init -+>>>' myKinematics = kinematics_active('thermal_expansion',kinematics_length) Ninstances = count(myKinematics) @@ -77,7 +77,7 @@ module function kinematics_thermal_expansion_init(kinematics_length) result(myKi enddo -end function kinematics_thermal_expansion_init +end function thermalexpansion_init !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index ed1dc64fb..843273c72 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -1,4 +1,4 @@ -submodule(phase:mechanics) plastic +submodule(phase:mechanical) plastic interface @@ -226,7 +226,7 @@ contains module subroutine plastic_init - print'(/,a)', ' <<<+- phase:mechanics:plastic init -+>>>' + print'(/,a)', ' <<<+- phase:mechanical:plastic init -+>>>' where(plastic_none_init()) phase_plasticity = PLASTICITY_NONE_ID where(plastic_isotropic_init()) phase_plasticity = PLASTICITY_ISOTROPIC_ID diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index a9946b6e7..6eb1e4abe 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -97,13 +97,14 @@ module function plastic_dislotungsten_init() result(myPlasticity) mech, & pl - print'(/,a)', ' <<<+- phase:mechanics:plastic:dislotungsten init -+>>>' - myPlasticity = plastic_active('dislotungsten') Ninstances = count(myPlasticity) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) if(Ninstances == 0) return + print'(/,a)', ' <<<+- phase:mechanical:plastic:dislotungsten init -+>>>' + print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT) + + print*, 'Cereceda et al., International Journal of Plasticity 78:242–256, 2016' print*, 'https://dx.doi.org/10.1016/j.ijplas.2015.09.002' diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index e1259fbd2..ca385f417 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -144,13 +144,13 @@ module function plastic_dislotwin_init() result(myPlasticity) mech, & pl - print'(/,a)', ' <<<+- phase:mechanics:plastic:dislotwin init -+>>>' - myPlasticity = plastic_active('dislotwin') Ninstances = count(myPlasticity) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) if(Ninstances == 0) return + print'(/,a)', ' <<<+- phase:mechanical:plastic:dislotwin init -+>>>' + print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT) + print*, 'Ma and Roters, Acta Materialia 52(12):3603–3612, 2004' print*, 'https://doi.org/10.1016/j.actamat.2004.04.012'//IO_EOL diff --git a/src/phase_mechanical_plastic_isotropic.f90 b/src/phase_mechanical_plastic_isotropic.f90 index 1a4e1449a..e315a38fb 100644 --- a/src/phase_mechanical_plastic_isotropic.f90 +++ b/src/phase_mechanical_plastic_isotropic.f90 @@ -68,13 +68,14 @@ module function plastic_isotropic_init() result(myPlasticity) mech, & pl - print'(/,a)', ' <<<+- phase:mechanics:plastic:isotropic init -+>>>' - myPlasticity = plastic_active('isotropic') Ninstances = count(myPlasticity) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) if(Ninstances == 0) return + print'(/,a)', ' <<<+- phase:mechanical:plastic:isotropic init -+>>>' + print'(a,i0)', ' # phses: ',Ninstances; flush(IO_STDOUT) + + print*, 'Maiti and Eisenlohr, Scripta Materialia 145:37–40, 2018' print*, 'https://doi.org/10.1016/j.scriptamat.2017.09.047' diff --git a/src/phase_mechanical_plastic_kinehardening.f90 b/src/phase_mechanical_plastic_kinehardening.f90 index f00e4171b..1244988de 100644 --- a/src/phase_mechanical_plastic_kinehardening.f90 +++ b/src/phase_mechanical_plastic_kinehardening.f90 @@ -80,13 +80,14 @@ module function plastic_kinehardening_init() result(myPlasticity) mech, & pl - print'(/,a)', ' <<<+- phase:mechanics:plastic:kinehardening init -+>>>' - myPlasticity = plastic_active('kinehardening') Ninstances = count(myPlasticity) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) if(Ninstances == 0) return + print'(/,a)', ' <<<+- phase:mechanical:plastic:kinehardening init -+>>>' + print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT) + + allocate(param(Ninstances)) allocate(state(Ninstances)) allocate(dotState(Ninstances)) diff --git a/src/phase_mechanical_plastic_none.f90 b/src/phase_mechanical_plastic_none.f90 index b95eaa039..1510262cc 100644 --- a/src/phase_mechanical_plastic_none.f90 +++ b/src/phase_mechanical_plastic_none.f90 @@ -16,32 +16,20 @@ module function plastic_none_init() result(myPlasticity) logical, dimension(:), allocatable :: myPlasticity integer :: & - Ninstances, & p, & Nconstituents class(tNode), pointer :: & - phases, & - phase, & - mech, & - pl + phases - print'(/,a)', ' <<<+- phase:mechanics:plastic:none init -+>>>' + + myPlasticity = plastic_active('none') + if(count(myPlasticity) == 0) return + + print'(/,a)', ' <<<+- phase:mechanical:plastic:none init -+>>>' + print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) phases => config_material%get('phase') - allocate(myPlasticity(phases%length), source = .false.) do p = 1, phases%length - phase => phases%get(p) - mech => phase%get('mechanics') - pl => mech%get ('plasticity') - if(pl%get_asString('type') == 'none') myPlasticity(p) = .true. - enddo - - Ninstances = count(myPlasticity) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) - if(Ninstances == 0) return - - do p = 1, phases%length - phase => phases%get(p) if(.not. myPlasticity(p)) cycle Nconstituents = count(material_phaseAt2 == p) call phase_allocateState(plasticState(p),Nconstituents,0,0,0) diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 724086916..4e94f83ae 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -187,16 +187,16 @@ module function plastic_nonlocal_init() result(myPlasticity) mech, & pl - print'(/,a)', ' <<<+- phase:mechanics:plastic:nonlocal init -+>>>' - myPlasticity = plastic_active('nonlocal') Ninstances = count(myPlasticity) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) if(Ninstances == 0) then call geometry_plastic_nonlocal_disable return endif + print'(/,a)', ' <<<+- phase:mechanical:plastic:nonlocal init -+>>>' + print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT) + print*, 'Reuber et al., Acta Materialia 71:333–348, 2014' print*, 'https://doi.org/10.1016/j.actamat.2014.03.012'//IO_EOL diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90 index f9791faa6..bea7ef06b 100644 --- a/src/phase_mechanical_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanical_plastic_phenopowerlaw.f90 @@ -89,13 +89,15 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) mech, & pl - print'(/,a)', ' <<<+- phase:mechanics:plastic:phenopowerlaw init -+>>>' myPlasticity = plastic_active('phenopowerlaw') Ninstances = count(myPlasticity) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) if(Ninstances == 0) return + print'(/,a)', ' <<<+- phase:mechanical:plastic:phenopowerlaw init -+>>>' + print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT) + + allocate(param(Ninstances)) allocate(state(Ninstances)) allocate(dotState(Ninstances)) diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90 index 85ca2f7a5..012554aef 100644 --- a/src/phase_thermal.f90 +++ b/src/phase_thermal.f90 @@ -3,6 +3,12 @@ !---------------------------------------------------------------------------------------------------- submodule(phase) thermal + integer, dimension(:), allocatable :: & + thermal_Nsources + + type(tSourceState), allocatable, dimension(:) :: & + thermalState + enum, bind(c); enumerator :: & THERMAL_UNDEFINED_ID ,& THERMAL_DISSIPATION_ID, & @@ -33,8 +39,6 @@ submodule(phase) thermal end function externalheat_init - - module subroutine externalheat_dotState(ph, me) integer, intent(in) :: & ph, & diff --git a/src/phase_thermal_dissipation.f90 b/src/phase_thermal_dissipation.f90 index 3c0095401..9c75763dc 100644 --- a/src/phase_thermal_dissipation.f90 +++ b/src/phase_thermal_dissipation.f90 @@ -31,15 +31,14 @@ module function dissipation_init(source_length) result(mySources) phase, & sources, thermal, & src - integer :: Ninstances,so,Nconstituents,ph + integer :: so,Nconstituents,ph - print'(/,a)', ' <<<+- phase:thermal:dissipation init -+>>>' mySources = thermal_active('dissipation',source_length) + if(count(mySources) == 0) return + print'(/,a)', ' <<<+- phase:thermal:dissipation init -+>>>' + print'(a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT) - Ninstances = count(mySources) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) - if(Ninstances == 0) return phases => config_material%get('phase') allocate(param(phases%length)) diff --git a/src/phase_thermal_externalheat.f90 b/src/phase_thermal_externalheat.f90 index d2874ded0..c7f894215 100644 --- a/src/phase_thermal_externalheat.f90 +++ b/src/phase_thermal_externalheat.f90 @@ -8,7 +8,7 @@ submodule(phase:thermal) externalheat integer, dimension(:), allocatable :: & - source_thermal_externalheat_offset !< which source is my current thermal dissipation mechanism? + source_thermal_externalheat_offset !< which source is my current thermal dissipation mechanism? type :: tParameters !< container type for internal constitutive parameters real(pReal), dimension(:), allocatable :: & @@ -38,15 +38,14 @@ module function externalheat_init(source_length) result(mySources) phase, & sources, thermal, & src - integer :: Ninstances,so,Nconstituents,ph + integer :: so,Nconstituents,ph - print'(/,a)', ' <<<+- phase:thermal:externalheat init -+>>>' mySources = thermal_active('externalheat',source_length) + if(count(mySources) == 0) return + print'(/,a)', ' <<<+- phase:thermal:externalheat init -+>>>' + print'(a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT) - Ninstances = count(mySources) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) - if(Ninstances == 0) return phases => config_material%get('phase') allocate(param(phases%length))