diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index afbae027f..8d3f913fa 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -160,7 +160,7 @@ function grid_damage_spectral_solution(timeinc) result(solution) real(pReal), intent(in) :: & timeinc !< increment in time for current solution - integer :: i, j, k, cell + integer :: i, j, k, ce type(tSolutionState) :: solution PetscInt :: devNull PetscReal :: phi_min, phi_max, stagNorm, solnNorm @@ -193,10 +193,10 @@ function grid_damage_spectral_solution(timeinc) result(solution) !-------------------------------------------------------------------------------------------------- ! updating damage state - cell = 0 + ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 - call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),1,cell) + ce = ce + 1 + call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),ce) enddo; enddo; enddo call VecMin(solution_vec,devNull,phi_min,ierr); CHKERRQ(ierr) @@ -216,7 +216,7 @@ end function grid_damage_spectral_solution subroutine grid_damage_spectral_forward(cutBack) logical, intent(in) :: cutBack - integer :: i, j, k, cell + integer :: i, j, k, ce DM :: dm_local PetscScalar, dimension(:,:,:), pointer :: x_scal PetscErrorCode :: ierr @@ -226,14 +226,14 @@ subroutine grid_damage_spectral_forward(cutBack) phi_stagInc = phi_lastInc !-------------------------------------------------------------------------------------------------- ! reverting damage field state - cell = 0 + ce = 0 call SNESGetDM(damage_snes,dm_local,ierr); CHKERRQ(ierr) call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with x_scal(xstart:xend,ystart:yend,zstart:zend) = phi_current 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) - cell = cell + 1 - call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),1,cell) + ce = ce + 1 + call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),ce) enddo; enddo; enddo else phi_lastInc = phi_current @@ -258,7 +258,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) f_scal PetscObject :: dummy PetscErrorCode :: ierr - integer :: i, j, k, cell + integer :: i, j, k, ce real(pReal) :: phiDot, dPhiDot_dPhi, mobility phi_current = x_scal @@ -269,20 +269,20 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) call utilities_FFTscalarForward call utilities_fourierScalarGradient !< calculate gradient of damage field call utilities_FFTvectorBackward - cell = 0 + ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 - vectorField_real(1:3,i,j,k) = matmul(damage_nonlocal_getDiffusion(1,cell) - K_ref, & + ce = ce + 1 + vectorField_real(1:3,i,j,k) = matmul(damage_nonlocal_getDiffusion(ce) - K_ref, & vectorField_real(1:3,i,j,k)) enddo; enddo; enddo call utilities_FFTvectorForward call utilities_fourierVectorDivergence !< calculate damage divergence in fourier field call utilities_FFTscalarBackward - cell = 0 + ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 - call damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi_current(i,j,k), 1, cell) - mobility = damage_nonlocal_getMobility(1,cell) + ce = ce + 1 + call damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, 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)) & + mu_ref*phi_current(i,j,k) @@ -310,15 +310,15 @@ end subroutine formResidual !-------------------------------------------------------------------------------------------------- subroutine updateReference - integer :: i,j,k,cell,ierr + integer :: i,j,k,ce,ierr - cell = 0 + ce = 0 K_ref = 0.0_pReal mu_ref = 0.0_pReal do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 - K_ref = K_ref + damage_nonlocal_getDiffusion(1,cell) - mu_ref = mu_ref + damage_nonlocal_getMobility(1,cell) + ce = ce + 1 + K_ref = K_ref + damage_nonlocal_getDiffusion(ce) + mu_ref = mu_ref + damage_nonlocal_getMobility(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/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 866914f5b..6403ee955 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -268,7 +268,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 - vectorField_real(1:3,i,j,k) = matmul(thermal_conduction_getConductivity(1,ce) - K_ref, & + vectorField_real(1:3,i,j,k) = matmul(thermal_conduction_getConductivity(ce) - K_ref, & vectorField_real(1:3,i,j,k)) enddo; enddo; enddo call utilities_FFTvectorForward @@ -277,7 +277,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 thermal_conduction_getSource(Tdot, 1,ce) + 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) - & @@ -310,7 +310,7 @@ subroutine updateReference mu_ref = 0.0_pReal do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 - K_ref = K_ref + thermal_conduction_getConductivity(1,ce) + K_ref = K_ref + thermal_conduction_getConductivity(ce) mu_ref = mu_ref + thermal_conduction_getMassDensity(ce)* thermal_conduction_getSpecificHeat(ce) enddo; enddo; enddo K_ref = K_ref*wgt diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 49210806b..3c424585b 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -133,11 +133,9 @@ module homogenization end function mechanical_updateState - module function thermal_conduction_getConductivity(ip,el) result(K) + module function thermal_conduction_getConductivity(ce) result(K) - integer, intent(in) :: & - ip, & !< integration point number - el !< element number + integer, intent(in) :: ce real(pReal), dimension(3,3) :: K end function thermal_conduction_getConductivity @@ -167,25 +165,21 @@ module homogenization real(pReal) :: T end function homogenization_thermal_T - module subroutine thermal_conduction_getSource(Tdot, ip,el) + module subroutine thermal_conduction_getSource(Tdot, ip, el) integer, intent(in) :: & - ip, & !< integration point number - el !< element number + ip, & + el real(pReal), intent(out) :: Tdot end subroutine thermal_conduction_getSource - module function damage_nonlocal_getMobility(ip,el) result(M) - integer, intent(in) :: & - ip, & !< integration point number - el !< element number + module function damage_nonlocal_getMobility(ce) result(M) + integer, intent(in) :: ce real(pReal) :: M end function damage_nonlocal_getMobility - module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el) + module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ce) - integer, intent(in) :: & - ip, & !< integration point number - el !< element number + integer, intent(in) :: ce real(pReal), intent(in) :: & phi real(pReal) :: & @@ -193,11 +187,9 @@ module homogenization end subroutine damage_nonlocal_getSourceAndItsTangent - module subroutine damage_nonlocal_putNonLocalDamage(phi,ip,el) + module subroutine damage_nonlocal_putNonLocalDamage(phi,ce) - integer, intent(in) :: & - ip, & !< integration point number - el !< element number + integer, intent(in) :: ce real(pReal), intent(in) :: & phi @@ -521,28 +513,25 @@ end subroutine damage_nonlocal_init !-------------------------------------------------------------------------------------------------- !> @brief returns homogenized non local damage diffusion tensor in reference configuration !-------------------------------------------------------------------------------------------------- -function damage_nonlocal_getDiffusion(ip,el) +function damage_nonlocal_getDiffusion(ce) - integer, intent(in) :: & - ip, & !< integration point number - el !< element number + integer, intent(in) :: ce real(pReal), dimension(3,3) :: & damage_nonlocal_getDiffusion integer :: & - homog, & - grain, & - ce + ho, & + grain - homog = material_homogenizationAt(el) + ho = material_homogenizationAt2(ce) damage_nonlocal_getDiffusion = 0.0_pReal - ce = (el-1)*discretization_nIPs + ip - do grain = 1, homogenization_Nconstituents(homog) + + do grain = 1, homogenization_Nconstituents(ho) damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + & crystallite_push33ToRef(grain,ce,lattice_D(1:3,1:3,material_phaseAt2(grain,ce))) enddo damage_nonlocal_getDiffusion = & - num_damage%charLength**2*damage_nonlocal_getDiffusion/real(homogenization_Nconstituents(homog),pReal) + num_damage%charLength**2*damage_nonlocal_getDiffusion/real(homogenization_Nconstituents(ho),pReal) end function damage_nonlocal_getDiffusion diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90 index 128725839..993ea84f6 100644 --- a/src/homogenization_damage.f90 +++ b/src/homogenization_damage.f90 @@ -85,22 +85,20 @@ end subroutine damage_partition !-------------------------------------------------------------------------------------------------- !> @brief Returns homogenized nonlocal damage mobility !-------------------------------------------------------------------------------------------------- -module function damage_nonlocal_getMobility(ip,el) result(M) +module function damage_nonlocal_getMobility(ce) result(M) - integer, intent(in) :: & - ip, & !< integration point number - el !< element number + integer, intent(in) :: ce integer :: & co real(pReal) :: M M = 0.0_pReal - do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - M = M + lattice_M(material_phaseAt(co,el)) + do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) + M = M + lattice_M(material_phaseAt2(co,ce)) enddo - M = M/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) + M = M/real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) end function damage_nonlocal_getMobility @@ -108,11 +106,9 @@ end function damage_nonlocal_getMobility !-------------------------------------------------------------------------------------------------- !> @brief calculates homogenized damage driving forces !-------------------------------------------------------------------------------------------------- -module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el) +module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ce) - integer, intent(in) :: & - ip, & !< integration point number - el !< element number + integer, intent(in) :: ce real(pReal), intent(in) :: & phi real(pReal) :: & @@ -121,9 +117,9 @@ module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, p phiDot = 0.0_pReal dPhiDot_dPhi = 0.0_pReal - call phase_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) + call phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ce) + phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) + dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) end subroutine damage_nonlocal_getSourceAndItsTangent @@ -131,20 +127,18 @@ end subroutine damage_nonlocal_getSourceAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief updated nonlocal damage field with solution from damage phase field PDE !-------------------------------------------------------------------------------------------------- -module subroutine damage_nonlocal_putNonLocalDamage(phi,ip,el) +module subroutine damage_nonlocal_putNonLocalDamage(phi,ce) - integer, intent(in) :: & - ip, & !< integration point number - el !< element number + integer, intent(in) :: ce real(pReal), intent(in) :: & phi integer :: & - homog, & - offset + ho, & + me - homog = material_homogenizationAt(el) - offset = material_homogenizationMemberAt(ip,el) - damagestate_h(homog)%state(1,offset) = phi + ho = material_homogenizationAt2(ce) + me = material_homogenizationMemberAt2(ce) + damagestate_h(ho)%state(1,me) = phi end subroutine damage_nonlocal_putNonLocalDamage diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90 index a02f798cd..928bd547b 100644 --- a/src/homogenization_thermal.f90 +++ b/src/homogenization_thermal.f90 @@ -99,20 +99,16 @@ end subroutine thermal_homogenize !-------------------------------------------------------------------------------------------------- !> @brief return homogenized thermal conductivity in reference configuration !-------------------------------------------------------------------------------------------------- -module function thermal_conduction_getConductivity(ip,el) result(K) +module function thermal_conduction_getConductivity(ce) result(K) - integer, intent(in) :: & - ip, & !< integration point number - el !< element number + integer, intent(in) :: ce real(pReal), dimension(3,3) :: K integer :: & - co, & - ce + co K = 0.0_pReal - ce = (el-1)*discretization_nIPs + ip do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) K = K + crystallite_push33ToRef(co,ce,lattice_K(:,:,material_phaseAt2(co,ce))) enddo @@ -220,11 +216,11 @@ end function homogenization_thermal_T !-------------------------------------------------------------------------------------------------- !> @brief return heat generation rate !-------------------------------------------------------------------------------------------------- -module subroutine thermal_conduction_getSource(Tdot, ip,el) +module subroutine thermal_conduction_getSource(Tdot, ip, el) integer, intent(in) :: & - ip, & !< integration point number - el !< element number + ip, & + el real(pReal), intent(out) :: & Tdot diff --git a/src/phase.f90 b/src/phase.f90 index 9ab5c386a..add26ae0b 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -231,10 +231,8 @@ module phase end function phase_homogenizedC - module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) - integer, intent(in) :: & - ip, & !< integration point number - el !< element number + module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ce) + integer, intent(in) :: ce real(pReal), intent(in) :: & phi !< damage parameter real(pReal), intent(inout) :: & diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 5e4e43017..c6fd7bfab 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -179,11 +179,9 @@ end subroutine damage_init !---------------------------------------------------------------------------------------------- !< @brief returns local part of nonlocal damage driving force !---------------------------------------------------------------------------------------------- -module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) +module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ce) - integer, intent(in) :: & - ip, & !< integration point number - el !< element number + integer, intent(in) :: ce real(pReal), intent(in) :: & phi !< damage parameter real(pReal), intent(inout) :: & @@ -201,9 +199,9 @@ module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, phiDot = 0.0_pReal dPhiDot_dPhi = 0.0_pReal - do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - ph = material_phaseAt(co,el) - me = material_phasememberAt(co,ip,el) + do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) + ph = material_phaseAt2(co,ce) + me = material_phasememberAt2(co,ce) select case(phase_source(ph)) case (DAMAGE_ISOBRITTLE_ID)