diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index 0f9d37ddb..9464595b1 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -365,7 +365,7 @@ subroutine flux(f,ts,n,time) f f(2) = 0.0_pReal - call thermal_conduction_getSource(f(1), ts(3), n(3),mesh_FEM2DAMASK_elem(n(1))) + call thermal_conduction_getSource(f(1), n(3),mesh_FEM2DAMASK_elem(n(1))) end subroutine flux diff --git a/src/constitutive.f90 b/src/constitutive.f90 index c6415a883..29796cfbc 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -120,19 +120,15 @@ module constitutive integer, intent(in) :: ph, me end subroutine mech_initializeRestorationPoints - module subroutine thermal_initializeRestorationPoints(ph,me) + module subroutine constitutive_thermal_initializeRestorationPoints(ph,me) integer, intent(in) :: ph, me - end subroutine thermal_initializeRestorationPoints + end subroutine constitutive_thermal_initializeRestorationPoints module subroutine mech_windForward(ph,me) integer, intent(in) :: ph, me end subroutine mech_windForward - module subroutine thermal_windForward(ph,me) - integer, intent(in) :: ph, me - end subroutine thermal_windForward - module subroutine mech_forward() end subroutine mech_forward @@ -146,10 +142,6 @@ module constitutive logical, intent(in) :: includeL end subroutine mech_restore - module subroutine thermal_restore(ip,el) - integer, intent(in) :: ip, el - end subroutine thermal_restore - module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF) real(pReal), intent(in) :: dt @@ -207,21 +199,20 @@ module constitutive integer, intent(in) :: co, ip, el end subroutine constitutive_mech_setF - module subroutine constitutive_thermal_setT(T,co,ip,el) + module subroutine constitutive_thermal_setT(T,co,ce) real(pReal), intent(in) :: T - integer, intent(in) :: co, ip, el + integer, intent(in) :: co, ce end subroutine constitutive_thermal_setT ! == cleaned:end =================================================================================== - module function integrateThermalState(Delta_t,co,ip,el) result(broken) + module function thermal_stress(Delta_t,ph,me) result(converged_) + real(pReal), intent(in) :: Delta_t - integer, intent(in) :: & - el, & !< element index in element loop - ip, & !< integration point index in ip loop - co !< grain index in grain loop - logical :: broken - end function integrateThermalState + integer, intent(in) :: ph, me + logical :: converged_ + + end function thermal_stress module function integrateDamageState(dt,co,ip,el) result(broken) real(pReal), intent(in) :: dt @@ -283,12 +274,10 @@ module constitutive dPhiDot_dPhi end subroutine constitutive_damage_getRateAndItsTangents - module subroutine constitutive_thermal_getRate(TDot, T,ip,el) + module subroutine constitutive_thermal_getRate(TDot, ip,el) integer, intent(in) :: & ip, & !< integration point number el !< element number - real(pReal), intent(in) :: & - T real(pReal), intent(out) :: & TDot end subroutine constitutive_thermal_getRate @@ -394,18 +383,19 @@ module constitutive converged, & crystallite_init, & crystallite_stress, & + thermal_stress, & constitutive_mech_dPdF, & crystallite_orientations, & crystallite_push33ToRef, & constitutive_restartWrite, & constitutive_restartRead, & - integrateThermalState, & integrateDamageState, & constitutive_thermal_setT, & constitutive_mech_getP, & constitutive_mech_setF, & constitutive_mech_getF, & constitutive_initializeRestorationPoints, & + constitutive_thermal_initializeRestorationPoints, & constitutive_windForward, & KINEMATICS_UNDEFINED_ID ,& KINEMATICS_CLEAVAGE_OPENING_ID, & @@ -553,7 +543,6 @@ subroutine constitutive_restore(ip,el,includeL) enddo call mech_restore(ip,el,includeL) - call thermal_restore(ip,el) end subroutine constitutive_restore @@ -720,7 +709,6 @@ subroutine constitutive_initializeRestorationPoints(ip,el) me = material_phaseMemberAt(co,ip,el) call mech_initializeRestorationPoints(ph,me) - call thermal_initializeRestorationPoints(ph,me) do so = 1, size(damageState(ph)%p) damageState(ph)%p(so)%partitionedState0(:,me) = damageState(ph)%p(so)%state0(:,me) @@ -750,7 +738,6 @@ subroutine constitutive_windForward(ip,el) me = material_phaseMemberAt(co,ip,el) call mech_windForward(ph,me) - call thermal_windForward(ph,me) do so = 1, phase_Nsources(material_phaseAt(co,el)) damageState(ph)%p(so)%partitionedState0(:,me) = damageState(ph)%p(so)%state(:,me) diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index 8bc85354f..50afe2e87 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -599,20 +599,18 @@ module subroutine constitutive_plastic_dependentState(co, ip, el) el !< element integer :: & - ho, & !< homogenization - tme, & !< thermal member position + ph, & instance, me - ho = material_homogenizationAt(el) - tme = material_homogenizationMemberAt(ip,el) - me = material_phasememberAt(co,ip,el) - instance = phase_plasticityInstance(material_phaseAt(co,el)) + ph = material_phaseAt(co,el) + me = material_phasememberAt(co,ip,el) + instance = phase_plasticityInstance(ph) plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) case (PLASTICITY_DISLOTWIN_ID) plasticityType - call plastic_dislotwin_dependentState(temperature(ho)%p(tme),instance,me) + call plastic_dislotwin_dependentState(thermal_T(ph,me),instance,me) case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType call plastic_dislotungsten_dependentState(instance,me) @@ -650,17 +648,13 @@ subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & real(pReal), dimension(3,3) :: & Mp !< Mandel stress work conjugate with Lp integer :: & - ho, & !< homogenization - tme !< thermal member position - integer :: & - i, j, instance, me + i, j, instance, me, ph - ho = material_homogenizationAt(el) - tme = material_homogenizationMemberAt(ip,el) Mp = matmul(matmul(transpose(Fi),Fi),S) me = material_phasememberAt(co,ip,el) - instance = phase_plasticityInstance(material_phaseAt(co,el)) + ph = material_phaseAt(co,el) + instance = phase_plasticityInstance(ph) plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) @@ -678,13 +672,13 @@ subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & call plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me) case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, temperature(ho)%p(tme),instance,me,ip,el) + call plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),instance,me,ip,el) case (PLASTICITY_DISLOTWIN_ID) plasticityType - call plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,me) + call plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),instance,me) case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType - call plastic_dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,me) + call plastic_dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),instance,me) end select plasticityType @@ -700,52 +694,49 @@ end subroutine constitutive_plastic_LpAndItsTangents !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -function mech_collectDotState(subdt,co,ip,el,ph,of) result(broken) +function mech_collectDotState(subdt,co,ip,el,ph,me) result(broken) integer, intent(in) :: & - co, & !< component-ID of integration point + co, & !< component-ID of integration point ip, & !< integration point el, & !< element ph, & - of + me real(pReal), intent(in) :: & subdt !< timestep real(pReal), dimension(3,3) :: & Mp integer :: & - ho, & !< homogenization - tme, & !< thermal member position instance logical :: broken - ho = material_homogenizationAt(el) - tme = material_homogenizationMemberAt(ip,el) + instance = phase_plasticityInstance(ph) - Mp = matmul(matmul(transpose(constitutive_mech_Fi(ph)%data(1:3,1:3,of)),& - constitutive_mech_Fi(ph)%data(1:3,1:3,of)),constitutive_mech_S(ph)%data(1:3,1:3,of)) + Mp = matmul(matmul(transpose(constitutive_mech_Fi(ph)%data(1:3,1:3,me)),& + constitutive_mech_Fi(ph)%data(1:3,1:3,me)),constitutive_mech_S(ph)%data(1:3,1:3,me)) plasticityType: select case (phase_plasticity(ph)) case (PLASTICITY_ISOTROPIC_ID) plasticityType - call plastic_isotropic_dotState(Mp,instance,of) + call plastic_isotropic_dotState(Mp,instance,me) case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType - call plastic_phenopowerlaw_dotState(Mp,instance,of) + call plastic_phenopowerlaw_dotState(Mp,instance,me) case (PLASTICITY_KINEHARDENING_ID) plasticityType - call plastic_kinehardening_dotState(Mp,instance,of) + call plastic_kinehardening_dotState(Mp,instance,me) case (PLASTICITY_DISLOTWIN_ID) plasticityType - call plastic_dislotwin_dotState(Mp,temperature(ho)%p(tme),instance,of) + call plastic_dislotwin_dotState(Mp,thermal_T(ph,me),instance,me) case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType - call plastic_disloTungsten_dotState(Mp,temperature(ho)%p(tme),instance,of) + call plastic_disloTungsten_dotState(Mp,thermal_T(ph,me),instance,me) case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_dotState(Mp,temperature(ho)%p(tme),subdt,instance,of,ip,el) + call plastic_nonlocal_dotState(Mp,thermal_T(ph,me),subdt,instance,me,ip,el) end select plasticityType - broken = any(IEEE_is_NaN(plasticState(ph)%dotState(:,of))) + broken = any(IEEE_is_NaN(plasticState(ph)%dotState(:,me))) end function mech_collectDotState @@ -1633,9 +1624,6 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) do so = 1, phase_Nsources(ph) damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state(:,me) enddo - do so = 1, thermal_Nsources(ph) - thermalState(ph)%p(so)%subState0(:,me) = thermalState(ph)%p(so)%state(:,me) - enddo endif !-------------------------------------------------------------------------------------------------- ! cut back (reduced time and restore) @@ -1652,9 +1640,6 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) do so = 1, phase_Nsources(ph) damageState(ph)%p(so)%state(:,me) = damageState(ph)%p(so)%subState0(:,me) enddo - do so = 1, thermal_Nsources(ph) - thermalState(ph)%p(so)%state(:,me) = thermalState(ph)%p(so)%subState0(:,me) - enddo todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair) endif @@ -1668,7 +1653,6 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) constitutive_mech_Fp(ph)%data(1:3,1:3,me)))) converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * dt,co,ip,el) converged_ = converged_ .and. .not. integrateDamageState(subStep * dt,co,ip,el) - converged_ = converged_ .and. .not. integrateThermalState(subStep * dt,co,ip,el) endif enddo cutbackLooping diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index 9787cb0e4..9007ef729 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -86,7 +86,7 @@ module subroutine thermal_init(phases) Nconstituents = count(material_phaseAt == ph) * discretization_nIPs - allocate(current(ph)%T(Nconstituents)) + allocate(current(ph)%T(Nconstituents),source=300.0_pReal) phase => phases%get(ph) if(phase%contains('thermal')) then thermal => phase%get('thermal') @@ -127,13 +127,11 @@ end subroutine thermal_init !---------------------------------------------------------------------------------------------- !< @brief calculates thermal dissipation rate !---------------------------------------------------------------------------------------------- -module subroutine constitutive_thermal_getRate(TDot, T, ip, el) +module subroutine constitutive_thermal_getRate(TDot, ip, el) integer, intent(in) :: & ip, & !< integration point number el !< element number - real(pReal), intent(in) :: & - T !< plastic velocity gradient real(pReal), intent(out) :: & TDot @@ -197,30 +195,37 @@ function constitutive_thermal_collectDotState(ph,me) result(broken) end function constitutive_thermal_collectDotState +module function thermal_stress(Delta_t,ph,me) result(converged_) + + real(pReal), intent(in) :: Delta_t + integer, intent(in) :: ph, me + logical :: converged_ + + integer :: so + + + do so = 1, thermal_Nsources(ph) + thermalState(ph)%p(so)%state(:,me) = thermalState(ph)%p(so)%subState0(:,me) + enddo + converged_ = .not. integrateThermalState(Delta_t,ph,me) + +end function thermal_stress + !-------------------------------------------------------------------------------------------------- !> @brief integrate state with 1st order explicit Euler method !-------------------------------------------------------------------------------------------------- -module function integrateThermalState(Delta_t,co,ip,el) result(broken) +function integrateThermalState(Delta_t, ph,me) result(broken) real(pReal), intent(in) :: Delta_t - integer, intent(in) :: & - el, & !< element index in element loop - ip, & !< integration point index in ip loop - co !< grain index in grain loop + integer, intent(in) :: ph, me logical :: & broken integer :: & - ph, & - me, & so, & sizeDotState - - ph = material_phaseAt(co,el) - me = material_phaseMemberAt(co,ip,el) - broken = constitutive_thermal_collectDotState(ph,me) if(broken) return @@ -233,7 +238,7 @@ module function integrateThermalState(Delta_t,co,ip,el) result(broken) end function integrateThermalState -module subroutine thermal_initializeRestorationPoints(ph,me) +module subroutine constitutive_thermal_initializeRestorationPoints(ph,me) integer, intent(in) :: ph, me @@ -244,24 +249,10 @@ module subroutine thermal_initializeRestorationPoints(ph,me) thermalState(ph)%p(so)%partitionedState0(:,me) = thermalState(ph)%p(so)%state0(:,me) enddo -end subroutine thermal_initializeRestorationPoints +end subroutine constitutive_thermal_initializeRestorationPoints -module subroutine thermal_windForward(ph,me) - - integer, intent(in) :: ph, me - - integer :: so - - - do so = 1, size(thermalState(ph)%p) - thermalState(ph)%p(so)%partitionedState0(:,me) = thermalState(ph)%p(so)%state(:,me) - enddo - -end subroutine thermal_windForward - - module subroutine thermal_forward() integer :: ph, so @@ -276,26 +267,6 @@ module subroutine thermal_forward() end subroutine thermal_forward -module subroutine thermal_restore(ip,el) - - integer, intent(in) :: ip, el - - integer :: co, ph, me, so - - - do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - ph = material_phaseAt(co,el) - me = material_phaseMemberAt(co,ip,el) - - do so = 1, size(thermalState(ph)%p) - thermalState(ph)%p(so)%state(:,me) = thermalState(ph)%p(so)%partitionedState0(:,me) - enddo - - enddo - -end subroutine thermal_restore - - !---------------------------------------------------------------------------------------------- !< @brief Get temperature (for use by non-thermal physics) !---------------------------------------------------------------------------------------------- @@ -313,13 +284,13 @@ end function thermal_T !---------------------------------------------------------------------------------------------- !< @brief Set temperature !---------------------------------------------------------------------------------------------- -module subroutine constitutive_thermal_setT(T,co,ip,el) +module subroutine constitutive_thermal_setT(T,co,ce) real(pReal), intent(in) :: T - integer, intent(in) :: co, ip, el + integer, intent(in) :: ce, co - current(material_phaseAt(co,el))%T(material_phaseMemberAt(co,ip,el)) = T + current(material_phaseAt2(co,ce))%T(material_phaseMemberAt2(co,ce)) = T end subroutine constitutive_thermal_setT diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 9d804ec56..d164f11e7 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -16,10 +16,11 @@ module grid_thermal_spectral use spectral_utilities use discretization_grid use thermal_conduction + use homogenization use YAML_types use config use material - + implicit none private @@ -45,11 +46,11 @@ module grid_thermal_spectral T_stagInc !< field of staggered temperature !-------------------------------------------------------------------------------------------------- -! reference diffusion tensor, mobility etc. +! reference diffusion tensor, mobility etc. integer :: totalIter = 0 !< total iteration in current increment real(pReal), dimension(3,3) :: K_ref real(pReal) :: mu_ref - + public :: & grid_thermal_spectral_init, & grid_thermal_spectral_solution, & @@ -63,8 +64,8 @@ contains !-------------------------------------------------------------------------------------------------- subroutine grid_thermal_spectral_init - PetscInt, dimension(0:worldsize-1) :: localK - integer :: i, j, k, cell + PetscInt, dimension(0:worldsize-1) :: localK + integer :: i, j, k, ce DM :: thermal_grid PetscScalar, dimension(:,:,:), pointer :: x_scal PetscErrorCode :: ierr @@ -93,11 +94,11 @@ subroutine grid_thermal_spectral_init CHKERRQ(ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),ierr) CHKERRQ(ierr) - + !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc call SNESCreate(PETSC_COMM_WORLD,thermal_snes,ierr); CHKERRQ(ierr) - call SNESSetOptionsPrefix(thermal_snes,'thermal_',ierr);CHKERRQ(ierr) + call SNESSetOptionsPrefix(thermal_snes,'thermal_',ierr);CHKERRQ(ierr) localK = 0 localK(worldrank) = grid3 call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) @@ -115,23 +116,23 @@ subroutine grid_thermal_spectral_init call DMsetUp(thermal_grid,ierr); CHKERRQ(ierr) call DMCreateGlobalVector(thermal_grid,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 1, i.e. every def grad tensor) call DMDASNESSetFunctionLocal(thermal_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector - CHKERRQ(ierr) + CHKERRQ(ierr) call SNESSetFromOptions(thermal_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments !-------------------------------------------------------------------------------------------------- -! init fields +! init fields call DMDAGetCorners(thermal_grid,xstart,ystart,zstart,xend,yend,zend,ierr) CHKERRQ(ierr) xend = xstart + xend - 1 yend = ystart + yend - 1 - zend = zstart + zend - 1 + zend = zstart + zend - 1 allocate(T_current(grid(1),grid(2),grid3), source=0.0_pReal) allocate(T_lastInc(grid(1),grid(2),grid3), source=0.0_pReal) allocate(T_stagInc(grid(1),grid(2),grid3), source=0.0_pReal) - cell = 0 + ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 - T_current(i,j,k) = temperature(material_homogenizationAt(cell))%p(material_homogenizationMemberAt(1,cell)) + ce = ce + 1 + T_current(i,j,k) = temperature(material_homogenizationAt(ce))%p(material_homogenizationMemberAt(1,ce)) T_lastInc(i,j,k) = T_current(i,j,k) T_stagInc(i,j,k) = T_current(i,j,k) enddo; enddo; enddo @@ -143,26 +144,26 @@ subroutine grid_thermal_spectral_init end subroutine grid_thermal_spectral_init - + !-------------------------------------------------------------------------------------------------- !> @brief solution for the spectral thermal scheme with internal iterations !-------------------------------------------------------------------------------------------------- function grid_thermal_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 :: T_min, T_max, stagNorm, solnNorm - PetscErrorCode :: ierr + PetscErrorCode :: ierr SNESConvergedReason :: reason solution%converged =.false. - + !-------------------------------------------------------------------------------------------------- -! set module wide availabe data +! set module wide availabe data params%timeinc = timeinc call SNESSolve(thermal_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) @@ -183,13 +184,14 @@ function grid_thermal_spectral_solution(timeinc) result(solution) solution%stagConverged = stagNorm < max(num%eps_thermal_atol, num%eps_thermal_rtol*solnNorm) !-------------------------------------------------------------------------------------------------- -! updating thermal state - cell = 0 +! updating thermal state + ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 + ce = ce + 1 call thermal_conduction_putTemperatureAndItsRate(T_current(i,j,k), & (T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc, & - 1,cell) + 1,ce) + homogenization_T(ce) = T_current(i,j,k) enddo; enddo; enddo call VecMin(solution_vec,devNull,T_min,ierr); CHKERRQ(ierr) @@ -198,7 +200,7 @@ function grid_thermal_spectral_solution(timeinc) result(solution) print'(/,a)', ' ... thermal conduction converged ..................................' print'(/,a,f8.4,2x,f8.4,2x,f8.4)', ' Minimum|Maximum|Delta Temperature / K = ', T_min, T_max, stagNorm print'(/,a)', ' ===========================================================================' - flush(IO_STDOUT) + flush(IO_STDOUT) end function grid_thermal_spectral_solution @@ -207,36 +209,37 @@ end function grid_thermal_spectral_solution !> @brief forwarding routine !-------------------------------------------------------------------------------------------------- subroutine grid_thermal_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 - - if (cutBack) then + + if (cutBack) then T_current = T_lastInc T_stagInc = T_lastInc !-------------------------------------------------------------------------------------------------- -! reverting thermal field state - cell = 0 +! reverting thermal field state + ce = 0 call SNESGetDM(thermal_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) = T_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 + ce = ce + 1 call thermal_conduction_putTemperatureAndItsRate(T_current(i,j,k), & (T_current(i,j,k) - & T_lastInc(i,j,k))/params%timeinc, & - 1,cell) + 1,ce) + homogenization_T(ce) = T_current(i,j,k) enddo; enddo; enddo else T_lastInc = T_current call updateReference endif - + end subroutine grid_thermal_spectral_forward @@ -244,7 +247,7 @@ end subroutine grid_thermal_spectral_forward !> @brief forms the spectral thermal residual vector !-------------------------------------------------------------------------------------------------- subroutine formResidual(in,x_scal,f_scal,dummy,ierr) - + DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & in PetscScalar, dimension( & @@ -255,33 +258,33 @@ 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) :: Tdot - T_current = x_scal + T_current = x_scal !-------------------------------------------------------------------------------------------------- ! evaluate polarization field scalarField_real = 0.0_pReal - scalarField_real(1:grid(1),1:grid(2),1:grid3) = T_current + scalarField_real(1:grid(1),1:grid(2),1:grid3) = T_current call utilities_FFTscalarForward call utilities_fourierScalarGradient !< calculate gradient of temperature 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(thermal_conduction_getConductivity(1,cell) - K_ref, & + 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)) enddo; enddo; enddo call utilities_FFTvectorForward call utilities_fourierVectorDivergence !< calculate temperature 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 thermal_conduction_getSource(Tdot, T_current(i,j,k), 1, cell) + 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 (1,cell)* & - thermal_conduction_getSpecificHeat(1,cell)*(T_lastInc(i,j,k) - & + + thermal_conduction_getMassDensity (1,ce)* & + thermal_conduction_getSpecificHeat(1,ce)*(T_lastInc(i,j,k) - & T_current(i,j,k))& + mu_ref*T_current(i,j,k) enddo; enddo; enddo @@ -291,7 +294,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) call utilities_FFTscalarForward call utilities_fourierGreenConvolution(K_ref, mu_ref, params%timeinc) call utilities_FFTscalarBackward - + !-------------------------------------------------------------------------------------------------- ! constructing residual f_scal = T_current - scalarField_real(1:grid(1),1:grid(2),1:grid3) @@ -304,15 +307,15 @@ end subroutine formResidual !-------------------------------------------------------------------------------------------------- subroutine updateReference - integer :: i,j,k,cell,ierr - - cell = 0 + integer :: i,j,k,ce,ierr + + 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 + thermal_conduction_getConductivity(1,cell) - mu_ref = mu_ref + thermal_conduction_getMassDensity(1,cell)* thermal_conduction_getSpecificHeat(1,cell) + ce = ce + 1 + K_ref = K_ref + thermal_conduction_getConductivity(1,ce) + mu_ref = mu_ref + thermal_conduction_getMassDensity(1,ce)* thermal_conduction_getSpecificHeat(1,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 9112562b9..a63487b6f 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -28,7 +28,8 @@ module homogenization !-------------------------------------------------------------------------------------------------- ! General variables for the homogenization at a material point real(pReal), dimension(:), allocatable, public :: & - homogenization_T + homogenization_T, & + homogenization_dot_T real(pReal), dimension(:,:,:), allocatable, public :: & homogenization_F0, & !< def grad of IP at start of FE increment homogenization_F !< def grad of IP to be reached at end of FE increment @@ -69,13 +70,15 @@ module homogenization el !< element number end subroutine mech_partition - module subroutine thermal_partition(T,ip,el) + module subroutine thermal_partition(T,ce) real(pReal), intent(in) :: T - integer, intent(in) :: & - ip, & !< integration point - el !< element number + integer, intent(in) :: ce end subroutine thermal_partition + module subroutine thermal_homogenize(ip,el) + integer, intent(in) :: ip,el + end subroutine thermal_homogenize + module subroutine mech_homogenize(dt,ip,el) real(pReal), intent(in) :: dt integer, intent(in) :: & @@ -161,7 +164,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE NiterationMPstate, & ip, & !< integration point number el, & !< element number - myNgrains, co, ce, ho, me + myNgrains, co, ce, ho, me, ph real(pReal) :: & subFrac, & subStep @@ -170,8 +173,8 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE logical, dimension(2) :: & doneAndHappy - - !$OMP PARALLEL DO PRIVATE(ce,me,ho,myNgrains,NiterationMPstate,subFrac,converged,subStep,doneAndHappy) + !$OMP PARALLEL + !$OMP DO PRIVATE(ce,me,ho,myNgrains,NiterationMPstate,subFrac,converged,subStep,doneAndHappy) do el = FEsolving_execElem(1),FEsolving_execElem(2) ho = material_homogenizationAt(el) myNgrains = homogenization_Nconstituents(ho) @@ -221,9 +224,8 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE if (subStep > num%subStepMinHomog) doneAndHappy = [.false.,.true.] NiterationMPstate = 0 - convergenceLooping: do while (.not. terminallyIll & - .and. .not. doneAndHappy(1) & - .and. NiterationMPstate < num%nMPstate) + convergenceLooping: do while (.not. (terminallyIll .or. doneAndHappy(1)) & + .and. NiterationMPstate < num%nMPstate) NiterationMPstate = NiterationMPstate + 1 !-------------------------------------------------------------------------------------------------- @@ -231,10 +233,9 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE if (.not. doneAndHappy(1)) then ce = (el-1)*discretization_nIPs + ip - call mech_partition(homogenization_F0(1:3,1:3,ce) & - + (homogenization_F(1:3,1:3,ce)-homogenization_F0(1:3,1:3,ce))& - *(subStep+subFrac), & - ip,el) + call mech_partition( homogenization_F0(1:3,1:3,ce) & + + (homogenization_F(1:3,1:3,ce)-homogenization_F0(1:3,1:3,ce))*(subStep+subFrac), & + ip,el) converged = .true. do co = 1, myNgrains converged = converged .and. crystallite_stress(dt*subStep,co,ip,el) @@ -257,24 +258,45 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE enddo cutBackLooping enddo enddo - !$OMP END PARALLEL DO + !$OMP END DO if (.not. terminallyIll ) then - !$OMP PARALLEL DO PRIVATE(ho,myNgrains) + !$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 thermal_partition(homogenization_T(ce),ce) + do co = 1, homogenization_Nconstituents(ho) + ph = material_phaseAt(co,el) + call constitutive_thermal_initializeRestorationPoints(ph,material_phaseMemberAt(co,ip,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) - myNgrains = homogenization_Nconstituents(ho) IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2) - do co = 1, myNgrains + do co = 1, homogenization_Nconstituents(ho) call crystallite_orientations(co,ip,el) enddo call mech_homogenize(dt,ip,el) enddo IpLooping3 enddo elementLooping3 - !$OMP END PARALLEL DO + !$OMP END DO else print'(/,a,/)', ' << HOMOG >> Material Point terminally ill' endif + !$OMP END PARALLEL end subroutine materialpoint_stressAndItsTangent diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90 index f181d97fa..de91a98d6 100644 --- a/src/homogenization_thermal.f90 +++ b/src/homogenization_thermal.f90 @@ -15,25 +15,38 @@ module subroutine thermal_init() print'(/,a)', ' <<<+- homogenization_thermal init -+>>>' allocate(homogenization_T(discretization_nIPs*discretization_Nelems)) - + allocate(homogenization_dot_T(discretization_nIPs*discretization_Nelems)) end subroutine thermal_init !-------------------------------------------------------------------------------------------------- -!> @brief Partition T onto the individual constituents. +!> @brief Partition temperature onto the individual constituents. !-------------------------------------------------------------------------------------------------- -module subroutine thermal_partition(T,ip,el) +module subroutine thermal_partition(T,ce) real(pReal), intent(in) :: T - integer, intent(in) :: & - ip, & !< integration point - el !< element number + integer, intent(in) :: ce + integer :: co - call constitutive_thermal_setT(T,1,ip,el) + do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) + call constitutive_thermal_setT(T,co,ce) + enddo end subroutine thermal_partition +!-------------------------------------------------------------------------------------------------- +!> @brief Homogenize temperature rates +!-------------------------------------------------------------------------------------------------- +module subroutine thermal_homogenize(ip,el) + + integer, intent(in) :: ip,el + + call constitutive_thermal_getRate(homogenization_dot_T((el-1)*discretization_nIPs+ip), ip,el) + +end subroutine thermal_homogenize + + end submodule homogenization_thermal diff --git a/src/material.f90 b/src/material.f90 index 16116ca91..e165bd466 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -51,11 +51,15 @@ module material thermal_initialT !< initial temperature per each homogenization integer, dimension(:), allocatable, public, protected :: & ! (elem) - material_homogenizationAt !< homogenization ID of each element + material_homogenizationAt, & !< homogenization ID of each element + material_homogenizationAt2, & !< per cell + material_homogenizationMemberAt2 !< cell integer, dimension(:,:), allocatable, public, protected :: & ! (ip,elem) material_homogenizationMemberAt !< position of the element within its homogenization instance integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem) - material_phaseAt !< phase ID of each element + material_phaseAt, & !< phase ID of each element + material_phaseAt2, & !< per constituent,cell + material_phaseMemberAt2 !< per constituent, cell integer, dimension(:,:,:), allocatable, public, protected :: & ! (constituent,IP,elem) material_phaseMemberAt !< position of the element within its phase instance @@ -215,8 +219,8 @@ subroutine material_parseMaterial real(pReal) :: & frac integer :: & - e, i, c, & - h + el, ip, co, & + h, ce materials => config_material%get('material') phases => config_material%get('phase') @@ -241,29 +245,41 @@ subroutine material_parseMaterial allocate(material_phaseAt(homogenization_maxNconstituents,discretization_Nelems),source=0) allocate(material_phaseMemberAt(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems),source=0) + + allocate(material_homogenizationAt2(discretization_nIPs*discretization_Nelems),source=0) + allocate(material_homogenizationMemberAt2(discretization_nIPs*discretization_Nelems),source=0) + allocate(material_phaseAt2(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0) + allocate(material_phaseMemberAt2(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0) + allocate(material_orientation0(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems)) - do e = 1, discretization_Nelems - material => materials%get(discretization_materialAt(e)) + do el = 1, discretization_Nelems + material => materials%get(discretization_materialAt(el)) constituents => material%get('constituents') - material_homogenizationAt(e) = homogenizations%getIndex(material%get_asString('homogenization')) - do i = 1, discretization_nIPs - counterHomogenization(material_homogenizationAt(e)) = counterHomogenization(material_homogenizationAt(e)) + 1 - material_homogenizationMemberAt(i,e) = counterHomogenization(material_homogenizationAt(e)) + material_homogenizationAt(el) = homogenizations%getIndex(material%get_asString('homogenization')) + do ip = 1, discretization_nIPs + ce = (el-1)*discretization_nIPs + ip + counterHomogenization(material_homogenizationAt(el)) = counterHomogenization(material_homogenizationAt(el)) + 1 + material_homogenizationMemberAt(ip,el) = counterHomogenization(material_homogenizationAt(el)) + material_homogenizationAt2(ce) = material_homogenizationAt(el) + material_homogenizationMemberAt2(ce) = material_homogenizationMemberAt(ip,el) enddo frac = 0.0_pReal - do c = 1, constituents%length - constituent => constituents%get(c) + do co = 1, constituents%length + constituent => constituents%get(co) frac = frac + constituent%get_asFloat('fraction') - material_phaseAt(c,e) = phases%getIndex(constituent%get_asString('phase')) - do i = 1, discretization_nIPs - counterPhase(material_phaseAt(c,e)) = counterPhase(material_phaseAt(c,e)) + 1 - material_phaseMemberAt(c,i,e) = counterPhase(material_phaseAt(c,e)) + material_phaseAt(co,el) = phases%getIndex(constituent%get_asString('phase')) + do ip = 1, discretization_nIPs + ce = (el-1)*discretization_nIPs + ip + counterPhase(material_phaseAt(co,el)) = counterPhase(material_phaseAt(co,el)) + 1 + material_phaseMemberAt(co,ip,el) = counterPhase(material_phaseAt(co,el)) - call material_orientation0(c,i,e)%fromQuaternion(constituent%get_asFloats('O',requiredSize=4)) ! should be done in crystallite + material_phaseAt2(co,ce) = material_phaseAt(co,el) + material_phaseMemberAt2(co,ce) = material_phaseMemberAt(co,ip,el) + call material_orientation0(co,ip,el)%fromQuaternion(constituent%get_asFloats('O',requiredSize=4)) ! should be done in crystallite enddo enddo diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index 79fe0d6cd..1711c5512 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -91,13 +91,11 @@ end subroutine thermal_conduction_init !-------------------------------------------------------------------------------------------------- !> @brief return heat generation rate !-------------------------------------------------------------------------------------------------- -subroutine thermal_conduction_getSource(Tdot, T,ip,el) +subroutine thermal_conduction_getSource(Tdot, ip,el) integer, intent(in) :: & ip, & !< integration point number el !< element number - real(pReal), intent(in) :: & - T real(pReal), intent(out) :: & Tdot @@ -105,7 +103,7 @@ subroutine thermal_conduction_getSource(Tdot, T,ip,el) homog homog = material_homogenizationAt(el) - call constitutive_thermal_getRate(TDot, T,ip,el) + call constitutive_thermal_getRate(TDot, ip,el) Tdot = Tdot/real(homogenization_Nconstituents(homog),pReal)