From a14070bad4729775f1a04479218cbadc50733f65 Mon Sep 17 00:00:00 2001 From: Pratheek Shanthraj Date: Mon, 1 Jun 2015 16:02:27 +0000 Subject: [PATCH] changed up handling of delta states, and some bug fixes --- code/DAMASK_spectral_solverAL.f90 | 70 ++++++++------------- code/DAMASK_spectral_solverPolarisation.f90 | 69 ++++++++------------ code/constitutive.f90 | 11 +--- code/crystallite.f90 | 45 ++++++++----- code/damage_nonlocal.f90 | 3 + code/plastic_disloKMC.f90 | 6 +- code/plastic_disloUCLA.f90 | 6 +- code/plastic_dislotwin.f90 | 6 +- code/plastic_j2.f90 | 6 +- code/plastic_none.f90 | 6 +- code/plastic_nonlocal.f90 | 6 +- code/plastic_phenopowerlaw.f90 | 5 +- code/plastic_titanmod.f90 | 6 +- code/prec.f90 | 1 + code/source_damage_anisoBrittle.f90 | 10 +-- code/source_damage_anisoDuctile.f90 | 10 +-- code/source_damage_isoBrittle.f90 | 10 +-- code/source_damage_isoDuctile.f90 | 12 ++-- code/source_thermal_dissipation.f90 | 6 +- code/source_vacancy_irradiation.f90 | 6 +- code/source_vacancy_phenoplasticity.f90 | 6 +- code/source_vacancy_thermalfluc.f90 | 6 +- 22 files changed, 162 insertions(+), 150 deletions(-) diff --git a/code/DAMASK_spectral_solverAL.f90 b/code/DAMASK_spectral_solverAL.f90 index f8836e401..4c5f9dd6d 100644 --- a/code/DAMASK_spectral_solverAL.f90 +++ b/code/DAMASK_spectral_solverAL.f90 @@ -13,7 +13,8 @@ module DAMASK_spectral_solverAL use math, only: & math_I3 use DAMASK_spectral_utilities, only: & - tSolutionState + tSolutionState, & + tSolutionParams implicit none private @@ -24,14 +25,6 @@ module DAMASK_spectral_solverAL !-------------------------------------------------------------------------------------------------- ! derived types - type tSolutionParams !< @todo use here the type definition for a full loadcase including mask - real(pReal), dimension(3,3) :: P_BC, rotation_BC - real(pReal) :: timeinc - real(pReal) :: timeincOld - real(pReal) :: temperature - real(pReal) :: density - end type tSolutionParams - type(tSolutionParams), private :: params real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal @@ -96,9 +89,7 @@ module DAMASK_spectral_solverAL SNESSetConvergenceTest, & SNESSetFromOptions, & SNESCreate, & - MPI_Abort, & - MPI_Bcast, & - MPI_Allreduce + MPI_Abort contains @@ -124,12 +115,9 @@ subroutine AL_init(temperature) use DAMASK_interface, only: & getSolverJobName use DAMASK_spectral_Utilities, only: & - Utilities_init, & Utilities_constitutiveResponse, & Utilities_updateGamma, & - Utilities_updateIPcoords, & - grid1Red, & - wgt + Utilities_updateIPcoords use mesh, only: & gridLocal, & gridGlobal @@ -150,7 +138,6 @@ subroutine AL_init(temperature) integer(pInt) :: proc character(len=1024) :: rankStr - call Utilities_init() if (worldrank == 0_pInt) then write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverAL init -+>>>' write(6,'(a)') ' $Id$' @@ -169,6 +156,7 @@ subroutine AL_init(temperature) !-------------------------------------------------------------------------------------------------- ! PETSc Init call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) + call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) allocate(localK(worldsize), source = 0); localK(worldrank+1) = gridLocal(3) do proc = 1, worldsize call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr) @@ -182,6 +170,7 @@ subroutine AL_init(temperature) gridLocal (1),gridLocal (2),localK, & ! local grid da,ierr) ! handle, error CHKERRQ(ierr) + call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) call DMDASNESSetFunctionLocal(da,INSERT_VALUES,AL_formResidual,dummy,ierr) CHKERRQ(ierr) @@ -266,7 +255,7 @@ end subroutine AL_init !-------------------------------------------------------------------------------------------------- type(tSolutionState) function & AL_solution(incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,temperature_bc, & - rotation_BC,density) + rotation_BC) use numerics, only: & update_gamma use math, only: & @@ -287,8 +276,7 @@ type(tSolutionState) function & timeinc, & !< increment in time for current solution timeinc_old, & !< increment in time of last increment loadCaseTime, & !< remaining time of current load case - temperature_bc, & - density + temperature_bc logical, intent(in) :: & guess type(tBoundaryCondition), intent(in) :: & @@ -324,7 +312,6 @@ type(tSolutionState) function & params%timeinc = timeinc params%timeincOld = timeinc_old params%temperature = temperature_bc - params%density = density !-------------------------------------------------------------------------------------------------- ! solve BVP @@ -363,16 +350,14 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) math_transpose33, & math_mul3333xx33, & math_invSym3333, & - math_mul33x33, & - PI + math_mul33x33 use DAMASK_spectral_Utilities, only: & wgt, & - field_realMPI, & - field_fourierMPI, & - Utilities_FFTforward, & - Utilities_fourierConvolution, & + tensorField_realMPI, & + utilities_FFTtensorForward, & + utilities_fourierGammaConvolution, & Utilities_inverseLaplace, & - Utilities_FFTbackward, & + utilities_FFTtensorBackward, & Utilities_constitutiveResponse, & Utilities_divergenceRMS, & Utilities_curlRMS @@ -444,9 +429,9 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! - field_realMPI = 0.0_pReal + tensorField_realMPI = 0.0_pReal do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt, gridLocal(1) - field_realMPI(1:3,1:3,i,j,k) = & + tensorField_realMPI(1:3,1:3,i,j,k) = & polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& polarAlpha*math_mul33x33(F(1:3,1:3,i,j,k), & math_mul3333xx33(C_scale,F_lambda(1:3,1:3,i,j,k) - math_I3)) @@ -455,13 +440,13 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! doing convolution in Fourier space - call Utilities_FFTforward() - call Utilities_fourierConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) - call Utilities_FFTbackward() + call utilities_FFTtensorForward() + call utilities_fourierGammaConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) + call utilities_FFTtensorBackward() !-------------------------------------------------------------------------------------------------- ! constructing residual - residual_F_lambda = polarBeta*F - field_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) + residual_F_lambda = polarBeta*F - tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response @@ -473,11 +458,11 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! calculate divergence - field_realMPI = 0.0_pReal - field_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = residual_F - call Utilities_FFTforward() + tensorField_realMPI = 0.0_pReal + tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = residual_F + call utilities_FFTtensorForward() err_div = Utilities_divergenceRMS() - call Utilities_FFTbackward() + call utilities_FFTtensorBackward() !-------------------------------------------------------------------------------------------------- ! constructing residual @@ -493,11 +478,11 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! calculating curl - field_realMPI = 0.0_pReal - field_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = F - call Utilities_FFTforward() + tensorField_realMPI = 0.0_pReal + tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = F + call utilities_FFTtensorForward() err_curl = Utilities_curlRMS() - call Utilities_FFTbackward() + call utilities_FFTtensorBackward() end subroutine AL_formResidual @@ -727,7 +712,6 @@ subroutine AL_destroy() call VecDestroy(solution_vec,ierr); CHKERRQ(ierr) call SNESDestroy(snes,ierr); CHKERRQ(ierr) call DMDestroy(da,ierr); CHKERRQ(ierr) - call Utilities_destroy() end subroutine AL_destroy diff --git a/code/DAMASK_spectral_solverPolarisation.f90 b/code/DAMASK_spectral_solverPolarisation.f90 index fc65e7053..874dea863 100644 --- a/code/DAMASK_spectral_solverPolarisation.f90 +++ b/code/DAMASK_spectral_solverPolarisation.f90 @@ -13,7 +13,8 @@ module DAMASK_spectral_solverPolarisation use math, only: & math_I3 use DAMASK_spectral_utilities, only: & - tSolutionState + tSolutionState, & + tSolutionParams implicit none private @@ -24,14 +25,6 @@ module DAMASK_spectral_solverPolarisation !-------------------------------------------------------------------------------------------------- ! derived types - type tSolutionParams !< @todo use here the type definition for a full loadcase including mask - real(pReal), dimension(3,3) :: P_BC, rotation_BC - real(pReal) :: timeinc - real(pReal) :: timeincOld - real(pReal) :: temperature - real(pReal) :: density - end type tSolutionParams - type(tSolutionParams), private :: params real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal @@ -124,12 +117,9 @@ subroutine Polarisation_init(temperature) use DAMASK_interface, only: & getSolverJobName use DAMASK_spectral_Utilities, only: & - Utilities_init, & Utilities_constitutiveResponse, & Utilities_updateGamma, & - Utilities_updateIPcoords, & - grid1Red, & - wgt + Utilities_updateIPcoords use mesh, only: & gridLocal, & gridGlobal @@ -150,7 +140,6 @@ subroutine Polarisation_init(temperature) integer(pInt) :: proc character(len=1024) :: rankStr - call Utilities_init() if (worldrank == 0_pInt) then write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverPolarisation init -+>>>' write(6,'(a)') ' $Id$' @@ -169,6 +158,7 @@ subroutine Polarisation_init(temperature) !-------------------------------------------------------------------------------------------------- ! PETSc Init call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) + call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) allocate(localK(worldsize), source = 0); localK(worldrank+1) = gridLocal(3) do proc = 1, worldsize call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr) @@ -182,10 +172,10 @@ subroutine Polarisation_init(temperature) gridLocal (1),gridLocal (2),localK, & ! local grid da,ierr) ! handle, error CHKERRQ(ierr) + call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) call DMDASNESSetFunctionLocal(da,INSERT_VALUES,Polarisation_formResidual,dummy,ierr) CHKERRQ(ierr) - call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) call SNESSetConvergenceTest(snes,Polarisation_converged,dummy,PETSC_NULL_FUNCTION,ierr) CHKERRQ(ierr) call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr) @@ -265,7 +255,7 @@ end subroutine Polarisation_init !-------------------------------------------------------------------------------------------------- type(tSolutionState) function & Polarisation_solution(incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,temperature_bc, & - rotation_BC,density) + rotation_BC) use numerics, only: & update_gamma use math, only: & @@ -277,8 +267,6 @@ type(tSolutionState) function & use FEsolving, only: & restartWrite, & terminallyIll - use numerics, only: & - worldrank implicit none @@ -288,8 +276,7 @@ type(tSolutionState) function & timeinc, & !< increment in time for current solution timeinc_old, & !< increment in time of last increment loadCaseTime, & !< remaining time of current load case - temperature_bc, & - density + temperature_bc logical, intent(in) :: & guess type(tBoundaryCondition), intent(in) :: & @@ -324,7 +311,6 @@ type(tSolutionState) function & params%timeinc = timeinc params%timeincOld = timeinc_old params%temperature = temperature_bc - params%density = density !-------------------------------------------------------------------------------------------------- ! solve BVP @@ -363,16 +349,14 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) math_transpose33, & math_mul3333xx33, & math_invSym3333, & - math_mul33x33, & - PI + math_mul33x33 use DAMASK_spectral_Utilities, only: & wgt, & - field_realMPI, & - field_fourierMPI, & - Utilities_FFTforward, & - Utilities_fourierConvolution, & + tensorField_realMPI, & + utilities_FFTtensorForward, & + utilities_fourierGammaConvolution, & Utilities_inverseLaplace, & - Utilities_FFTbackward, & + utilities_FFTtensorBackward, & Utilities_constitutiveResponse, & Utilities_divergenceRMS, & Utilities_curlRMS @@ -444,9 +428,9 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! - field_realMPI = 0.0_pReal + tensorField_realMPI = 0.0_pReal do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt, gridLocal(1) - field_realMPI(1:3,1:3,i,j,k) = & + tensorField_realMPI(1:3,1:3,i,j,k) = & polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& polarAlpha*math_mul33x33(F(1:3,1:3,i,j,k), & math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3)) @@ -454,13 +438,13 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! doing convolution in Fourier space - call Utilities_FFTforward() - call Utilities_fourierConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) - call Utilities_FFTbackward() + call utilities_FFTtensorForward() + call utilities_fourierGammaConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) + call utilities_FFTtensorBackward() !-------------------------------------------------------------------------------------------------- ! constructing residual - residual_F_tau = polarBeta*F - field_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) + residual_F_tau = polarBeta*F - tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response @@ -472,11 +456,11 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! calculate divergence - field_realMPI = 0.0_pReal - field_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = residual_F - call Utilities_FFTforward() + tensorField_realMPI = 0.0_pReal + tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = residual_F + call utilities_FFTtensorForward() err_div = Utilities_divergenceRMS() - call Utilities_FFTbackward() + call utilities_FFTtensorBackward() !-------------------------------------------------------------------------------------------------- ! constructing residual @@ -492,11 +476,11 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! calculating curl - field_realMPI = 0.0_pReal - field_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = F - call Utilities_FFTforward() + tensorField_realMPI = 0.0_pReal + tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = F + call utilities_FFTtensorForward() err_curl = Utilities_curlRMS() - call Utilities_FFTbackward() + call utilities_FFTtensorBackward() end subroutine Polarisation_formResidual @@ -727,7 +711,6 @@ subroutine Polarisation_destroy() call VecDestroy(solution_vec,ierr); CHKERRQ(ierr) call SNESDestroy(snes,ierr); CHKERRQ(ierr) call DMDestroy(da,ierr); CHKERRQ(ierr) - call Utilities_destroy() end subroutine Polarisation_destroy diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 00a74ef5b..201f7d199 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -950,7 +950,7 @@ end subroutine constitutive_collectDotState !> @brief for constitutive models having an instantaneous change of state (so far, only nonlocal) !> will return false if delta state is not needed/supported by the constitutive model !-------------------------------------------------------------------------------------------------- -logical function constitutive_collectDeltaState(Tstar_v, Fe, ipc, ip, el) +subroutine constitutive_collectDeltaState(Tstar_v, Fe, ipc, ip, el) use prec, only: & pReal, & pLongInt @@ -998,26 +998,19 @@ logical function constitutive_collectDeltaState(Tstar_v, Fe, ipc, ip, el) call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) select case (phase_plasticity(material_phase(ipc,ip,el))) - case (PLASTICITY_NONLOCAL_ID) - constitutive_collectDeltaState = .true. call plastic_nonlocal_deltaState(Tstar_v,ip,el) - case default - constitutive_collectDeltaState = .false. end select do mySource = 1_pInt, phase_Nsources(material_phase(ipc,ip,el)) select case (phase_source(mySource,material_phase(ipc,ip,el))) case (SOURCE_damage_isoBrittle_ID) - constitutive_collectDeltaState = constitutive_collectDeltaState .and. .true. call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(ipc,ip,el), Fe, & ipc, ip, el) case (SOURCE_vacancy_irradiation_ID) - constitutive_collectDeltaState = constitutive_collectDeltaState .and. .true. call source_vacancy_irradiation_deltaState(ipc, ip, el) case (SOURCE_vacancy_thermalfluc_ID) - constitutive_collectDeltaState = constitutive_collectDeltaState .and. .true. call source_vacancy_thermalfluc_deltaState(ipc, ip, el) end select @@ -1033,7 +1026,7 @@ logical function constitutive_collectDeltaState(Tstar_v, Fe, ipc, ip, el) !$OMP END CRITICAL (debugTimingDeltaState) endif -end function constitutive_collectDeltaState +end subroutine constitutive_collectDeltaState !-------------------------------------------------------------------------------------------------- diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 8bcc539fc..019d1f96a 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -3337,6 +3337,8 @@ logical function crystallite_stateJump(g,i,e) debug_g use material, only: & plasticState, & + sourceState, & + phase_Nsources, & mappingConstitutive use constitutive, only: & constitutive_collectDeltaState @@ -3350,32 +3352,41 @@ logical function crystallite_stateJump(g,i,e) integer(pInt) :: & c, & p, & - mySizePlasticDotState + mySource, & + mySizePlasticDeltaState, & + mySizeSourceDeltaState c = mappingConstitutive(1,g,i,e) p = mappingConstitutive(2,g,i,e) - if (constitutive_collectDeltaState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe(1:3,1:3,g,i,e), g,i,e)) then - mySizePlasticDotState = plasticState(p)%sizeDotState - if( any(plasticState(p)%deltaState(:,c) /= plasticState(p)%deltaState(:,c))) then ! NaN occured in deltaState + call constitutive_collectDeltaState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe(1:3,1:3,g,i,e), g,i,e) + mySizePlasticDeltaState = plasticState(p)%sizeDeltaState + if( any(plasticState(p)%deltaState(:,c) /= plasticState(p)%deltaState(:,c))) then ! NaN occured in deltaState + crystallite_stateJump = .false. + return + endif + plasticState(p)%state(1:mySizePlasticDeltaState,c) = plasticState(p)%state(1:mySizePlasticDeltaState,c) + & + plasticState(p)%deltaState(1:mySizePlasticDeltaState,c) + do mySource = 1_pInt, phase_Nsources(p) + mySizeSourceDeltaState = sourceState(p)%p(mySource)%sizeDeltaState + if( any(sourceState(p)%p(mySource)%deltaState(:,c) /= sourceState(p)%p(mySource)%deltaState(:,c))) then ! NaN occured in deltaState crystallite_stateJump = .false. return endif - plasticState(p)%state(1:mySizePlasticDotState,c) = plasticState(p)%state(1:mySizePlasticDotState,c) + & - plasticState(p)%deltaState(1:mySizePlasticDotState,c) + sourceState(p)%p(mySource)%state(1:mySizeSourceDeltaState,c) = & + sourceState(p)%p(mySource)%state(1:mySizeSourceDeltaState,c) + & + sourceState(p)%p(mySource)%deltaState(1:mySizeSourceDeltaState,c) + enddo #ifndef _OPENMP - p = mappingConstitutive(2,g,i,e) - c = mappingConstitutive(1,g,i,e) - if (any(plasticState(p)%deltaState(1:mySizePlasticDotState,c) /= 0.0_pReal) & - .and. iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & - .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) & - .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then - write(6,'(a,i8,1x,i2,1x,i3, /)') '<< CRYST >> update state at el ip g ',e,i,g - write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> deltaState', plasticState(p)%deltaState(1:mySizePlasticDotState,c) - write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> new state', plasticState(p)%state (1:mySizePlasticDotState,c) - endif -#endif + if (any(plasticState(p)%deltaState(1:mySizePlasticDeltaState,c) /= 0.0_pReal) & + .and. iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & + .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) & + .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then + write(6,'(a,i8,1x,i2,1x,i3, /)') '<< CRYST >> update state at el ip g ',e,i,g + write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> deltaState', plasticState(p)%deltaState(1:mySizePlasticDeltaState,c) + write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> new state', plasticState(p)%state (1:mySizePlasticDeltaState,c) endif +#endif crystallite_stateJump = .true. diff --git a/code/damage_nonlocal.f90 b/code/damage_nonlocal.f90 index d0b327234..23828f654 100644 --- a/code/damage_nonlocal.f90 +++ b/code/damage_nonlocal.f90 @@ -253,6 +253,8 @@ end subroutine damage_nonlocal_getSourceAndItsTangent !> @brief returns homogenized non local damage diffusion tensor in reference configuration !-------------------------------------------------------------------------------------------------- function damage_nonlocal_getDiffusion33(ip,el) + use numerics, only: & + charLength use lattice, only: & lattice_DamageDiffusion33 use material, only: & @@ -280,6 +282,7 @@ function damage_nonlocal_getDiffusion33(ip,el) enddo damage_nonlocal_getDiffusion33 = & + charLength*charLength* & damage_nonlocal_getDiffusion33/ & homogenization_Ngrains(homog) diff --git a/code/plastic_disloKMC.f90 b/code/plastic_disloKMC.f90 index 5bb51c137..c277bba42 100644 --- a/code/plastic_disloKMC.f90 +++ b/code/plastic_disloKMC.f90 @@ -201,7 +201,7 @@ subroutine plastic_disloKMC_init(fileUnit) Nchunks_TwinSlip = 0_pInt, Nchunks_TwinTwin = 0_pInt, & Nchunks_SlipFamilies = 0_pInt, Nchunks_TwinFamilies = 0_pInt, Nchunks_nonSchmid = 0_pInt, & offset_slip, index_myFamily, index_otherFamily - integer(pInt) :: sizeState, sizeDotState + integer(pInt) :: sizeState, sizeDotState, sizeDeltaState integer(pInt) :: NofMyPhase character(len=65536) :: & tag = '', & @@ -658,12 +658,14 @@ subroutine plastic_disloKMC_init(fileUnit) ! allocate state arrays sizeDotState = int(size(plastic_disloKMC_listBasicSlipStates),pInt) * ns & + int(size(plastic_disloKMC_listBasicTwinStates),pInt) * nt + sizeDeltaState = sizeDotState sizeState = sizeDotState & + int(size(plastic_disloKMC_listDependentSlipStates),pInt) * ns & + int(size(plastic_disloKMC_listDependentTwinStates),pInt) * nt plasticState(phase)%sizeState = sizeState plasticState(phase)%sizeDotState = sizeDotState + plasticState(phase)%sizeDeltaState = sizeDeltaState plasticState(phase)%sizePostResults = plastic_disloKMC_sizePostResults(instance) plasticState(phase)%nSlip = plastic_disloKMC_totalNslip(instance) plasticState(phase)%nTwin = 0_pInt @@ -676,7 +678,7 @@ subroutine plastic_disloKMC_init(fileUnit) allocate(plasticState(phase)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) - allocate(plasticState(phase)%deltaState (sizeDotState,NofMyPhase), source=0.0_pReal) + allocate(plasticState(phase)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) if (any(numerics_integrator == 1_pInt)) then allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) diff --git a/code/plastic_disloUCLA.f90 b/code/plastic_disloUCLA.f90 index f0e168f4f..2522500d4 100644 --- a/code/plastic_disloUCLA.f90 +++ b/code/plastic_disloUCLA.f90 @@ -206,7 +206,7 @@ subroutine plastic_disloUCLA_init(fileUnit) Nchunks_TwinSlip = 0_pInt, Nchunks_TwinTwin = 0_pInt, & Nchunks_SlipFamilies = 0_pInt, Nchunks_TwinFamilies = 0_pInt, Nchunks_nonSchmid = 0_pInt, & offset_slip, index_myFamily, index_otherFamily - integer(pInt) :: sizeState, sizeDotState + integer(pInt) :: sizeState, sizeDotState, sizeDeltaState integer(pInt) :: NofMyPhase character(len=65536) :: & tag = '', & @@ -677,12 +677,14 @@ subroutine plastic_disloUCLA_init(fileUnit) ! allocate state arrays sizeDotState = int(size(plastic_disloUCLA_listBasicSlipStates),pInt) * ns & + int(size(plastic_disloUCLA_listBasicTwinStates),pInt) * nt + sizeDeltaState = 0_pInt sizeState = sizeDotState & + int(size(plastic_disloUCLA_listDependentSlipStates),pInt) * ns & + int(size(plastic_disloUCLA_listDependentTwinStates),pInt) * nt plasticState(phase)%sizeState = sizeState plasticState(phase)%sizeDotState = sizeDotState + plasticState(phase)%sizeDeltaState = sizeDeltaState plasticState(phase)%sizePostResults = plastic_disloUCLA_sizePostResults(instance) plasticState(phase)%nSlip = plastic_disloucla_totalNslip(instance) plasticState(phase)%nTwin = plastic_disloucla_totalNtwin(instance) @@ -695,7 +697,7 @@ subroutine plastic_disloUCLA_init(fileUnit) allocate(plasticState(phase)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) - allocate(plasticState(phase)%deltaState (sizeDotState,NofMyPhase), source=0.0_pReal) + allocate(plasticState(phase)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) if (any(numerics_integrator == 1_pInt)) then allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) diff --git a/code/plastic_dislotwin.f90 b/code/plastic_dislotwin.f90 index 34c7ed83c..bd2135f23 100644 --- a/code/plastic_dislotwin.f90 +++ b/code/plastic_dislotwin.f90 @@ -233,7 +233,7 @@ subroutine plastic_dislotwin_init(fileUnit) Nchunks_TwinSlip = 0_pInt, Nchunks_TwinTwin = 0_pInt, & Nchunks_SlipFamilies = 0_pInt, Nchunks_TwinFamilies = 0_pInt, Nchunks_TransFamilies = 0_pInt, & offset_slip, index_myFamily, index_otherFamily - integer(pInt) :: sizeState, sizeDotState + integer(pInt) :: sizeState, sizeDotState, sizeDeltaState integer(pInt) :: NofMyPhase character(len=65536) :: & tag = '', & @@ -822,12 +822,14 @@ subroutine plastic_dislotwin_init(fileUnit) sizeDotState = int(size(plastic_dislotwin_listBasicSlipStates),pInt) * ns & + int(size(plastic_dislotwin_listBasicTwinStates),pInt) * nt & + int(size(plastic_dislotwin_listBasicTransStates),pInt) * nr + sizeDeltaState = 0_pInt sizeState = sizeDotState & + int(size(plastic_dislotwin_listDependentSlipStates),pInt) * ns & + int(size(plastic_dislotwin_listDependentTwinStates),pInt) * nt plasticState(phase)%sizeState = sizeState plasticState(phase)%sizeDotState = sizeDotState + plasticState(phase)%sizeDeltaState = sizeDeltaState plasticState(phase)%sizePostResults = plastic_dislotwin_sizePostResults(instance) plasticState(phase)%nSlip = plastic_dislotwin_totalNslip(instance) plasticState(phase)%nTwin = plastic_dislotwin_totalNtwin(instance) @@ -840,7 +842,7 @@ subroutine plastic_dislotwin_init(fileUnit) allocate(plasticState(phase)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) - allocate(plasticState(phase)%deltaState (sizeDotState,NofMyPhase), source=0.0_pReal) + allocate(plasticState(phase)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) if (any(numerics_integrator == 1_pInt)) then allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) diff --git a/code/plastic_j2.f90 b/code/plastic_j2.f90 index db4e13f9f..4a3ec1141 100644 --- a/code/plastic_j2.f90 +++ b/code/plastic_j2.f90 @@ -141,7 +141,8 @@ subroutine plastic_j2_init(fileUnit) instance, & mySize, & sizeDotState, & - sizeState + sizeState, & + sizeDeltaState character(len=65536) :: & tag = '', & line = '' @@ -324,8 +325,10 @@ subroutine plastic_j2_init(fileUnit) ! allocate state arrays sizeState = 2_pInt sizeDotState = sizeState + sizeDeltaState = 0_pInt plasticState(phase)%sizeState = sizeState plasticState(phase)%sizeDotState = sizeDotState + plasticState(phase)%sizeDeltaState = sizeDeltaState plasticState(phase)%sizePostResults = plastic_j2_sizePostResults(instance) plasticState(phase)%nSlip = 1 plasticState(phase)%nTwin = 0 @@ -341,6 +344,7 @@ subroutine plastic_j2_init(fileUnit) allocate(plasticState(phase)%state ( sizeState,NofMyPhase),source=0.0_pReal) allocate(plasticState(phase)%state_backup ( sizeState,NofMyPhase),source=0.0_pReal) allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase),source=0.0_pReal) + allocate(plasticState(phase)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%dotState_backup (sizeDotState,NofMyPhase),source=0.0_pReal) if (any(numerics_integrator == 1_pInt)) then allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase),source=0.0_pReal) diff --git a/code/plastic_none.f90 b/code/plastic_none.f90 index db7bbeec5..667c3b038 100644 --- a/code/plastic_none.f90 +++ b/code/plastic_none.f90 @@ -52,7 +52,8 @@ subroutine plastic_none_init phase, & NofMyPhase, & sizeState, & - sizeDotState + sizeDotState, & + sizeDeltaState mainProcess: if (worldrank == 0) then write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_NONE_label//' init -+>>>' @@ -75,6 +76,8 @@ subroutine plastic_none_init plasticState(phase)%sizeState = sizeState sizeDotState = sizeState plasticState(phase)%sizeDotState = sizeDotState + sizeDeltaState = 0_pInt + plasticState(phase)%sizeDeltaState = sizeDeltaState plasticState(phase)%sizePostResults = 0_pInt plasticState(phase)%nSlip = 0_pInt plasticState(phase)%nTwin = 0_pInt @@ -87,6 +90,7 @@ subroutine plastic_none_init allocate(plasticState(phase)%state_backup (sizeState,NofMyPhase)) allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase)) + allocate(plasticState(phase)%deltaState (sizeDeltaState,NofMyPhase)) allocate(plasticState(phase)%dotState_backup (sizeDotState,NofMyPhase)) if (any(numerics_integrator == 1_pInt)) then allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase)) diff --git a/code/plastic_nonlocal.f90 b/code/plastic_nonlocal.f90 index 69de63fc5..ce3336272 100644 --- a/code/plastic_nonlocal.f90 +++ b/code/plastic_nonlocal.f90 @@ -330,7 +330,7 @@ integer(pInt) :: phase, & tag = '', & line = '' - integer(pInt) :: sizeState, sizeDotState,sizeDependentState + integer(pInt) :: sizeState, sizeDotState,sizeDependentState, sizeDeltaState integer(pInt) :: NofMyPhase @@ -1144,6 +1144,7 @@ allocate(nonSchmidProjection(3,3,4,maxTotalNslip,maxNinstances), sizeDependentState = int(size(DEPENDENTSTATES),pInt) * ns sizeState = sizeDotState + sizeDependentState & + int(size(OTHERSTATES),pInt) * ns + sizeDeltaState = sizeDotState !*** determine indices to state array @@ -1299,6 +1300,7 @@ allocate(nonSchmidProjection(3,3,4,maxTotalNslip,maxNinstances), plasticState(phase)%sizeState = sizeState plasticState(phase)%sizeDotState = sizeDotState + plasticState(phase)%sizeDeltaState = sizeDeltaState plasticState(phase)%sizePostResults = plastic_nonlocal_sizePostResults(instance) plasticState(phase)%nonlocal = .true. plasticState(phase)%nSlip = totalNslip(instance) @@ -1312,7 +1314,7 @@ allocate(nonSchmidProjection(3,3,4,maxTotalNslip,maxNinstances), allocate(plasticState(phase)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) - allocate(plasticState(phase)%deltaState (sizeDotState,NofMyPhase), source=0.0_pReal) + allocate(plasticState(phase)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) if (any(numerics_integrator == 1_pInt)) then allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) diff --git a/code/plastic_phenopowerlaw.f90 b/code/plastic_phenopowerlaw.f90 index d4b48a01c..1260483f6 100644 --- a/code/plastic_phenopowerlaw.f90 +++ b/code/plastic_phenopowerlaw.f90 @@ -158,7 +158,7 @@ subroutine plastic_phenopowerlaw_init(fileUnit) Nchunks_TransFamilies = 0_pInt, Nchunks_nonSchmid = 0_pInt, & NipcMyPhase, & offset_slip, index_myFamily, index_otherFamily, & - mySize=0_pInt,sizeState,sizeDotState + mySize=0_pInt,sizeState,sizeDotState, sizeDeltaState character(len=65536) :: & tag = '', & line = '' @@ -552,8 +552,10 @@ subroutine plastic_phenopowerlaw_init(fileUnit) + plastic_phenopowerlaw_totalNtwin(instance) ! accshear_twin sizeDotState = sizeState + sizeDeltaState = 0_pInt plasticState(phase)%sizeState = sizeState plasticState(phase)%sizeDotState = sizeDotState + plasticState(phase)%sizeDeltaState = sizeDeltaState plasticState(phase)%sizePostResults = plastic_phenopowerlaw_sizePostResults(instance) plasticState(phase)%nSlip =plastic_phenopowerlaw_totalNslip(instance) plasticState(phase)%nTwin =plastic_phenopowerlaw_totalNtwin(instance) @@ -565,6 +567,7 @@ subroutine plastic_phenopowerlaw_init(fileUnit) allocate(plasticState(phase)%state ( sizeState,NipcMyPhase), source=0.0_pReal) allocate(plasticState(phase)%state_backup ( sizeState,NipcMyPhase), source=0.0_pReal) allocate(plasticState(phase)%dotState (sizeDotState,NipcMyPhase), source=0.0_pReal) + allocate(plasticState(phase)%deltaState (sizeDeltaState,NipcMyPhase), source=0.0_pReal) allocate(plasticState(phase)%dotState_backup (sizeDotState,NipcMyPhase), source=0.0_pReal) if (any(numerics_integrator == 1_pInt)) then allocate(plasticState(phase)%previousDotState (sizeDotState,NipcMyPhase),source=0.0_pReal) diff --git a/code/plastic_titanmod.f90 b/code/plastic_titanmod.f90 index 895a117c4..f563ba312 100644 --- a/code/plastic_titanmod.f90 +++ b/code/plastic_titanmod.f90 @@ -237,7 +237,7 @@ subroutine plastic_titanmod_init(fileUnit) Nchunks_SlipFamilies = 0_pInt, Nchunks_TwinFamilies = 0_pInt, & offset_slip, mySize, & maxTotalNslip,maxTotalNtwin, maxNinstance - integer(pInt) :: sizeState, sizeDotState + integer(pInt) :: sizeState, sizeDotState, sizeDeltaState integer(pInt) :: NofMyPhase = 0_pInt character(len=65536) :: & tag = '', & @@ -812,6 +812,7 @@ subroutine plastic_titanmod_init(fileUnit) sizeState = sizeDotState+ & size(plastic_titanmod_listDependentSlipStates)*ns + & size(plastic_titanmod_listDependentTwinStates)*nt + sizeDeltaState = 0_pInt !-------------------------------------------------------------------------------------------------- ! determine size of postResults array @@ -848,6 +849,7 @@ subroutine plastic_titanmod_init(fileUnit) ! Determine size of state array plasticState(phase)%sizeState = sizeState plasticState(phase)%sizeDotState = sizeDotState + plasticState(phase)%sizeDeltaState = sizeDeltaState plasticState(phase)%sizePostResults = plastic_titanmod_sizePostResults(instance) plasticState(phase)%nSlip =plastic_titanmod_totalNslip(instance) plasticState(phase)%nTwin = 0_pInt @@ -860,7 +862,7 @@ subroutine plastic_titanmod_init(fileUnit) allocate(plasticState(phase)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) - allocate(plasticState(phase)%deltaState (sizeDotState,NofMyPhase), source=0.0_pReal) + allocate(plasticState(phase)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) if (any(numerics_integrator == 1_pInt)) then allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) diff --git a/code/prec.f90 b/code/prec.f90 index 97252c60e..fac90015b 100644 --- a/code/prec.f90 +++ b/code/prec.f90 @@ -70,6 +70,7 @@ use ifport integer(pInt) :: & sizeState = 0_pInt , & !< size of state sizeDotState = 0_pInt, & !< size of dot state, i.e. parts of the state that are integrated + sizeDeltaState = 0_pInt, & !< size of delta state, i.e. parts of the state that have discontinuous rates sizePostResults = 0_pInt !< size of output data real(pReal), allocatable, dimension(:) :: & atolState diff --git a/code/source_damage_anisoBrittle.f90 b/code/source_damage_anisoBrittle.f90 index 6b05876cd..99289411e 100644 --- a/code/source_damage_anisoBrittle.f90 +++ b/code/source_damage_anisoBrittle.f90 @@ -106,7 +106,7 @@ subroutine source_damage_anisoBrittle_init(fileUnit) integer(pInt), parameter :: MAXNCHUNKS = 7_pInt integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions integer(pInt) :: maxNinstance,mySize=0_pInt,phase,instance,source,sourceOffset,o - integer(pInt) :: sizeState, sizeDotState + integer(pInt) :: sizeState, sizeDotState, sizeDeltaState integer(pInt) :: NofMyPhase integer(pInt) :: Nchunks_CleavageFamilies = 0_pInt, j character(len=65536) :: & @@ -255,10 +255,12 @@ subroutine source_damage_anisoBrittle_init(fileUnit) !-------------------------------------------------------------------------------------------------- ! Determine size of state array sizeDotState = 1_pInt + sizeDeltaState = 0_pInt sizeState = 1_pInt - sourceState(phase)%p(sourceOffset)%sizeState = sizeState - sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState + sourceState(phase)%p(sourceOffset)%sizeState = sizeState + sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState + sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState sourceState(phase)%p(sourceOffset)%sizePostResults = source_damage_anisoBrittle_sizePostResults(instance) allocate(sourceState(phase)%p(sourceOffset)%aTolState (sizeState), & source=source_damage_anisoBrittle_aTol(instance)) @@ -269,7 +271,7 @@ subroutine source_damage_anisoBrittle_init(fileUnit) allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDotState,NofMyPhase), source=0.0_pReal) + allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) if (any(numerics_integrator == 1_pInt)) then allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) diff --git a/code/source_damage_anisoDuctile.f90 b/code/source_damage_anisoDuctile.f90 index b78bcc08d..213eb8ebd 100644 --- a/code/source_damage_anisoDuctile.f90 +++ b/code/source_damage_anisoDuctile.f90 @@ -110,7 +110,7 @@ subroutine source_damage_anisoDuctile_init(fileUnit) integer(pInt), parameter :: MAXNCHUNKS = 7_pInt integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions integer(pInt) :: maxNinstance,mySize=0_pInt,phase,instance,source,sourceOffset,o - integer(pInt) :: sizeState, sizeDotState + integer(pInt) :: sizeState, sizeDotState, sizeDeltaState integer(pInt) :: NofMyPhase integer(pInt) :: Nchunks_SlipFamilies = 0_pInt, j character(len=65536) :: & @@ -258,9 +258,11 @@ subroutine source_damage_anisoDuctile_init(fileUnit) !-------------------------------------------------------------------------------------------------- ! Determine size of state array sizeDotState = 1_pInt + sizeDeltaState = 0_pInt sizeState = 1_pInt - sourceState(phase)%p(sourceOffset)%sizeState = sizeState - sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState + sourceState(phase)%p(sourceOffset)%sizeState = sizeState + sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState + sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState sourceState(phase)%p(sourceOffset)%sizePostResults = source_damage_anisoDuctile_sizePostResults(instance) allocate(sourceState(phase)%p(sourceOffset)%aTolState (sizeState), & source=source_damage_anisoDuctile_aTol(instance)) @@ -271,7 +273,7 @@ subroutine source_damage_anisoDuctile_init(fileUnit) allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDotState,NofMyPhase), source=0.0_pReal) + allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) if (any(numerics_integrator == 1_pInt)) then allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) diff --git a/code/source_damage_isoBrittle.f90 b/code/source_damage_isoBrittle.f90 index cf560e1f4..36227b707 100644 --- a/code/source_damage_isoBrittle.f90 +++ b/code/source_damage_isoBrittle.f90 @@ -92,7 +92,7 @@ subroutine source_damage_isoBrittle_init(fileUnit) integer(pInt), parameter :: MAXNCHUNKS = 7_pInt integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions integer(pInt) :: maxNinstance,mySize=0_pInt,phase,instance,source,sourceOffset,o - integer(pInt) :: sizeState, sizeDotState + integer(pInt) :: sizeState, sizeDotState, sizeDeltaState integer(pInt) :: NofMyPhase character(len=65536) :: & tag = '', & @@ -204,10 +204,12 @@ subroutine source_damage_isoBrittle_init(fileUnit) enddo outputsLoop ! Determine size of state array sizeDotState = 1_pInt + sizeDeltaState = 1_pInt sizeState = 1_pInt - sourceState(phase)%p(sourceOffset)%sizeState = sizeState - sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState + sourceState(phase)%p(sourceOffset)%sizeState = sizeState + sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState + sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState sourceState(phase)%p(sourceOffset)%sizePostResults = source_damage_isoBrittle_sizePostResults(instance) allocate(sourceState(phase)%p(sourceOffset)%aTolState (sizeState), & source=source_damage_isoBrittle_aTol(instance)) @@ -218,7 +220,7 @@ subroutine source_damage_isoBrittle_init(fileUnit) allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDotState,NofMyPhase), source=0.0_pReal) + allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) if (any(numerics_integrator == 1_pInt)) then allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) diff --git a/code/source_damage_isoDuctile.f90 b/code/source_damage_isoDuctile.f90 index af547aaa7..39d5d72bc 100644 --- a/code/source_damage_isoDuctile.f90 +++ b/code/source_damage_isoDuctile.f90 @@ -93,7 +93,7 @@ subroutine source_damage_isoDuctile_init(fileUnit) integer(pInt), parameter :: MAXNCHUNKS = 7_pInt integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions integer(pInt) :: maxNinstance,mySize=0_pInt,phase,instance,source,sourceOffset,o - integer(pInt) :: sizeState, sizeDotState + integer(pInt) :: sizeState, sizeDotState, sizeDeltaState integer(pInt) :: NofMyPhase character(len=65536) :: & tag = '', & @@ -208,11 +208,13 @@ subroutine source_damage_isoDuctile_init(fileUnit) endif enddo outputsLoop ! Determine size of state array - sizeDotState = 0_pInt + sizeDotState = 1_pInt + sizeDeltaState = 0_pInt sizeState = 1_pInt - sourceState(phase)%p(sourceOffset)%sizeState = sizeState - sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState + sourceState(phase)%p(sourceOffset)%sizeState = sizeState + sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState + sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState sourceState(phase)%p(sourceOffset)%sizePostResults = source_damage_isoDuctile_sizePostResults(instance) allocate(sourceState(phase)%p(sourceOffset)%aTolState (sizeState), & source=source_damage_isoDuctile_aTol(instance)) @@ -223,7 +225,7 @@ subroutine source_damage_isoDuctile_init(fileUnit) allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDotState,NofMyPhase), source=0.0_pReal) + allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) if (any(numerics_integrator == 1_pInt)) then allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) diff --git a/code/source_thermal_dissipation.f90 b/code/source_thermal_dissipation.f90 index 04e8c8690..34ed9a749 100644 --- a/code/source_thermal_dissipation.f90 +++ b/code/source_thermal_dissipation.f90 @@ -79,7 +79,7 @@ subroutine source_thermal_dissipation_init(fileUnit) integer(pInt), parameter :: MAXNCHUNKS = 7_pInt integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset - integer(pInt) :: sizeState, sizeDotState + integer(pInt) :: sizeState, sizeDotState, sizeDeltaState integer(pInt) :: NofMyPhase character(len=65536) :: & tag = '', & @@ -152,9 +152,11 @@ subroutine source_thermal_dissipation_init(fileUnit) sourceOffset = source_thermal_dissipation_offset(phase) sizeDotState = 0_pInt + sizeDeltaState = 0_pInt sizeState = 0_pInt sourceState(phase)%p(sourceOffset)%sizeState = sizeState sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState + sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState sourceState(phase)%p(sourceOffset)%sizePostResults = source_thermal_dissipation_sizePostResults(instance) allocate(sourceState(phase)%p(sourceOffset)%aTolState (sizeState), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%state0 (sizeState,NofMyPhase), source=0.0_pReal) @@ -164,7 +166,7 @@ subroutine source_thermal_dissipation_init(fileUnit) allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDotState,NofMyPhase), source=0.0_pReal) + allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) if (any(numerics_integrator == 1_pInt)) then allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) diff --git a/code/source_vacancy_irradiation.f90 b/code/source_vacancy_irradiation.f90 index 07798df29..ae1dd02fa 100644 --- a/code/source_vacancy_irradiation.f90 +++ b/code/source_vacancy_irradiation.f90 @@ -81,7 +81,7 @@ subroutine source_vacancy_irradiation_init(fileUnit) integer(pInt), parameter :: MAXNCHUNKS = 7_pInt integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset - integer(pInt) :: sizeState, sizeDotState + integer(pInt) :: sizeState, sizeDotState, sizeDeltaState integer(pInt) :: NofMyPhase character(len=65536) :: & tag = '', & @@ -158,9 +158,11 @@ subroutine source_vacancy_irradiation_init(fileUnit) sourceOffset = source_vacancy_irradiation_offset(phase) sizeDotState = 2_pInt + sizeDeltaState = 2_pInt sizeState = 2_pInt sourceState(phase)%p(sourceOffset)%sizeState = sizeState sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState + sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState sourceState(phase)%p(sourceOffset)%sizePostResults = source_vacancy_irradiation_sizePostResults(instance) allocate(sourceState(phase)%p(sourceOffset)%aTolState (sizeState), source=0.1_pReal) allocate(sourceState(phase)%p(sourceOffset)%state0 (sizeState,NofMyPhase), source=0.0_pReal) @@ -170,7 +172,7 @@ subroutine source_vacancy_irradiation_init(fileUnit) allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDotState,NofMyPhase), source=0.0_pReal) + allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) if (any(numerics_integrator == 1_pInt)) then allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) diff --git a/code/source_vacancy_phenoplasticity.f90 b/code/source_vacancy_phenoplasticity.f90 index 8fca355ea..f780a1d2c 100644 --- a/code/source_vacancy_phenoplasticity.f90 +++ b/code/source_vacancy_phenoplasticity.f90 @@ -79,7 +79,7 @@ subroutine source_vacancy_phenoplasticity_init(fileUnit) integer(pInt), parameter :: MAXNCHUNKS = 7_pInt integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset - integer(pInt) :: sizeState, sizeDotState + integer(pInt) :: sizeState, sizeDotState, sizeDeltaState integer(pInt) :: NofMyPhase character(len=65536) :: & tag = '', & @@ -152,9 +152,11 @@ subroutine source_vacancy_phenoplasticity_init(fileUnit) sourceOffset = source_vacancy_phenoplasticity_offset(phase) sizeDotState = 0_pInt + sizeDeltaState = 0_pInt sizeState = 0_pInt sourceState(phase)%p(sourceOffset)%sizeState = sizeState sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState + sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState sourceState(phase)%p(sourceOffset)%sizePostResults = source_vacancy_phenoplasticity_sizePostResults(instance) allocate(sourceState(phase)%p(sourceOffset)%aTolState (sizeState), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%state0 (sizeState,NofMyPhase), source=0.0_pReal) @@ -164,7 +166,7 @@ subroutine source_vacancy_phenoplasticity_init(fileUnit) allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDotState,NofMyPhase), source=0.0_pReal) + allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) if (any(numerics_integrator == 1_pInt)) then allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) diff --git a/code/source_vacancy_thermalfluc.f90 b/code/source_vacancy_thermalfluc.f90 index 3bfc05884..1150d1567 100644 --- a/code/source_vacancy_thermalfluc.f90 +++ b/code/source_vacancy_thermalfluc.f90 @@ -80,7 +80,7 @@ subroutine source_vacancy_thermalfluc_init(fileUnit) integer(pInt), parameter :: MAXNCHUNKS = 7_pInt integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset - integer(pInt) :: sizeState, sizeDotState + integer(pInt) :: sizeState, sizeDotState, sizeDeltaState integer(pInt) :: NofMyPhase character(len=65536) :: & tag = '', & @@ -153,9 +153,11 @@ subroutine source_vacancy_thermalfluc_init(fileUnit) sourceOffset = source_vacancy_thermalfluc_offset(phase) sizeDotState = 1_pInt + sizeDeltaState = 1_pInt sizeState = 1_pInt sourceState(phase)%p(sourceOffset)%sizeState = sizeState sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState + sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState sourceState(phase)%p(sourceOffset)%sizePostResults = source_vacancy_thermalfluc_sizePostResults(instance) allocate(sourceState(phase)%p(sourceOffset)%aTolState (sizeState), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%state0 (sizeState,NofMyPhase), source=0.0_pReal) @@ -165,7 +167,7 @@ subroutine source_vacancy_thermalfluc_init(fileUnit) allocate(sourceState(phase)%p(sourceOffset)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDotState,NofMyPhase), source=0.0_pReal) + allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) if (any(numerics_integrator == 1_pInt)) then allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)