diff --git a/src/IO.f90 b/src/IO.f90 index b6dbf5664..9ea35503b 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -604,7 +604,7 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg) msg = 'read only the first document' case default msg = 'unknown warning number' - end select + end select !$OMP CRITICAL (write2out) write(IO_STDERR,'(/,a)') ' ┌'//IO_DIVIDER//'┐' @@ -658,7 +658,7 @@ subroutine selfTest if(any([1,1,1] /= IO_stringPos('a'))) error stop 'IO_stringPos' if(any([2,2,3,5,5] /= IO_stringPos(' aa b'))) error stop 'IO_stringPos' - str=' 1.0 xxx' + str = ' 1.0 xxx' chunkPos = IO_stringPos(str) if(dNeq(1.0_pReal,IO_floatValue(str,chunkPos,1))) error stop 'IO_floatValue' diff --git a/src/YAML_types.f90 b/src/YAML_types.f90 index 3d83831e6..2a6bd64f9 100644 --- a/src/YAML_types.f90 +++ b/src/YAML_types.f90 @@ -266,7 +266,7 @@ subroutine selfTest if(any(dNeq(l2%get_as1dFloat(1),[2.0_pReal,3.0_pReal]))) error stop 'byIndex_as1dFloat' call l2%append(l3) x = l2%as2dFloat() - if(x(2,1)/= 4.0_pReal) error stop 'byKey_as2dFloat' + if(dNeq(x(2,1),4.0_pReal)) error stop 'byKey_as2dFloat' if(any(dNeq(pack(l2%as2dFloat(),.true.),& [2.0_pReal,4.0_pReal,3.0_pReal,5.0_pReal]))) error stop 'byKey_as2dFloat' n => l2 diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 6f99365f7..002a15708 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -196,7 +196,7 @@ function grid_damage_spectral_solution(timeinc) result(solution) ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 - call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),ce) + call homogenization_set_phi(phi_current(i,j,k),ce) enddo; enddo; enddo call VecMin(solution_vec,devNull,phi_min,ierr); CHKERRQ(ierr) @@ -233,7 +233,7 @@ subroutine grid_damage_spectral_forward(cutBack) call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 - call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),ce) + call homogenization_set_phi(phi_current(i,j,k),ce) enddo; enddo; enddo else phi_lastInc = phi_current @@ -259,7 +259,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) PetscObject :: dummy PetscErrorCode :: ierr integer :: i, j, k, ce - real(pReal) :: phiDot, dPhiDot_dPhi, mobility + real(pReal) :: phiDot, mobility phi_current = x_scal !-------------------------------------------------------------------------------------------------- @@ -281,7 +281,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 - call damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi_current(i,j,k),ce) + call damage_nonlocal_getSourceAndItsTangent(phiDot, phi_current(i,j,k),ce) mobility = damage_nonlocal_getMobility(ce) scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + phiDot) & + mobility*(phi_lastInc(i,j,k) - phi_current(i,j,k)) & diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index dd431ad9b..14e2affb9 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -282,9 +282,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) ce = ce + 1 call thermal_conduction_getSource(Tdot,1,ce) scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + Tdot) & - + thermal_conduction_getMassDensity (ce)* & - thermal_conduction_getSpecificHeat(ce)*(T_lastInc(i,j,k) - & - T_current(i,j,k))& + + homogenization_thermal_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) & + mu_ref*T_current(i,j,k) enddo; enddo; enddo @@ -314,7 +312,7 @@ subroutine updateReference do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 K_ref = K_ref + thermal_conduction_getConductivity(ce) - mu_ref = mu_ref + thermal_conduction_getMassDensity(ce)* thermal_conduction_getSpecificHeat(ce) + mu_ref = mu_ref + homogenization_thermal_mu_T(ce) enddo; enddo; enddo K_ref = K_ref*wgt call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 3e4e18f56..4ce75feca 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -39,8 +39,6 @@ module homogenization thermal_type !< thermal transport model integer(kind(DAMAGE_none_ID)), dimension(:), allocatable :: & damage_type !< nonlocal damage model - integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:), allocatable :: & - homogenization_type !< type of each homogenization type, private :: tNumerics_damage real(pReal) :: & @@ -117,6 +115,16 @@ module homogenization integer, intent(in) :: ho end subroutine mechanical_results + module subroutine damage_results(ho,group) + integer, intent(in) :: ho + character(len=*), intent(in) :: group + end subroutine damage_results + + module subroutine thermal_results(ho,group) + integer, intent(in) :: ho + character(len=*), intent(in) :: group + end subroutine thermal_results + module function mechanical_updateState(subdt,subF,ce) result(doneAndHappy) real(pReal), intent(in) :: & subdt !< current time step @@ -133,26 +141,16 @@ module homogenization real(pReal), dimension(3,3) :: K end function thermal_conduction_getConductivity - module function thermal_conduction_getSpecificHeat(ce) result(c_P) + module function homogenization_thermal_mu_T(ce) result(mu_T) integer, intent(in) :: ce - real(pReal) :: c_P - end function thermal_conduction_getSpecificHeat - - module function thermal_conduction_getMassDensity(ce) result(rho) - integer, intent(in) :: ce - real(pReal) :: rho - end function thermal_conduction_getMassDensity + real(pReal) :: mu_T + end function homogenization_thermal_mu_T module subroutine homogenization_thermal_setField(T,dot_T, ce) integer, intent(in) :: ce real(pReal), intent(in) :: T, dot_T end subroutine homogenization_thermal_setField - module subroutine thermal_conduction_results(ho,group) - integer, intent(in) :: ho - character(len=*), intent(in) :: group - end subroutine thermal_conduction_results - module function homogenization_thermal_T(ce) result(T) integer, intent(in) :: ce real(pReal) :: T @@ -170,37 +168,31 @@ module homogenization real(pReal) :: M end function damage_nonlocal_getMobility - module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ce) + module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, phi, ce) integer, intent(in) :: ce - real(pReal), intent(in) :: & + real(pReal), intent(in) :: & phi - real(pReal) :: & - phiDot, dPhiDot_dPhi + real(pReal), intent(out) :: & + phiDot end subroutine damage_nonlocal_getSourceAndItsTangent - module subroutine damage_nonlocal_putNonLocalDamage(phi,ce) + module subroutine homogenization_set_phi(phi,ce) integer, intent(in) :: ce real(pReal), intent(in) :: & phi - end subroutine damage_nonlocal_putNonLocalDamage - - module subroutine damage_nonlocal_results(ho,group) - integer, intent(in) :: ho - character(len=*), intent(in) :: group - end subroutine damage_nonlocal_results + end subroutine homogenization_set_phi end interface public :: & homogenization_init, & materialpoint_stressAndItsTangent, & - thermal_conduction_getSpecificHeat, & + homogenization_thermal_mu_T, & thermal_conduction_getConductivity, & - thermal_conduction_getMassDensity, & thermal_conduction_getSource, & damage_nonlocal_getMobility, & damage_nonlocal_getSourceAndItsTangent, & - damage_nonlocal_putNonLocalDamage, & + homogenization_set_phi, & homogenization_thermal_setfield, & homogenization_thermal_T, & homogenization_forward, & @@ -211,7 +203,6 @@ module homogenization DAMAGE_NONLOCAL_ID public :: & - damage_nonlocal_init, & damage_nonlocal_getDiffusion contains @@ -242,8 +233,6 @@ subroutine homogenization_init() call mechanical_init(num_homog) call thermal_init() call damage_init() - call damage_nonlocal_init() - end subroutine homogenization_init @@ -371,14 +360,14 @@ subroutine homogenization_results case(DAMAGE_NONLOCAL_ID) group = trim(group_base)//'/damage' call results_closeGroup(results_addGroup(group)) - call damage_nonlocal_results(ho,group) + call damage_results(ho,group) end select select case(thermal_type(ho)) case(THERMAL_CONDUCTION_ID) group = trim(group_base)//'/thermal' call results_closeGroup(results_addGroup(group)) - call thermal_conduction_results(ho,group) + call thermal_results(ho,group) end select enddo @@ -458,41 +447,6 @@ 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 !-------------------------------------------------------------------------------------------------- @@ -521,14 +475,12 @@ 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 @@ -536,23 +488,11 @@ subroutine material_parseHomogenization 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) do h=1, size(material_name_homogenization) homog => material_homogenization%get(h) - homogMech => homog%get('mechanical') - 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 if (homog%contains('thermal')) then homogThermal => homog%get('thermal') diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90 index 94978ccc8..a067dde43 100644 --- a/src/homogenization_damage.f90 +++ b/src/homogenization_damage.f90 @@ -37,12 +37,15 @@ module subroutine damage_init() class(tNode), pointer :: & configHomogenizations, & configHomogenization, & - configHomogenizationDamage + configHomogenizationDamage, & + num_generic, & + material_homogenization integer :: ho + integer :: Ninstances,Nmaterialpoints,h print'(/,a)', ' <<<+- homogenization:damage init -+>>>' - print'(/,a)', ' <<<+- homogenization:damage:isodamage init -+>>>' + print'(/,a)', ' <<<+- homogenization:damage:pass init -+>>>' configHomogenizations => config_material%get('homogenization') allocate(param(configHomogenizations%length)) @@ -65,6 +68,24 @@ module subroutine damage_init() end associate enddo +!------------------------------------------------------------------------------------ +! 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_init @@ -76,14 +97,10 @@ module subroutine damage_partition(ce) real(pReal) :: phi integer, intent(in) :: ce - integer :: co - if(damageState_h(material_homogenizationID(ce))%sizeState < 1) return phi = damagestate_h(material_homogenizationID(ce))%state(1,material_homogenizationEntry(ce)) - do co = 1, homogenization_Nconstituents(material_homogenizationID(ce)) - call phase_damage_set_phi(phi,co,ce) - enddo + call phase_damage_set_phi(phi,1,ce) end subroutine damage_partition @@ -95,17 +112,9 @@ end subroutine damage_partition module function damage_nonlocal_getMobility(ce) result(M) integer, intent(in) :: ce - integer :: & - co real(pReal) :: M - M = 0.0_pReal - - do co = 1, homogenization_Nconstituents(material_homogenizationID(ce)) - M = M + lattice_M(material_phaseID(co,ce)) - enddo - - M = M/real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal) + M = lattice_M(material_phaseID(1,ce)) end function damage_nonlocal_getMobility @@ -113,20 +122,15 @@ end function damage_nonlocal_getMobility !-------------------------------------------------------------------------------------------------- !> @brief calculates homogenized damage driving forces !-------------------------------------------------------------------------------------------------- -module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ce) +module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, phi, ce) integer, intent(in) :: ce - real(pReal), intent(in) :: & + real(pReal), intent(in) :: & phi - real(pReal) :: & - phiDot, dPhiDot_dPhi + real(pReal), intent(out) :: & + phiDot - phiDot = 0.0_pReal - dPhiDot_dPhi = 0.0_pReal - - call phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ce) - phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal) - dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal) + phiDot = phase_damage_phi_dot(phi, 1, ce) end subroutine damage_nonlocal_getSourceAndItsTangent @@ -134,7 +138,7 @@ end subroutine damage_nonlocal_getSourceAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief updated nonlocal damage field with solution from damage phase field PDE !-------------------------------------------------------------------------------------------------- -module subroutine damage_nonlocal_putNonLocalDamage(phi,ce) +module subroutine homogenization_set_phi(phi,ce) integer, intent(in) :: ce real(pReal), intent(in) :: & @@ -146,14 +150,15 @@ module subroutine damage_nonlocal_putNonLocalDamage(phi,ce) ho = material_homogenizationID(ce) en = material_homogenizationEntry(ce) damagestate_h(ho)%state(1,en) = phi + current(ho)%phi(en) = phi -end subroutine damage_nonlocal_putNonLocalDamage +end subroutine homogenization_set_phi !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -module subroutine damage_nonlocal_results(ho,group) +module subroutine damage_results(ho,group) integer, intent(in) :: ho character(len=*), intent(in) :: group @@ -170,6 +175,6 @@ module subroutine damage_nonlocal_results(ho,group) enddo outputsLoop end associate -end subroutine damage_nonlocal_results +end subroutine damage_results end submodule damage diff --git a/src/homogenization_mechanical.f90 b/src/homogenization_mechanical.f90 index 954aac403..7d1c64445 100644 --- a/src/homogenization_mechanical.f90 +++ b/src/homogenization_mechanical.f90 @@ -71,6 +71,9 @@ submodule(homogenization) mechanical end interface + integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:), allocatable :: & + homogenization_type !< type of each homogenization + contains !-------------------------------------------------------------------------------------------------- @@ -86,6 +89,8 @@ module subroutine mechanical_init(num_homog) print'(/,a)', ' <<<+- homogenization:mechanical init -+>>>' + call material_parseHomogenization2() + 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 homogenization_F = homogenization_F0 ! initialize to identity @@ -127,7 +132,7 @@ module subroutine mechanical_partition(subF,ce) end select chosenHomogenization do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) - call phase_mechanical_setF(Fs(1:3,1:3,co),co,ce) + call phase_set_F(Fs(1:3,1:3,co),co,ce) enddo @@ -150,13 +155,13 @@ module subroutine mechanical_homogenize(dt,ce) chosenHomogenization: select case(homogenization_type(material_homogenizationID(ce))) case (HOMOGENIZATION_NONE_ID) chosenHomogenization - homogenization_P(1:3,1:3,ce) = phase_mechanical_getP(1,ce) + homogenization_P(1:3,1:3,ce) = phase_P(1,ce) homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(dt,1,ce) case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization do co = 1, homogenization_Nconstituents(material_homogenizationID(ce)) dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ce) - Ps(:,:,co) = phase_mechanical_getP(co,ce) + Ps(:,:,co) = phase_P(co,ce) enddo call isostrain_averageStressAndItsTangent(& homogenization_P(1:3,1:3,ce), & @@ -167,7 +172,7 @@ module subroutine mechanical_homogenize(dt,ce) case (HOMOGENIZATION_RGC_ID) chosenHomogenization do co = 1, homogenization_Nconstituents(material_homogenizationID(ce)) dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ce) - Ps(:,:,co) = phase_mechanical_getP(co,ce) + Ps(:,:,co) = phase_P(co,ce) enddo call RGC_averageStressAndItsTangent(& homogenization_P(1:3,1:3,ce), & @@ -203,8 +208,8 @@ module function mechanical_updateState(subdt,subF,ce) result(doneAndHappy) if (homogenization_type(material_homogenizationID(ce)) == HOMOGENIZATION_RGC_ID) then do co = 1, homogenization_Nconstituents(material_homogenizationID(ce)) dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(subdt,co,ce) - Fs(:,:,co) = phase_mechanical_getF(co,ce) - Ps(:,:,co) = phase_mechanical_getP(co,ce) + Fs(:,:,co) = phase_F(co,ce) + Ps(:,:,co) = phase_P(co,ce) enddo doneAndHappy = RGC_updateState(Ps,Fs,subF,subdt,dPdFs,ce) else @@ -244,4 +249,38 @@ module subroutine mechanical_results(group_base,ho) end subroutine mechanical_results +!-------------------------------------------------------------------------------------------------- +!> @brief parses the homogenization part from the material configuration +!-------------------------------------------------------------------------------------------------- +subroutine material_parseHomogenization2() + + class(tNode), pointer :: & + material_homogenization, & + homog, & + homogMech + + integer :: h + + material_homogenization => config_material%get('homogenization') + + allocate(homogenization_type(size(material_name_homogenization)), source=HOMOGENIZATION_undefined_ID) + + do h=1, size(material_name_homogenization) + homog => material_homogenization%get(h) + homogMech => homog%get('mechanical') + 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 + enddo + +end subroutine material_parseHomogenization2 + + end submodule mechanical diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90 index 7bb81d24a..8dcddda7a 100644 --- a/src/homogenization_thermal.f90 +++ b/src/homogenization_thermal.f90 @@ -45,7 +45,7 @@ module subroutine thermal_init() print'(/,a)', ' <<<+- homogenization:thermal init -+>>>' - print'(/,a)', ' <<<+- homogenization:thermal:isotemperature init -+>>>' + print'(/,a)', ' <<<+- homogenization:thermal:pass init -+>>>' @@ -128,10 +128,20 @@ module function thermal_conduction_getConductivity(ce) result(K) end function thermal_conduction_getConductivity +module function homogenization_thermal_mu_T(ce) result(mu_T) + + integer, intent(in) :: ce + real(pReal) :: mu_T + + mu_T = c_P(ce) * rho(ce) + +end function homogenization_thermal_mu_T + + !-------------------------------------------------------------------------------------------------- !> @brief returns homogenized specific heat capacity !-------------------------------------------------------------------------------------------------- -module function thermal_conduction_getSpecificHeat(ce) result(c_P) +function c_P(ce) integer, intent(in) :: ce real(pReal) :: c_P @@ -139,21 +149,20 @@ module function thermal_conduction_getSpecificHeat(ce) result(c_P) integer :: co - c_P = 0.0_pReal - - do co = 1, homogenization_Nconstituents(material_homogenizationID(ce)) + c_P = lattice_c_p(material_phaseID(1,ce)) + do co = 2, homogenization_Nconstituents(material_homogenizationID(ce)) c_P = c_P + lattice_c_p(material_phaseID(co,ce)) enddo c_P = c_P / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal) -end function thermal_conduction_getSpecificHeat +end function c_P !-------------------------------------------------------------------------------------------------- !> @brief returns homogenized mass density !-------------------------------------------------------------------------------------------------- -module function thermal_conduction_getMassDensity(ce) result(rho) +function rho(ce) integer, intent(in) :: ce real(pReal) :: rho @@ -161,15 +170,14 @@ module function thermal_conduction_getMassDensity(ce) result(rho) integer :: co - rho = 0.0_pReal - - do co = 1, homogenization_Nconstituents(material_homogenizationID(ce)) + rho = lattice_rho(material_phaseID(1,ce)) + do co = 2, homogenization_Nconstituents(material_homogenizationID(ce)) rho = rho + lattice_rho(material_phaseID(co,ce)) enddo rho = rho / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal) -end function thermal_conduction_getMassDensity +end function rho @@ -193,7 +201,7 @@ end subroutine homogenization_thermal_setField !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -module subroutine thermal_conduction_results(ho,group) +module subroutine thermal_results(ho,group) integer, intent(in) :: ho character(len=*), intent(in) :: group @@ -209,7 +217,7 @@ module subroutine thermal_conduction_results(ho,group) enddo outputsLoop end associate -end subroutine thermal_conduction_results +end subroutine thermal_results module function homogenization_thermal_T(ce) result(T) diff --git a/src/phase.f90 b/src/phase.f90 index d10df84ba..1adba03d9 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -145,20 +145,20 @@ module phase real(pReal), dimension(3,3) :: L_p end function mechanical_L_p - module function phase_mechanical_getF(co,ce) result(F) + module function phase_F(co,ce) result(F) integer, intent(in) :: co, ce real(pReal), dimension(3,3) :: F - end function phase_mechanical_getF + end function phase_F module function mechanical_F_e(ph,me) result(F_e) integer, intent(in) :: ph,me real(pReal), dimension(3,3) :: F_e end function mechanical_F_e - module function phase_mechanical_getP(co,ce) result(P) + module function phase_P(co,ce) result(P) integer, intent(in) :: co, ce real(pReal), dimension(3,3) :: P - end function phase_mechanical_getP + end function phase_P module function phase_damage_get_phi(co,ip,el) result(phi) integer, intent(in) :: co, ip, el @@ -181,10 +181,10 @@ module phase end function damage_phi - module subroutine phase_mechanical_setF(F,co,ce) + module subroutine phase_set_F(F,co,ce) real(pReal), dimension(3,3), intent(in) :: F integer, intent(in) :: co, ce - end subroutine phase_mechanical_setF + end subroutine phase_set_F module subroutine phase_thermal_setField(T,dot_T, co,ce) real(pReal), intent(in) :: T, dot_T @@ -227,14 +227,13 @@ module phase end function phase_homogenizedC - module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ce) - integer, intent(in) :: ce + module function phase_damage_phi_dot(phi,co,ce) result(phi_dot) + integer, intent(in) :: ce,co real(pReal), intent(in) :: & phi !< damage parameter - real(pReal), intent(inout) :: & - phiDot, & - dPhiDot_dPhi - end subroutine phase_damage_getRateAndItsTangents + real(pReal) :: & + phi_dot + end function phase_damage_phi_dot module subroutine phase_thermal_getRate(TDot, ph,me) integer, intent(in) :: ph, me @@ -301,7 +300,7 @@ module phase public :: & phase_init, & phase_homogenizedC, & - phase_damage_getRateAndItsTangents, & + phase_damage_phi_dot, & phase_thermal_getRate, & phase_results, & phase_allocateState, & @@ -321,9 +320,9 @@ module phase phase_thermal_setField, & phase_damage_set_phi, & phase_damage_get_phi, & - phase_mechanical_getP, & - phase_mechanical_setF, & - phase_mechanical_getF + phase_P, & + phase_set_F, & + phase_F contains @@ -333,8 +332,7 @@ contains subroutine phase_init integer :: & - ph, & !< counter in phase loop - so !< counter in source loop + ph class (tNode), pointer :: & debug_constitutive, & materials, & @@ -476,7 +474,6 @@ subroutine crystallite_init() co, & !< counter in integration point component loop ip, & !< counter in integration point loop el, & !< counter in element loop - so, & cMax, & !< maximum number of integration point components iMax, & !< maximum number of integration points eMax !< maximum number of elements @@ -591,7 +588,7 @@ function crystallite_push33ToRef(co,ce, tensor33) ph = material_phaseID(co,ce) en = material_phaseEntry(co,ce) - T = matmul(material_orientation0(co,ph,en)%asMatrix(),transpose(math_inv33(phase_mechanical_getF(co,ce)))) ! ToDo: initial orientation correct? + T = matmul(material_orientation0(co,ph,en)%asMatrix(),transpose(math_inv33(phase_F(co,ce)))) ! ToDo: initial orientation correct? crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T)) diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 4397d5cad..62ce70368 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -65,43 +65,6 @@ submodule(phase) damage integer, intent(in) :: ph,me end subroutine isoductile_dotState - - 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 anisobrittle_getRateAndItsTangent - - 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 anisoductile_getRateAndItsTangent - - 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 isobrittle_getRateAndItsTangent - - 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 isoductile_getRateAndItsTangent - module subroutine anisobrittle_results(phase,group) integer, intent(in) :: phase character(len=*), intent(in) :: group @@ -179,53 +142,30 @@ end subroutine damage_init !---------------------------------------------------------------------------------------------- !< @brief returns local part of nonlocal damage driving force !---------------------------------------------------------------------------------------------- -module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ce) +module function phase_damage_phi_dot(phi,co,ce) result(phi_dot) - integer, intent(in) :: ce + integer, intent(in) :: ce,co real(pReal), intent(in) :: & phi !< damage parameter - real(pReal), intent(inout) :: & - phiDot, & - dPhiDot_dPhi - real(pReal) :: & - localphiDot, & - dLocalphiDot_dPhi + phi_dot + integer :: & ph, & - co, & - me + en - phiDot = 0.0_pReal - dPhiDot_dPhi = 0.0_pReal + ph = material_phaseID(co,ce) + en = material_phaseEntry(co,ce) - do co = 1, homogenization_Nconstituents(material_homogenizationID(ce)) - ph = material_phaseID(co,ce) - me = material_phaseEntry(co,ce) + select case(phase_source(ph)) + case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ISODUCTILE_ID,DAMAGE_ANISOBRITTLE_ID,DAMAGE_ANISODUCTILE_ID) + phi_dot = 1.0_pReal & + - phi*damageState(ph)%state(1,en) + case default + phi_dot = 0.0_pReal + end select - select case(phase_source(ph)) - case (DAMAGE_ISOBRITTLE_ID) - call isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, ph, me) - - case (DAMAGE_ISODUCTILE_ID) - call isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, ph, me) - - case (DAMAGE_ANISOBRITTLE_ID) - call anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) - - case (DAMAGE_ANISODUCTILE_ID) - call anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) - - case default - localphiDot = 0.0_pReal - dLocalphiDot_dPhi = 0.0_pReal - - end select - phiDot = phiDot + localphiDot - dPhiDot_dPhi = dPhiDot_dPhi + dLocalphiDot_dPhi - enddo - -end subroutine phase_damage_getRateAndItsTangents +end function phase_damage_phi_dot diff --git a/src/phase_damage_anisobrittle.f90 b/src/phase_damage_anisobrittle.f90 index e9672dd3b..ddcca9044 100644 --- a/src/phase_damage_anisobrittle.f90 +++ b/src/phase_damage_anisobrittle.f90 @@ -120,9 +120,6 @@ module subroutine anisobrittle_dotState(S, ph,me) S integer :: & - sourceOffset, & - damageOffset, & - homog, & i real(pReal) :: & traction_d, traction_t, traction_n, traction_crit @@ -148,29 +145,6 @@ module subroutine anisobrittle_dotState(S, ph,me) end subroutine anisobrittle_dotState -!-------------------------------------------------------------------------------------------------- -!> @brief returns local part of nonlocal damage driving force -!-------------------------------------------------------------------------------------------------- -module subroutine anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) - - integer, intent(in) :: & - ph, & - me - real(pReal), intent(in) :: & - phi - real(pReal), intent(out) :: & - localphiDot, & - dLocalphiDot_dPhi - - - dLocalphiDot_dPhi = -damageState(ph)%state(1,me) - - localphiDot = 1.0_pReal & - + dLocalphiDot_dPhi*phi - -end subroutine anisobrittle_getRateAndItsTangent - - !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_damage_anisoductile.f90 b/src/phase_damage_anisoductile.f90 index 54b63278f..049f148bc 100644 --- a/src/phase_damage_anisoductile.f90 +++ b/src/phase_damage_anisoductile.f90 @@ -113,29 +113,6 @@ module subroutine anisoductile_dotState(ph,me) end subroutine anisoductile_dotState -!-------------------------------------------------------------------------------------------------- -!> @brief returns local part of nonlocal damage driving force -!-------------------------------------------------------------------------------------------------- -module subroutine anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me) - - integer, intent(in) :: & - ph, & - me - real(pReal), intent(in) :: & - phi - real(pReal), intent(out) :: & - localphiDot, & - dLocalphiDot_dPhi - - - dLocalphiDot_dPhi = -damageState(ph)%state(1,me) - - localphiDot = 1.0_pReal & - + dLocalphiDot_dPhi*phi - -end subroutine anisoductile_getRateAndItsTangent - - !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_damage_isobrittle.f90 b/src/phase_damage_isobrittle.f90 index 59cedb554..e88a87352 100644 --- a/src/phase_damage_isobrittle.f90 +++ b/src/phase_damage_isobrittle.f90 @@ -113,29 +113,6 @@ module subroutine isobrittle_deltaState(C, Fe, ph,me) end subroutine isobrittle_deltaState -!-------------------------------------------------------------------------------------------------- -!> @brief returns local part of nonlocal damage driving force -!-------------------------------------------------------------------------------------------------- -module subroutine isobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) - - integer, intent(in) :: & - ph, me - real(pReal), intent(in) :: & - phi - real(pReal), intent(out) :: & - localphiDot, & - dLocalphiDot_dPhi - - - 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 isobrittle_getRateAndItsTangent - - !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_damage_isoductile.f90 b/src/phase_damage_isoductile.f90 index 1f1bba847..997b948fe 100644 --- a/src/phase_damage_isoductile.f90 +++ b/src/phase_damage_isoductile.f90 @@ -103,29 +103,6 @@ module subroutine isoductile_dotState(ph, me) end subroutine isoductile_dotState -!-------------------------------------------------------------------------------------------------- -!> @brief returns local part of nonlocal damage driving force -!-------------------------------------------------------------------------------------------------- -module subroutine isoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) - - integer, intent(in) :: & - ph, & - me - real(pReal), intent(in) :: & - phi - real(pReal), intent(out) :: & - localphiDot, & - dLocalphiDot_dPhi - - - dLocalphiDot_dPhi = -damageState(ph)%state(1,me) - - localphiDot = 1.0_pReal & - + dLocalphiDot_dPhi*phi - -end subroutine isoductile_getRateAndItsTangent - - !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index ddd6e1fcb..9f2bc4df2 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -1421,20 +1421,6 @@ module function mechanical_L_p(ph,me) result(L_p) end function mechanical_L_p -!---------------------------------------------------------------------------------------------- -!< @brief Get deformation gradient (for use by homogenization) -!---------------------------------------------------------------------------------------------- -module function phase_mechanical_getF(co,ce) result(F) - - integer, intent(in) :: co, ce - real(pReal), dimension(3,3) :: F - - - F = phase_mechanical_F(material_phaseID(co,ce))%data(1:3,1:3,material_phaseEntry(co,ce)) - -end function phase_mechanical_getF - - !---------------------------------------------------------------------------------------------- !< @brief Get elastic deformation gradient (for use by non-mech physics) !---------------------------------------------------------------------------------------------- @@ -1449,11 +1435,10 @@ module function mechanical_F_e(ph,me) result(F_e) end function mechanical_F_e - !---------------------------------------------------------------------------------------------- !< @brief Get second Piola-Kichhoff stress (for use by homogenization) !---------------------------------------------------------------------------------------------- -module function phase_mechanical_getP(co,ce) result(P) +module function phase_P(co,ce) result(P) integer, intent(in) :: co, ce real(pReal), dimension(3,3) :: P @@ -1461,11 +1446,27 @@ module function phase_mechanical_getP(co,ce) result(P) P = phase_mechanical_P(material_phaseID(co,ce))%data(1:3,1:3,material_phaseEntry(co,ce)) -end function phase_mechanical_getP +end function phase_P -! setter for homogenization -module subroutine phase_mechanical_setF(F,co,ce) +!---------------------------------------------------------------------------------------------- +!< @brief Get deformation gradient (for use by homogenization) +!---------------------------------------------------------------------------------------------- +module function phase_F(co,ce) result(F) + + integer, intent(in) :: co, ce + real(pReal), dimension(3,3) :: F + + + F = phase_mechanical_F(material_phaseID(co,ce))%data(1:3,1:3,material_phaseEntry(co,ce)) + +end function phase_F + + +!---------------------------------------------------------------------------------------------- +!< @brief Set deformation gradient (for use by homogenization) +!---------------------------------------------------------------------------------------------- +module subroutine phase_set_F(F,co,ce) real(pReal), dimension(3,3), intent(in) :: F integer, intent(in) :: co, ce @@ -1473,7 +1474,7 @@ module subroutine phase_mechanical_setF(F,co,ce) phase_mechanical_F(material_phaseID(co,ce))%data(1:3,1:3,material_phaseEntry(co,ce)) = F -end subroutine phase_mechanical_setF +end subroutine phase_set_F end submodule mechanical diff --git a/src/phase_mechanical_eigen.f90 b/src/phase_mechanical_eigen.f90 index 3f2aec58c..7c4e70c1c 100644 --- a/src/phase_mechanical_eigen.f90 +++ b/src/phase_mechanical_eigen.f90 @@ -46,7 +46,6 @@ module subroutine eigendeformation_init(phases) class(tNode), pointer :: & phase, & kinematics, & - damage, & mechanics print'(/,a)', ' <<<+- phase:mechanical:eigen init -+>>>'