From ea664688f8011cbaab8e311cf6c32c4d2d98b3a4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 18 Oct 2013 20:56:10 +0000 Subject: [PATCH] introduced dummy temperature calculation. --- code/DAMASK_spectral_driver.f90 | 5 ++- code/DAMASK_spectral_utilities.f90 | 68 ++++++++++++++++++++++++++++-- code/constitutive.f90 | 20 --------- code/homogenization.f90 | 5 ++- code/homogenization_RGC.f90 | 2 +- 5 files changed, 72 insertions(+), 28 deletions(-) diff --git a/code/DAMASK_spectral_driver.f90 b/code/DAMASK_spectral_driver.f90 index e6bfa3196..ce7745dbd 100644 --- a/code/DAMASK_spectral_driver.f90 +++ b/code/DAMASK_spectral_driver.f90 @@ -75,7 +75,8 @@ program DAMASK_spectral_Driver geomSize, & tBoundaryCondition, & tSolutionState, & - cutBack + cutBack, & + utilities_temperatureUpdate use DAMASK_spectral_SolverBasic #ifdef PETSc use DAMASK_spectral_SolverBasicPETSC @@ -522,7 +523,7 @@ program DAMASK_spectral_Driver ' increment ', totalIncsCounter, ' NOT converged' notConvergedCounter = notConvergedCounter + 1_pInt endif; flush(6) - + call utilities_temperatureUpdate() if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency write(6,'(1/,a)') ' ... writing results to file ......................................' write(resUnit) materialpoint_results ! write result to file diff --git a/code/DAMASK_spectral_utilities.f90 b/code/DAMASK_spectral_utilities.f90 index ab9fab231..200ef5b87 100644 --- a/code/DAMASK_spectral_utilities.f90 +++ b/code/DAMASK_spectral_utilities.f90 @@ -109,7 +109,8 @@ module DAMASK_spectral_utilities utilities_constitutiveResponse, & utilities_calculateRate, & utilities_forwardField, & - utilities_destroy + utilities_destroy, & + utilities_temperatureUpdate private :: & utilities_getFilter @@ -861,11 +862,10 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,& endif call CPFEM_general(collectMode,usePingPong,F_lastInc(1:3,1:3,1,1,1),F(1:3,1:3,1,1,1), & ! collect mode handles Jacobian backup / restoration - temperature,timeinc,1_pInt,1_pInt) + crystallite_temperature(1,1),timeinc,1_pInt,1_pInt) materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid)]) materialpoint_F = reshape(F, [3,3,1,product(grid)]) - crystallite_temperature = temperature call debug_reset() @@ -1004,6 +1004,68 @@ real(pReal) function utilities_getFilter(k) end function utilities_getFilter +!-------------------------------------------------------------------------------------------------- +!> @brief calculates filter for fourier convolution depending on type given in numerics.config +!-------------------------------------------------------------------------------------------------- +subroutine utilities_temperatureUpdate() + use crystallite, only: & + crystallite_temperature + use homogenization, only: & + materialpoint_heat + + implicit none + integer x,y,z,e + real(pReal) :: a + do e=1_pInt, product(grid) + if (e==1) print*, materialpoint_heat(1,e)*0.0001_pReal + crystallite_temperature(1,e) = crystallite_temperature(1,e) + materialpoint_heat(1,e)*0.0001_pReal + enddo + e = 0_pInt + z = 0 + y = 0 + x = 0 + +!< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me] + do z = 0_pInt,grid(3)-1_pInt + do y = 0_pInt,grid(2)-1_pInt + do x = 0_pInt,grid(1)-1_pInt + e = e + 1_pInt + a = 0.0_pReal + a = a+ crystallite_temperature(1, z * grid(1) * grid(2) & + + y * grid(1) & + + modulo(x+1_pInt,grid(1)) & + + 1_pInt) + a = a+ crystallite_temperature(1,z * grid(1) * grid(2) & + + y * grid(1) & + + modulo(x-1_pInt,grid(1)) & + + 1_pInt) + a = a+ crystallite_temperature(1,z * grid(1) * grid(2) & + + modulo(y+1_pInt,grid(2)) * grid(1) & + + x & + + 1_pInt) + a = a+ crystallite_temperature(1,z * grid(1) * grid(2) & + + modulo(y-1_pInt,grid(2)) * grid(1) & + + x & + + 1_pInt) + a = a+ crystallite_temperature(1, modulo(z+1_pInt,grid(3)) * grid(1) * grid(2) & + + y * grid(1) & + + x & + + 1_pInt) + a = a+ crystallite_temperature(1,modulo(z-1_pInt,grid(3)) * grid(1) * grid(2) & + + y * grid(1) & + + x & + + 1_pInt) + if (e==1) print*, a + if (e==1) print*, crystallite_temperature(1,e) + if (e==1) print*,(crystallite_temperature(1,e)+a)/7.0_pReal + crystallite_temperature(1,e) = (crystallite_temperature(1,e)+a)/7.0_pReal + enddo + enddo + enddo + +end subroutine utilities_temperatureUpdate + + !-------------------------------------------------------------------------------------------------- !> @brief cleans up !-------------------------------------------------------------------------------------------------- diff --git a/code/constitutive.f90 b/code/constitutive.f90 index c0d72feab..242587a35 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -66,7 +66,6 @@ module constitutive constitutive_TandItsTangent, & constitutive_collectDotState, & constitutive_collectDeltaState, & - constitutive_heat, & constitutive_postResults private :: & @@ -846,25 +845,6 @@ subroutine constitutive_collectDeltaState(Tstar_v, ipc, ip, el) end subroutine constitutive_collectDeltaState -!-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the rate of change of microstructure -!-------------------------------------------------------------------------------------------------- -real(pReal) function constitutive_heat(Tstar_v,Temperature,ipc,ip,el) - implicit none - integer(pInt), intent(in) :: & - ipc, & !< grain number - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - Temperature - real(pReal), intent(in), dimension(6) :: & - Tstar_v !< 2nd Piola-Kirchhoff stress - - constitutive_heat = 0.0_pReal - -end function constitutive_heat - - !-------------------------------------------------------------------------------------------------- !> @brief returns array of constitutive results !-------------------------------------------------------------------------------------------------- diff --git a/code/homogenization.f90 b/code/homogenization.f90 index 3f0ec47e4..4fc99c05a 100644 --- a/code/homogenization.f90 +++ b/code/homogenization.f90 @@ -52,6 +52,8 @@ module homogenization integer(pInt), public, protected :: & materialpoint_sizeResults, & homogenization_maxSizePostResults + real(pReal), dimension(:,:), allocatable, public, protected :: & + materialpoint_heat type(p_vec), dimension(:,:), allocatable, private :: & homogenization_subState0 !< pointer array to homogenization state at start of homogenization increment @@ -61,8 +63,7 @@ module homogenization real(pReal), dimension(:,:), allocatable, private :: & materialpoint_subFrac, & materialpoint_subStep, & - materialpoint_subdt, & - materialpoint_heat + materialpoint_subdt integer(pInt), dimension(:,:), allocatable, private :: & homogenization_sizePostResults !< size of postResults array per material point integer(pInt), private :: & diff --git a/code/homogenization_RGC.f90 b/code/homogenization_RGC.f90 index ed0e64301..587d52123 100644 --- a/code/homogenization_RGC.f90 +++ b/code/homogenization_RGC.f90 @@ -211,7 +211,7 @@ subroutine homogenization_RGC_init(myUnit) do i = 1_pInt,maxNinstance do j = 1_pInt,maxval(homogenization_Noutput) select case(homogenization_RGC_output(j,i)) - case('temperature','constitutivework','penaltyenergy','volumediscrepancy'& + case('temperature','constitutivework','penaltyenergy','volumediscrepancy', & 'averagerelaxrate','maximumrelaxrate') mySize = 1_pInt case('ipcoords','magnitudemismatch')