diff --git a/PRIVATE b/PRIVATE index 313dd5de6..45ef93dbf 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 313dd5de618c996cdf9ace95a096f25e757386d9 +Subproject commit 45ef93dbfa3e0e6fa830914b3632e188c308a099 diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index beabfcae1..a5bbe69ca 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -46,10 +46,8 @@ #include "kinematics_slipplane_opening.f90" #include "crystallite.f90" #include "thermal_isothermal.f90" -#include "thermal_adiabatic.f90" #include "thermal_conduction.f90" #include "damage_none.f90" -#include "damage_local.f90" #include "damage_nonlocal.f90" #include "homogenization.f90" #include "homogenization_mech.f90" diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 31a6bde2d..e594ef5a1 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1629,7 +1629,6 @@ subroutine crystallite_forward enddo; enddo do i = 1,size(material_name_homogenization) homogState (i)%state0 = homogState (i)%state - thermalState(i)%state0 = thermalState(i)%state damageState (i)%state0 = damageState (i)%state enddo diff --git a/src/damage_local.f90 b/src/damage_local.f90 deleted file mode 100644 index 97eaf9a8c..000000000 --- a/src/damage_local.f90 +++ /dev/null @@ -1,172 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine for locally evolving damage field -!-------------------------------------------------------------------------------------------------- -module damage_local - use prec - use IO - use material - use config - use YAML_types - use constitutive - use results - - implicit none - private - - type :: tParameters - character(len=pStringLen), allocatable, dimension(:) :: & - output - end type tParameters - - type, private :: tNumerics - real(pReal) :: & - residualStiffness !< non-zero residual damage - end type tNumerics - - type(tparameters), dimension(:), allocatable :: & - param - - type(tNumerics), private :: num - - public :: & - damage_local_init, & - damage_local_updateState, & - damage_local_results - -contains - -!-------------------------------------------------------------------------------------------------- -!> @brief module initialization -!> @details reads in material parameters, allocates arrays, and does sanity checks -!-------------------------------------------------------------------------------------------------- -subroutine damage_local_init - - integer :: Ninstances,Nmaterialpoints,h - class(tNode), pointer :: & - num_generic, & - material_homogenization, & - homog, & - homogDamage - - print'(/,a)', ' <<<+- damage_local init -+>>>'; flush(IO_STDOUT) - -!---------------------------------------------------------------------------------------------- -! read numerics parameter and do sanity check - num_generic => config_numerics%get('generic',defaultVal=emptyDict) - num%residualStiffness = num_generic%get_asFloat('residualStiffness', defaultVal=1.0e-6_pReal) - if (num%residualStiffness < 0.0_pReal) call IO_error(301,ext_msg='residualStiffness') - - Ninstances = count(damage_type == DAMAGE_local_ID) - allocate(param(Ninstances)) - - material_homogenization => config_material%get('homogenization') - do h = 1, material_homogenization%length - if (damage_type(h) /= DAMAGE_LOCAL_ID) cycle - homog => material_homogenization%get(h) - homogDamage => homog%get('damage') - associate(prm => param(damage_typeInstance(h))) - -#if defined (__GFORTRAN__) - prm%output = output_asStrings(homogDamage) -#else - prm%output = homogDamage%get_asStrings('output',defaultVal=emptyStringArray) -#endif - - Nmaterialpoints = count(material_homogenizationAt == h) - damageState(h)%sizeState = 1 - allocate(damageState(h)%state0 (1,Nmaterialpoints), source=1.0_pReal) - allocate(damageState(h)%subState0(1,Nmaterialpoints), source=1.0_pReal) - allocate(damageState(h)%state (1,Nmaterialpoints), source=1.0_pReal) - - damage(h)%p => damageState(h)%state(1,:) - - end associate - enddo - -end subroutine damage_local_init - - -!-------------------------------------------------------------------------------------------------- -!> @brief calculates local change in damage field -!-------------------------------------------------------------------------------------------------- -function damage_local_updateState(subdt, ip, el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - subdt - logical, dimension(2) :: & - damage_local_updateState - integer :: & - homog, & - offset - real(pReal) :: & - phi, phiDot, dPhiDot_dPhi - - homog = material_homogenizationAt(el) - offset = material_homogenizationMemberAt(ip,el) - phi = damageState(homog)%subState0(1,offset) - call damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el) - phi = max(num%residualStiffness,min(1.0_pReal,phi + subdt*phiDot)) - - damage_local_updateState = [ abs(phi - damageState(homog)%state(1,offset)) & - <= 1.0e-2_pReal & - .or. abs(phi - damageState(homog)%state(1,offset)) & - <= 1.0e-6_pReal*abs(damageState(homog)%state(1,offset)), & - .true.] - - damageState(homog)%state(1,offset) = phi - -end function damage_local_updateState - - -!-------------------------------------------------------------------------------------------------- -!> @brief calculates homogenized local damage driving forces -!-------------------------------------------------------------------------------------------------- -subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - phi - real(pReal) :: & - phiDot, dPhiDot_dPhi - - phiDot = 0.0_pReal - dPhiDot_dPhi = 0.0_pReal - - call constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) - - phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) - dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) - -end subroutine damage_local_getSourceAndItsTangent - - -!-------------------------------------------------------------------------------------------------- -!> @brief writes results to HDF5 output file -!-------------------------------------------------------------------------------------------------- -subroutine damage_local_results(homog,group) - - integer, intent(in) :: homog - character(len=*), intent(in) :: group - - integer :: o - - associate(prm => param(damage_typeInstance(homog))) - outputsLoop: do o = 1,size(prm%output) - select case(prm%output(o)) - case ('phi') - call results_writeDataset(group,damage(homog)%p,prm%output(o),& - 'damage indicator','-') - end select - enddo outputsLoop - end associate - -end subroutine damage_local_results - - -end module damage_local diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 57478e039..30d9cfb90 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -15,10 +15,8 @@ module homogenization use FEsolving use discretization use thermal_isothermal - use thermal_adiabatic use thermal_conduction use damage_none - use damage_local use damage_nonlocal use results @@ -162,11 +160,9 @@ subroutine homogenization_init call mech_init(num_homog) if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init - if (any(thermal_type == THERMAL_adiabatic_ID)) call thermal_adiabatic_init if (any(thermal_type == THERMAL_conduction_ID)) call thermal_conduction_init if (any(damage_type == DAMAGE_none_ID)) call damage_none_init - if (any(damage_type == DAMAGE_local_ID)) call damage_local_init if (any(damage_type == DAMAGE_nonlocal_ID)) call damage_nonlocal_init @@ -212,10 +208,6 @@ subroutine materialpoint_stressAndItsTangent(dt) homogState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = & homogState(material_homogenizationAt(e))%State0( :,material_homogenizationMemberAt(i,e)) - if (thermalState(material_homogenizationAt(e))%sizeState > 0) & - thermalState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = & - thermalState(material_homogenizationAt(e))%State0( :,material_homogenizationMemberAt(i,e)) - if (damageState(material_homogenizationAt(e))%sizeState > 0) & damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = & damageState(material_homogenizationAt(e))%State0( :,material_homogenizationMemberAt(i,e)) @@ -245,9 +237,6 @@ subroutine materialpoint_stressAndItsTangent(dt) if(homogState(material_homogenizationAt(e))%sizeState > 0) & homogState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = & homogState(material_homogenizationAt(e))%State (:,material_homogenizationMemberAt(i,e)) - if(thermalState(material_homogenizationAt(e))%sizeState > 0) & - thermalState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = & - thermalState(material_homogenizationAt(e))%State (:,material_homogenizationMemberAt(i,e)) if(damageState(material_homogenizationAt(e))%sizeState > 0) & damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = & damageState(material_homogenizationAt(e))%State (:,material_homogenizationMemberAt(i,e)) @@ -270,9 +259,6 @@ subroutine materialpoint_stressAndItsTangent(dt) if(homogState(material_homogenizationAt(e))%sizeState > 0) & homogState(material_homogenizationAt(e))%State( :,material_homogenizationMemberAt(i,e)) = & homogState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) - if(thermalState(material_homogenizationAt(e))%sizeState > 0) & - thermalState(material_homogenizationAt(e))%State( :,material_homogenizationMemberAt(i,e)) = & - thermalState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) if(damageState(material_homogenizationAt(e))%sizeState > 0) & damageState(material_homogenizationAt(e))%State( :,material_homogenizationMemberAt(i,e)) = & damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) @@ -400,24 +386,6 @@ function updateState(subdt,subF,ip,el) el) end select chosenHomogenization - chosenThermal: select case (thermal_type(material_homogenizationAt(el))) - case (THERMAL_adiabatic_ID) chosenThermal - updateState = & - updateState .and. & - thermal_adiabatic_updateState(subdt, & - ip, & - el) - end select chosenThermal - - chosenDamage: select case (damage_type(material_homogenizationAt(el))) - case (DAMAGE_local_ID) chosenDamage - updateState = & - updateState .and. & - damage_local_updateState(subdt, & - ip, & - el) - end select chosenDamage - end function updateState @@ -441,8 +409,6 @@ subroutine homogenization_results group = trim(group_base)//'/damage' call results_closeGroup(results_addGroup(group)) select case(damage_type(p)) - case(DAMAGE_LOCAL_ID) - call damage_local_results(p,group) case(DAMAGE_NONLOCAL_ID) call damage_nonlocal_results(p,group) end select @@ -450,8 +416,6 @@ subroutine homogenization_results group = trim(group_base)//'/thermal' call results_closeGroup(results_addGroup(group)) select case(thermal_type(p)) - case(THERMAL_ADIABATIC_ID) - call thermal_adiabatic_results(p,group) case(THERMAL_CONDUCTION_ID) call thermal_conduction_results(p,group) end select diff --git a/src/material.f90 b/src/material.f90 index bb5f484f6..1f2437ad3 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -41,10 +41,8 @@ module material STIFFNESS_DEGRADATION_UNDEFINED_ID, & STIFFNESS_DEGRADATION_DAMAGE_ID, & THERMAL_ISOTHERMAL_ID, & - THERMAL_ADIABATIC_ID, & THERMAL_CONDUCTION_ID, & DAMAGE_NONE_ID, & - DAMAGE_LOCAL_ID, & DAMAGE_NONLOCAL_ID, & HOMOGENIZATION_UNDEFINED_ID, & HOMOGENIZATION_NONE_ID, & @@ -86,7 +84,6 @@ module material type(tState), allocatable, dimension(:), public :: & homogState, & - thermalState, & damageState type(Rotation), dimension(:,:,:), allocatable, public, protected :: & @@ -123,10 +120,8 @@ module material STIFFNESS_DEGRADATION_UNDEFINED_ID, & STIFFNESS_DEGRADATION_DAMAGE_ID, & THERMAL_ISOTHERMAL_ID, & - THERMAL_ADIABATIC_ID, & THERMAL_CONDUCTION_ID, & DAMAGE_NONE_ID, & - DAMAGE_LOCAL_ID, & DAMAGE_NONLOCAL_ID, & HOMOGENIZATION_NONE_ID, & HOMOGENIZATION_ISOSTRAIN_ID, & @@ -152,7 +147,6 @@ subroutine material_init(restart) allocate(homogState (size(material_name_homogenization))) - allocate(thermalState (size(material_name_homogenization))) allocate(damageState (size(material_name_homogenization))) allocate(temperature (size(material_name_homogenization))) @@ -218,8 +212,6 @@ subroutine material_parseHomogenization select case (homogThermal%get_asString('type')) case('isothermal') thermal_type(h) = THERMAL_isothermal_ID - case('adiabatic') - thermal_type(h) = THERMAL_adiabatic_ID case('conduction') thermal_type(h) = THERMAL_conduction_ID case default @@ -232,8 +224,6 @@ subroutine material_parseHomogenization select case (homogDamage%get_asString('type')) case('none') damage_type(h) = DAMAGE_none_ID - case('local') - damage_type(h) = DAMAGE_local_ID case('nonlocal') damage_type(h) = DAMAGE_nonlocal_ID case default diff --git a/src/thermal_adiabatic.f90 b/src/thermal_adiabatic.f90 deleted file mode 100644 index c67d004bf..000000000 --- a/src/thermal_adiabatic.f90 +++ /dev/null @@ -1,226 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine for adiabatic temperature evolution -!-------------------------------------------------------------------------------------------------- -module thermal_adiabatic - use prec - use config - use material - use results - use constitutive - use YAML_types - use crystallite - use lattice - - implicit none - private - - type :: tParameters - character(len=pStringLen), allocatable, dimension(:) :: & - output - end type tParameters - - type(tparameters), dimension(:), allocatable :: & - param - - public :: & - thermal_adiabatic_init, & - thermal_adiabatic_updateState, & - thermal_adiabatic_getSourceAndItsTangent, & - thermal_adiabatic_getSpecificHeat, & - thermal_adiabatic_getMassDensity, & - thermal_adiabatic_results - -contains - - -!-------------------------------------------------------------------------------------------------- -!> @brief module initialization -!> @details reads in material parameters, allocates arrays, and does sanity checks -!-------------------------------------------------------------------------------------------------- -subroutine thermal_adiabatic_init - - integer :: maxNinstances,h,Nmaterialpoints - class(tNode), pointer :: & - material_homogenization, & - homog, & - homogThermal - - print'(/,a)', ' <<<+- thermal_adiabatic init -+>>>'; flush(6) - - maxNinstances = count(thermal_type == THERMAL_adiabatic_ID) - if (maxNinstances == 0) return - - allocate(param(maxNinstances)) - - material_homogenization => config_material%get('homogenization') - do h = 1, size(material_name_homogenization) - if (thermal_type(h) /= THERMAL_adiabatic_ID) cycle - homog => material_homogenization%get(h) - homogThermal => homog%get('thermal') - - associate(prm => param(thermal_typeInstance(h))) - -#if defined (__GFORTRAN__) - prm%output = output_asStrings(homogThermal) -#else - prm%output = homogThermal%get_asStrings('output',defaultVal=emptyStringArray) -#endif - - Nmaterialpoints=count(material_homogenizationAt==h) - thermalState(h)%sizeState = 1 - allocate(thermalState(h)%state0 (1,Nmaterialpoints), source=thermal_initialT(h)) - allocate(thermalState(h)%subState0(1,Nmaterialpoints), source=thermal_initialT(h)) - allocate(thermalState(h)%state (1,Nmaterialpoints), source=thermal_initialT(h)) - - temperature(h)%p => thermalState(h)%state(1,:) - allocate(temperatureRate(h)%p(Nmaterialpoints),source = 0.0_pReal) - - end associate - enddo - -end subroutine thermal_adiabatic_init - - -!-------------------------------------------------------------------------------------------------- -!> @brief calculates adiabatic change in temperature based on local heat generation model -!-------------------------------------------------------------------------------------------------- -function thermal_adiabatic_updateState(subdt, ip, el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - subdt - - logical, dimension(2) :: & - thermal_adiabatic_updateState - integer :: & - homog, & - offset - real(pReal) :: & - T, Tdot, dTdot_dT - - homog = material_homogenizationAt(el) - offset = material_homogenizationMemberAt(ip,el) - - T = thermalState(homog)%subState0(1,offset) - call thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) - T = T + subdt*Tdot/(thermal_adiabatic_getSpecificHeat(ip,el)*thermal_adiabatic_getMassDensity(ip,el)) - - thermal_adiabatic_updateState = [ abs(T - thermalState(homog)%state(1,offset)) & - <= 1.0e-2_pReal & - .or. abs(T - thermalState(homog)%state(1,offset)) & - <= 1.0e-6_pReal*abs(thermalState(homog)%state(1,offset)), & - .true.] - - temperature (homog)%p(material_homogenizationMemberAt(ip,el)) = T - temperatureRate(homog)%p(material_homogenizationMemberAt(ip,el)) = & - (thermalState(homog)%state(1,offset) - thermalState(homog)%subState0(1,offset))/(subdt+tiny(0.0_pReal)) - -end function thermal_adiabatic_updateState - - -!-------------------------------------------------------------------------------------------------- -!> @brief returns heat generation rate -!-------------------------------------------------------------------------------------------------- -subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - T - real(pReal), intent(out) :: & - Tdot, dTdot_dT - integer :: & - homog - - Tdot = 0.0_pReal - dTdot_dT = 0.0_pReal - - homog = material_homogenizationAt(el) - call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, crystallite_S, crystallite_Lp, ip, el) - - Tdot = Tdot/real(homogenization_Nconstituents(homog),pReal) - dTdot_dT = dTdot_dT/real(homogenization_Nconstituents(homog),pReal) - -end subroutine thermal_adiabatic_getSourceAndItsTangent - - -!-------------------------------------------------------------------------------------------------- -!> @brief returns homogenized specific heat capacity -!-------------------------------------------------------------------------------------------------- -function thermal_adiabatic_getSpecificHeat(ip,el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - - real(pReal) :: & - thermal_adiabatic_getSpecificHeat - integer :: & - grain - - thermal_adiabatic_getSpecificHeat = 0.0_pReal - - do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - thermal_adiabatic_getSpecificHeat = thermal_adiabatic_getSpecificHeat & - + lattice_c_p(material_phaseAt(grain,el)) - enddo - - thermal_adiabatic_getSpecificHeat = thermal_adiabatic_getSpecificHeat & - / real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) - -end function thermal_adiabatic_getSpecificHeat - - -!-------------------------------------------------------------------------------------------------- -!> @brief returns homogenized mass density -!-------------------------------------------------------------------------------------------------- -function thermal_adiabatic_getMassDensity(ip,el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal) :: & - thermal_adiabatic_getMassDensity - integer :: & - grain - - thermal_adiabatic_getMassDensity = 0.0_pReal - - do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - thermal_adiabatic_getMassDensity = thermal_adiabatic_getMassDensity & - + lattice_rho(material_phaseAt(grain,el)) - enddo - - thermal_adiabatic_getMassDensity = thermal_adiabatic_getMassDensity & - / real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) - -end function thermal_adiabatic_getMassDensity - - -!-------------------------------------------------------------------------------------------------- -!> @brief writes results to HDF5 output file -!-------------------------------------------------------------------------------------------------- -subroutine thermal_adiabatic_results(homog,group) - - integer, intent(in) :: homog - character(len=*), intent(in) :: group - - integer :: o - - associate(prm => param(damage_typeInstance(homog))) - outputsLoop: do o = 1,size(prm%output) - select case(trim(prm%output(o))) - case('T') - call results_writeDataset(group,temperature(homog)%p,'T',& - 'temperature','K') - end select - enddo outputsLoop - end associate - -end subroutine thermal_adiabatic_results - -end module thermal_adiabatic diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index 602bdab35..37a407101 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -66,10 +66,6 @@ subroutine thermal_conduction_init #endif Nmaterialpoints=count(material_homogenizationAt==h) - thermalState(h)%sizeState = 0 - allocate(thermalState(h)%state0 (0,Nmaterialpoints)) - allocate(thermalState(h)%subState0(0,Nmaterialpoints)) - allocate(thermalState(h)%state (0,Nmaterialpoints)) allocate (temperature (h)%p(Nmaterialpoints), source=thermal_initialT(h)) allocate (temperatureRate(h)%p(Nmaterialpoints), source=0.0_pReal) diff --git a/src/thermal_isothermal.f90 b/src/thermal_isothermal.f90 index adf2257de..2a41ada49 100644 --- a/src/thermal_isothermal.f90 +++ b/src/thermal_isothermal.f90 @@ -25,10 +25,6 @@ subroutine thermal_isothermal_init if (thermal_type(h) /= THERMAL_isothermal_ID) cycle Nmaterialpoints = count(material_homogenizationAt == h) - thermalState(h)%sizeState = 0 - allocate(thermalState(h)%state0 (0,Nmaterialpoints)) - allocate(thermalState(h)%subState0(0,Nmaterialpoints)) - allocate(thermalState(h)%state (0,Nmaterialpoints)) allocate(temperature (h)%p(Nmaterialpoints),source=thermal_initialT(h)) allocate(temperatureRate(h)%p(Nmaterialpoints),source = 0.0_pReal)