diff --git a/src/IO.f90 b/src/IO.f90 index ba33a8a5c..9e0b4c40d 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -548,6 +548,8 @@ subroutine IO_error(error_ID,ext_msg,label1,ID1,label2,ID2) !-------------------------------------------------------------------------------------------------- ! user errors + case (600) + msg = 'only one source entry allowed' case (603) msg = 'invalid data for table' case (610) diff --git a/src/Marc/DAMASK_Marc.f90 b/src/Marc/DAMASK_Marc.f90 index e104ca345..ec6f34923 100644 --- a/src/Marc/DAMASK_Marc.f90 +++ b/src/Marc/DAMASK_Marc.f90 @@ -153,11 +153,10 @@ end module DAMASK_interface #include "../phase_mechanical_plastic_dislotungsten.f90" #include "../phase_mechanical_plastic_nonlocal.f90" #include "../phase_mechanical_eigen.f90" -#include "../phase_mechanical_eigen_cleavageopening.f90" #include "../phase_mechanical_eigen_thermalexpansion.f90" #include "../phase_thermal.f90" -#include "../phase_thermal_dissipation.f90" -#include "../phase_thermal_externalheat.f90" +#include "../phase_thermal_source_dissipation.f90" +#include "../phase_thermal_source_externalheat.f90" #include "../phase_damage.f90" #include "../phase_damage_isobrittle.f90" #include "../phase_damage_anisobrittle.f90" diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index d52a6f211..3cf378d7c 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -75,7 +75,7 @@ program DAMASK_grid cutBack = .false.,& sig integer :: & - i, j, m, field, & + i, j, field, & errorID = 0, & cutBackLevel = 0, & !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$ stepFraction = 0, & !< fraction of current time interval @@ -110,15 +110,8 @@ program DAMASK_grid load, & num_solver, & num_grid, & - load_step, & - solver, & - step_bc, & - step_mech, & - step_discretization + solver type(tList), pointer :: & -#ifdef __INTEL_LLVM_COMPILER - tensor, & -#endif load_steps character(len=:), allocatable :: & fileContent, fname @@ -210,13 +203,251 @@ program DAMASK_grid ID(field) = FIELD_DAMAGE_ID end if damageActive +!-------------------------------------------------------------------------------------------------- +! doing initialization depending on active solvers + call spectral_utilities_init() + + do field = 2, nActiveFields + select case (ID(field)) + + case (FIELD_THERMAL_ID) + call grid_thermal_spectral_init(num_grid) + + case (FIELD_DAMAGE_ID) + call grid_damage_spectral_init(num_grid) + + end select + end do + + call mechanical_init(num_grid) + call config_numerics_deallocate() !-------------------------------------------------------------------------------------------------- - load_steps => load%get_list('loadstep') - allocate(loadCases(load_steps%length)) ! array of load cases +! write header of output file + if (worldrank == 0) then + writeHeader: if (CLI_restartInc < 1) then + open(newunit=statUnit,file=trim(getSolverJobName())//'.sta',form='FORMATTED',status='REPLACE') + write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded StagIterationsNeeded' ! statistics file + else writeHeader + open(newunit=statUnit,file=trim(getSolverJobName())//& + '.sta',form='FORMATTED', position='APPEND', status='OLD') + end if writeHeader + end if + writeUndeformed: if (CLI_restartInc < 1) then + print'(/,1x,a)', '... saving initial configuration ..........................................' + flush(IO_STDOUT) + call materialpoint_result(0,0.0_pREAL) + end if writeUndeformed + + loadCases = parseLoadSteps(load%get_list('loadstep')) + + loadCaseLooping: do l = 1, size(loadCases) + t_0 = t ! load case start time + guess = loadCases(l)%estimate_rate ! change of load case? homogeneous guess for the first inc + + incLooping: do inc = 1, loadCases(l)%N + totalIncsCounter = totalIncsCounter + 1 + +!-------------------------------------------------------------------------------------------------- +! forwarding time + Delta_t_prev = Delta_t ! last time intervall that brought former inc to an end + if (dEq(loadCases(l)%r,1.0_pREAL,1.e-9_pREAL)) then ! linear scale + Delta_t = loadCases(l)%t/real(loadCases(l)%N,pREAL) + else + Delta_t = loadCases(l)%t * (loadCases(l)%r**(inc-1)-loadCases(l)%r**inc) & + / (1.0_pREAL-loadCases(l)%r**loadCases(l)%N) + end if + Delta_t = Delta_t * real(subStepFactor,pREAL)**real(-cutBackLevel,pREAL) ! depending on cut back level, decrease time step + + skipping: if (totalIncsCounter <= CLI_restartInc) then ! not yet at restart inc? + t = t + Delta_t ! just advance time, skip already performed calculation + guess = .true. ! QUESTION:why forced guessing instead of inheriting loadcase preference + else skipping + stepFraction = 0 ! fraction scaled by stepFactor**cutLevel + + subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel) + t_remaining = loadCases(l)%t + t_0 - t + t = t + Delta_t ! forward target time + stepFraction = stepFraction + 1 ! count step + +!-------------------------------------------------------------------------------------------------- +! report beginning of new step + print'(/,1x,a)', '###########################################################################' + print'(1x,a,1x,es12.5,6(a,i0))', & + 'Time', t, & + 's: Increment ', inc,'/',loadCases(l)%N,& + '-', stepFraction,'/',subStepFactor**cutBackLevel,& + ' of load case ', l,'/',size(loadCases) + write(incInfo,'(4(a,i0))') & + 'Increment ',totalIncsCounter,'/',sum(loadCases%N),& + '-', stepFraction,'/',subStepFactor**cutBackLevel + flush(IO_STDOUT) + +!-------------------------------------------------------------------------------------------------- +! forward fields + do field = 1, nActiveFields + select case(ID(field)) + case(FIELD_MECH_ID) + call mechanical_forward (& + cutBack,guess,Delta_t,Delta_t_prev,t_remaining, & + deformation_BC = loadCases(l)%deformation, & + stress_BC = loadCases(l)%stress, & + rotation_BC = loadCases(l)%rot) + + case(FIELD_THERMAL_ID); call grid_thermal_spectral_forward(cutBack) + case(FIELD_DAMAGE_ID); call grid_damage_spectral_forward(cutBack) + end select + end do + if (.not. cutBack) call materialpoint_forward + +!-------------------------------------------------------------------------------------------------- +! solve fields + stagIter = 1 + stagIterate = .true. + do while (stagIterate) + + if (nActiveFields > 1) print'(/,1x,a,i0)', 'Staggered Iteration ',stagIter + do field = 1, nActiveFields + select case(ID(field)) + case(FIELD_MECH_ID) + solres(field) = mechanical_solution(incInfo) + case(FIELD_THERMAL_ID) + solres(field) = grid_thermal_spectral_solution(Delta_t) + case(FIELD_DAMAGE_ID) + solres(field) = grid_damage_spectral_solution(Delta_t) + end select + + if (.not. solres(field)%converged) exit ! no solution found + + end do + stagIter = stagIter + 1 + stagIterate = stagIter <= stagItMax & + .and. all(solres(:)%converged) & + .and. .not. all(solres(:)%stagConverged) ! stationary with respect to staggered iteration + end do + +!-------------------------------------------------------------------------------------------------- +! check solution and either advance or retry with smaller timestep + + if ( (all(solres(:)%converged .and. solres(:)%stagConverged)) & ! converged + .and. .not. solres(1)%termIll) then ! and acceptable solution found + call mechanical_updateCoords() + Delta_t_prev = Delta_t + cutBack = .false. + guess = .true. ! start guessing after first converged (sub)inc + if (worldrank == 0) then + write(statUnit,*) totalIncsCounter, t, cutBackLevel, & + solres(1)%converged, solres(1)%iterationsNeeded, StagIter + flush(statUnit) + end if + elseif (cutBackLevel < maxCutBack) then ! further cutbacking tolerated? + cutBack = .true. + stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator + cutBackLevel = cutBackLevel + 1 + t = t - Delta_t + Delta_t = Delta_t/real(subStepFactor,pREAL) ! cut timestep + print'(/,1x,a)', 'cutting back ' + else ! no more options to continue + if (worldrank == 0) close(statUnit) + call IO_error(950) + end if + + end do subStepLooping + + cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc + + if (all(solres(:)%converged)) then + print'(/,1x,a,1x,i0,1x,a)', 'increment', totalIncsCounter, 'converged' + else + print'(/,1x,a,1x,i0,1x,a)', 'increment', totalIncsCounter, 'NOT converged' + end if; flush(IO_STDOUT) + + call MPI_Allreduce(signal_SIGUSR1,sig,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) + if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + if (mod(inc,loadCases(l)%f_out) == 0 .or. sig) then + print'(/,1x,a)', '... saving results ........................................................' + flush(IO_STDOUT) + call materialpoint_result(totalIncsCounter,t) + end if + if (sig) call signal_setSIGUSR1(.false.) + call MPI_Allreduce(signal_SIGUSR2,sig,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) + if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + if (mod(inc,loadCases(l)%f_restart) == 0 .or. sig) then + do field = 1, nActiveFields + select case (ID(field)) + case(FIELD_MECH_ID) + call mechanical_restartWrite() + case(FIELD_THERMAL_ID) + call grid_thermal_spectral_restartWrite() + case(FIELD_DAMAGE_ID) + call grid_damage_spectral_restartWrite() + end select + end do + call materialpoint_restartWrite() + end if + if (sig) call signal_setSIGUSR2(.false.) + call MPI_Allreduce(signal_SIGINT,sig,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) + if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + if (sig) exit loadCaseLooping + end if skipping + + end do incLooping + + end do loadCaseLooping + + +!-------------------------------------------------------------------------------------------------- +! report summary of whole calculation + print'(/,1x,a)', '###########################################################################' + if (worldrank == 0) close(statUnit) + + call quit(0) ! no complains ;) + + +contains + +subroutine getMaskedTensor(values,mask,tensor) + + real(pREAL), intent(out), dimension(3,3) :: values + logical, intent(out), dimension(3,3) :: mask + type(tList), pointer :: tensor + + type(tList), pointer :: row + integer :: i,j + + + values = 0.0_pREAL + do i = 1,3 + row => tensor%get_list(i) + do j = 1,3 + mask(i,j) = row%get_asStr(j) == 'x' + if (.not. mask(i,j)) values(i,j) = row%get_asReal(j) + end do + end do + +end subroutine getMaskedTensor + + +function parseLoadsteps(load_steps) result(loadCases) + + type(tList), intent(in), target :: load_steps + type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases + + integer :: l,m + type(tDict), pointer :: & + load_step, & + step_bc, & + step_mech, & + step_discretization +#ifdef __INTEL_LLVM_COMPILER + type(tList), pointer :: & + tensor +#endif + + + allocate(loadCases(load_steps%length)) do l = 1, load_steps%length - load_step => load_steps%get_dict(l) step_bc => load_step%get_dict('boundary_conditions') step_mech => step_bc%get_dict('mechanical') @@ -316,226 +547,6 @@ program DAMASK_grid end if reportAndCheck end do - -!-------------------------------------------------------------------------------------------------- -! doing initialization depending on active solvers - call spectral_utilities_init() - - do field = 2, nActiveFields - select case (ID(field)) - - case (FIELD_THERMAL_ID) - call grid_thermal_spectral_init(num_grid) - - case (FIELD_DAMAGE_ID) - call grid_damage_spectral_init(num_grid) - - end select - end do - - call mechanical_init(num_grid) - call config_numerics_deallocate() - -!-------------------------------------------------------------------------------------------------- -! write header of output file - if (worldrank == 0) then - writeHeader: if (CLI_restartInc < 1) then - open(newunit=statUnit,file=trim(getSolverJobName())//'.sta',form='FORMATTED',status='REPLACE') - write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded' ! statistics file - else writeHeader - open(newunit=statUnit,file=trim(getSolverJobName())//& - '.sta',form='FORMATTED', position='APPEND', status='OLD') - end if writeHeader - end if - - writeUndeformed: if (CLI_restartInc < 1) then - print'(/,1x,a)', '... saving initial configuration ..........................................' - flush(IO_STDOUT) - call materialpoint_result(0,0.0_pREAL) - end if writeUndeformed - - loadCaseLooping: do l = 1, size(loadCases) - t_0 = t ! load case start time - guess = loadCases(l)%estimate_rate ! change of load case? homogeneous guess for the first inc - - incLooping: do inc = 1, loadCases(l)%N - totalIncsCounter = totalIncsCounter + 1 - -!-------------------------------------------------------------------------------------------------- -! forwarding time - Delta_t_prev = Delta_t ! last time intervall that brought former inc to an end - if (dEq(loadCases(l)%r,1.0_pREAL,1.e-9_pREAL)) then ! linear scale - Delta_t = loadCases(l)%t/real(loadCases(l)%N,pREAL) - else - Delta_t = loadCases(l)%t * (loadCases(l)%r**(inc-1)-loadCases(l)%r**inc) & - / (1.0_pREAL-loadCases(l)%r**loadCases(l)%N) - end if - Delta_t = Delta_t * real(subStepFactor,pREAL)**real(-cutBackLevel,pREAL) ! depending on cut back level, decrease time step - - skipping: if (totalIncsCounter <= CLI_restartInc) then ! not yet at restart inc? - t = t + Delta_t ! just advance time, skip already performed calculation - guess = .true. ! QUESTION:why forced guessing instead of inheriting loadcase preference - else skipping - stepFraction = 0 ! fraction scaled by stepFactor**cutLevel - - subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel) - t_remaining = loadCases(l)%t + t_0 - t - t = t + Delta_t ! forward target time - stepFraction = stepFraction + 1 ! count step - -!-------------------------------------------------------------------------------------------------- -! report beginning of new step - print'(/,1x,a)', '###########################################################################' - print'(1x,a,1x,es12.5,6(a,i0))', & - 'Time', t, & - 's: Increment ', inc,'/',loadCases(l)%N,& - '-', stepFraction,'/',subStepFactor**cutBackLevel,& - ' of load case ', l,'/',size(loadCases) - write(incInfo,'(4(a,i0))') & - 'Increment ',totalIncsCounter,'/',sum(loadCases%N),& - '-', stepFraction,'/',subStepFactor**cutBackLevel - flush(IO_STDOUT) - -!-------------------------------------------------------------------------------------------------- -! forward fields - do field = 1, nActiveFields - select case(ID(field)) - case(FIELD_MECH_ID) - call mechanical_forward (& - cutBack,guess,Delta_t,Delta_t_prev,t_remaining, & - deformation_BC = loadCases(l)%deformation, & - stress_BC = loadCases(l)%stress, & - rotation_BC = loadCases(l)%rot) - - case(FIELD_THERMAL_ID); call grid_thermal_spectral_forward(cutBack) - case(FIELD_DAMAGE_ID); call grid_damage_spectral_forward(cutBack) - end select - end do - if (.not. cutBack) call materialpoint_forward - -!-------------------------------------------------------------------------------------------------- -! solve fields - stagIter = 0 - stagIterate = .true. - do while (stagIterate) - do field = 1, nActiveFields - select case(ID(field)) - case(FIELD_MECH_ID) - solres(field) = mechanical_solution(incInfo) - case(FIELD_THERMAL_ID) - solres(field) = grid_thermal_spectral_solution(Delta_t) - case(FIELD_DAMAGE_ID) - solres(field) = grid_damage_spectral_solution(Delta_t) - end select - - if (.not. solres(field)%converged) exit ! no solution found - - end do - stagIter = stagIter + 1 - stagIterate = stagIter < stagItMax & - .and. all(solres(:)%converged) & - .and. .not. all(solres(:)%stagConverged) ! stationary with respect to staggered iteration - end do - -!-------------------------------------------------------------------------------------------------- -! check solution and either advance or retry with smaller timestep - - if ( (all(solres(:)%converged .and. solres(:)%stagConverged)) & ! converged - .and. .not. solres(1)%termIll) then ! and acceptable solution found - call mechanical_updateCoords() - Delta_t_prev = Delta_t - cutBack = .false. - guess = .true. ! start guessing after first converged (sub)inc - if (worldrank == 0) then - write(statUnit,*) totalIncsCounter, t, cutBackLevel, & - solres(1)%converged, solres(1)%iterationsNeeded - flush(statUnit) - end if - elseif (cutBackLevel < maxCutBack) then ! further cutbacking tolerated? - cutBack = .true. - stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator - cutBackLevel = cutBackLevel + 1 - t = t - Delta_t - Delta_t = Delta_t/real(subStepFactor,pREAL) ! cut timestep - print'(/,1x,a)', 'cutting back ' - else ! no more options to continue - if (worldrank == 0) close(statUnit) - call IO_error(950) - end if - - end do subStepLooping - - cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc - - if (all(solres(:)%converged)) then - print'(/,1x,a,1x,i0,1x,a)', 'increment', totalIncsCounter, 'converged' - else - print'(/,1x,a,1x,i0,1x,a)', 'increment', totalIncsCounter, 'NOT converged' - end if; flush(IO_STDOUT) - - call MPI_Allreduce(signal_SIGUSR1,sig,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' - if (mod(inc,loadCases(l)%f_out) == 0 .or. sig) then - print'(/,1x,a)', '... saving results ........................................................' - flush(IO_STDOUT) - call materialpoint_result(totalIncsCounter,t) - end if - if (sig) call signal_setSIGUSR1(.false.) - call MPI_Allreduce(signal_SIGUSR2,sig,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' - if (mod(inc,loadCases(l)%f_restart) == 0 .or. sig) then - do field = 1, nActiveFields - select case (ID(field)) - case(FIELD_MECH_ID) - call mechanical_restartWrite() - case(FIELD_THERMAL_ID) - call grid_thermal_spectral_restartWrite() - case(FIELD_DAMAGE_ID) - call grid_damage_spectral_restartWrite() - end select - end do - call materialpoint_restartWrite() - end if - if (sig) call signal_setSIGUSR2(.false.) - call MPI_Allreduce(signal_SIGINT,sig,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' - if (sig) exit loadCaseLooping - end if skipping - - end do incLooping - - end do loadCaseLooping - - -!-------------------------------------------------------------------------------------------------- -! report summary of whole calculation - print'(/,1x,a)', '###########################################################################' - if (worldrank == 0) close(statUnit) - - call quit(0) ! no complains ;) - - -contains - -subroutine getMaskedTensor(values,mask,tensor) - - real(pREAL), intent(out), dimension(3,3) :: values - logical, intent(out), dimension(3,3) :: mask - type(tList), pointer :: tensor - - type(tList), pointer :: row - integer :: i,j - - - values = 0.0_pREAL - do i = 1,3 - row => tensor%get_list(i) - do j = 1,3 - mask(i,j) = row%get_asStr(j) == 'x' - if (.not. mask(i,j)) values(i,j) = row%get_asReal(j) - end do - end do - -end subroutine getMaskedTensor +end function parseLoadsteps end program DAMASK_grid diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 27295b76f..48023cc54 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -75,7 +75,6 @@ subroutine grid_damage_spectral_init(num_grid) type(tDict), pointer, intent(in) :: num_grid integer(MPI_INTEGER_KIND), dimension(0:worldsize-1) :: cells3_global - integer :: i, j, k, ce DM :: DM_damage real(pREAL), dimension(:,:,:), pointer :: phi ! 0-indexed Vec :: uBound, lBound @@ -92,7 +91,6 @@ subroutine grid_damage_spectral_init(num_grid) petsc_options - print'(/,1x,a)', '<<<+- grid_spectral_damage init -+>>>' print'(/,1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019' @@ -140,7 +138,7 @@ subroutine grid_damage_spectral_init(num_grid) 1_pPETSCINT, 1_pPETSCINT, int(worldsize,pPETSCINT), & 1_pPETSCINT, 0_pPETSCINT, & ! #dof (phi, scalar), ghost boundary width (domain overlap) [int(cells(1),pPetscInt)],[int(cells(2),pPetscInt)],int(cells3_global,pPETSCINT), & ! local cells - DM_damage,err_PETSc) ! handle, error + DM_damage,err_PETSc) ! handle, error CHKERRQ(err_PETSc) call DMsetFromOptions(DM_damage,err_PETSc) CHKERRQ(err_PETSc) @@ -194,11 +192,7 @@ subroutine grid_damage_spectral_init(num_grid) phi_stagInc = phi_lastInc end if restartRead - ce = 0 - do k = 0, cells3-1; do j = 0, cells(2)-1; do i = 0, cells(1)-1 - ce = ce + 1 - call homogenization_set_phi(phi(i,j,k),ce) - end do; end do; end do + call homogenization_set_phi(reshape(phi,[product(cells(1:2))*cells3])) call DMDAVecRestoreArrayF90(DM_damage,phi_PETSc,phi,err_PETSc) CHKERRQ(err_PETSc) @@ -216,7 +210,6 @@ function grid_damage_spectral_solution(Delta_t) result(solution) real(pREAL), intent(in) :: & Delta_t !< increment in time for current solution - integer :: i, j, k, ce type(tSolutionState) :: solution PetscInt :: devNull PetscReal :: phi_min, phi_max, stagNorm @@ -227,8 +220,6 @@ function grid_damage_spectral_solution(Delta_t) result(solution) SNESConvergedReason :: reason - solution%converged = .false. - !-------------------------------------------------------------------------------------------------- ! set module wide availabe data params%Delta_t = Delta_t @@ -256,13 +247,7 @@ function grid_damage_spectral_solution(Delta_t) result(solution) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' phi_stagInc = phi -!-------------------------------------------------------------------------------------------------- -! updating damage state - ce = 0 - do k = 0, cells3-1; do j = 0, cells(2)-1; do i = 0,cells(1)-1 - ce = ce + 1 - call homogenization_set_phi(phi(i,j,k),ce) - end do; end do; end do + call homogenization_set_phi(reshape(phi,[product(cells(1:2))*cells3])) call DMDAVecRestoreArrayF90(DM_damage,phi_PETSc,phi,err_PETSc) CHKERRQ(err_PETSc) @@ -283,7 +268,6 @@ subroutine grid_damage_spectral_forward(cutBack) logical, intent(in) :: cutBack - integer :: i, j, k, ce DM :: DM_damage real(pREAL), dimension(:,:,:), pointer :: phi ! 0-indexed PetscErrorCode :: err_PETSc @@ -295,11 +279,7 @@ subroutine grid_damage_spectral_forward(cutBack) CHKERRQ(err_PETSc) if (cutBack) then - ce = 0 - do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) - ce = ce + 1 - call homogenization_set_phi(phi_lastInc(i,j,k),ce) - end do; end do; end do + call homogenization_set_phi(reshape(phi_lastInc,[product(cells(1:2))*cells3])) phi = phi_lastInc phi_stagInc = phi_lastInc else diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 0db9db64b..ed1ff1e2f 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -365,7 +365,7 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_ if (stress_BC%myType=='dot_P') P_aim = P_aim & + merge(.0_pREAL,stress_BC%values,stress_BC%mask)*Delta_t - F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average + F = reshape(utilities_forwardTensorField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average rotation_BC%rotate(F_aim,active=.true.)),[9,cells(1),cells(2),cells3]) call DMDAVecRestoreArrayF90(DM_mech,F_PETSc,F,err_PETSc) CHKERRQ(err_PETSc) diff --git a/src/grid/grid_mech_spectral_polarization.f90 b/src/grid/grid_mech_spectral_polarization.f90 index 5342fa51a..c6c7c2fd5 100644 --- a/src/grid/grid_mech_spectral_polarization.f90 +++ b/src/grid/grid_mech_spectral_polarization.f90 @@ -408,11 +408,11 @@ subroutine grid_mechanical_spectral_polarization_forward(cutBack,guess,Delta_t,D if (stress_BC%myType=='dot_P') P_aim = P_aim & + merge(.0_pREAL,stress_BC%values,stress_BC%mask)*Delta_t - F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average - rotation_BC%rotate(F_aim,active=.true.)),& + F = reshape(utilities_forwardTensorField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average + rotation_BC%rotate(F_aim,active=.true.)),& [9,cells(1),cells(2),cells3]) if (guess) then - F_tau = reshape(Utilities_forwardField(Delta_t,F_tau_lastInc,F_taudot), & + F_tau = reshape(Utilities_forwardTensorField(Delta_t,F_tau_lastInc,F_taudot), & [9,cells(1),cells(2),cells3]) ! does not have any average value as boundary condition else do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1) diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index f6b5c19a0..86d8e04f3 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -74,7 +74,7 @@ subroutine grid_thermal_spectral_init(num_grid) type(tDict), pointer, intent(in) :: num_grid integer(MPI_INTEGER_KIND), dimension(0:worldsize-1) :: cells3_global - integer :: i, j, k, ce + integer :: ce DM :: DM_thermal real(pREAL), dimension(:,:,:), pointer :: T ! 0-indexed integer(MPI_INTEGER_KIND) :: err_MPI @@ -88,7 +88,6 @@ subroutine grid_thermal_spectral_init(num_grid) petsc_options - print'(/,1x,a)', '<<<+- grid_thermal_spectral init -+>>>' print'(/,1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019' @@ -171,11 +170,8 @@ subroutine grid_thermal_spectral_init(num_grid) dotT_lastInc = 0.0_pREAL * T_lastInc end if restartRead - ce = 0 - do k = 0, cells3-1; do j = 0, cells(2)-1; do i = 0, cells(1)-1 - ce = ce + 1 - call homogenization_thermal_setField(T(i,j,k),0.0_pREAL,ce) - end do; end do; end do + call homogenization_thermal_setField(reshape(T,[product(cells(1:2))*cells3]), & + [(0.0_pReal, ce = 1,product(cells(1:2))*cells3)]) call DMDAVecRestoreArrayF90(DM_thermal,T_PETSc,T,err_PETSc) CHKERRQ(err_PETSc) @@ -204,8 +200,6 @@ function grid_thermal_spectral_solution(Delta_t) result(solution) SNESConvergedReason :: reason - solution%converged = .false. - !-------------------------------------------------------------------------------------------------- ! set module wide availabe data params%Delta_t = Delta_t @@ -233,13 +227,8 @@ function grid_thermal_spectral_solution(Delta_t) result(solution) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' T_stagInc = T -!-------------------------------------------------------------------------------------------------- -! updating thermal state - ce = 0 - do k = 0, cells3-1; do j = 0, cells(2)-1; do i = 0, cells(1)-1 - ce = ce + 1 - call homogenization_thermal_setField(T(i,j,k),(T(i,j,k)-T_lastInc(i+1,j+1,k+1))/params%Delta_t,ce) - end do; end do; end do + call homogenization_thermal_setField(reshape(T,[product(cells(1:2))*cells3]), & + reshape(T-T_lastInc,[product(cells(1:2))*cells3])/params%Delta_t) call DMDAVecRestoreArrayF90(DM_thermal,T_PETSc,T,err_PETSc) CHKERRQ(err_PETSc) @@ -260,7 +249,6 @@ subroutine grid_thermal_spectral_forward(cutBack) logical, intent(in) :: cutBack - integer :: i, j, k, ce DM :: DM_thermal real(pREAL), dimension(:,:,:), pointer :: T ! 0-indexed PetscErrorCode :: err_PETSc @@ -272,11 +260,8 @@ subroutine grid_thermal_spectral_forward(cutBack) CHKERRQ(err_PETSc) if (cutBack) then - ce = 0 - do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) - ce = ce + 1 - call homogenization_thermal_setField(T_lastInc(i,j,k),dotT_lastInc(i,j,k),ce) - end do; end do; end do + call homogenization_thermal_setField(reshape(T_lastInc,[product(cells(1:2))*cells3]), & + reshape(dotT_lastInc,[product(cells(1:2))*cells3])) T = T_lastInc T_stagInc = T_lastInc else diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 453723a9b..03d945bee 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -124,7 +124,7 @@ module spectral_utilities utilities_maskedCompliance, & utilities_constitutiveResponse, & utilities_calculateRate, & - utilities_forwardField, & + utilities_forwardTensorField, & utilities_updateCoords contains @@ -864,7 +864,7 @@ end function utilities_calculateRate !> @brief forwards a field with a pointwise given rate, if aim is given, !> ensures that the average matches the aim !-------------------------------------------------------------------------------------------------- -function utilities_forwardField(Delta_t,field_lastInc,rate,aim) +function utilities_forwardTensorField(Delta_t,field_lastInc,rate,aim) real(pREAL), intent(in) :: & Delta_t !< Delta_t of current step @@ -875,22 +875,22 @@ function utilities_forwardField(Delta_t,field_lastInc,rate,aim) aim !< average field value aim real(pREAL), dimension(3,3,cells(1),cells(2),cells3) :: & - utilities_forwardField + utilities_forwardTensorField real(pREAL), dimension(3,3) :: fieldDiff !< - aim integer(MPI_INTEGER_KIND) :: err_MPI - utilities_forwardField = field_lastInc + rate*Delta_t + utilities_forwardTensorField = field_lastInc + rate*Delta_t if (present(aim)) then !< correct to match average - fieldDiff = sum(sum(sum(utilities_forwardField,dim=5),dim=4),dim=3)*wgt + fieldDiff = sum(sum(sum(utilities_forwardTensorField,dim=5),dim=4),dim=3)*wgt call MPI_Allreduce(MPI_IN_PLACE,fieldDiff,9_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' fieldDiff = fieldDiff - aim - utilities_forwardField = utilities_forwardField & - - spread(spread(spread(fieldDiff,3,cells(1)),4,cells(2)),5,cells3) + utilities_forwardTensorField = utilities_forwardTensorField & + - spread(spread(spread(fieldDiff,3,cells(1)),4,cells(2)),5,cells3) end if -end function utilities_forwardField +end function utilities_forwardTensorField !-------------------------------------------------------------------------------------------------- diff --git a/src/homogenization.f90 b/src/homogenization.f90 index c85547917..63eacdf9d 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -145,9 +145,8 @@ module homogenization real(pREAL) :: f end function homogenization_f_T - module subroutine homogenization_thermal_setField(T,dot_T, ce) - integer, intent(in) :: ce - real(pREAL), intent(in) :: T, dot_T + module subroutine homogenization_thermal_setField(T,dot_T) + real(pREAL), dimension(:), intent(in) :: T, dot_T end subroutine homogenization_thermal_setField module function homogenization_damage_active() result(active) @@ -170,10 +169,8 @@ module homogenization real(pREAL) :: f end function homogenization_f_phi - module subroutine homogenization_set_phi(phi,ce) - integer, intent(in) :: ce - real(pREAL), intent(in) :: & - phi + module subroutine homogenization_set_phi(phi) + real(pREAL), dimension(:), intent(in) :: phi end subroutine homogenization_set_phi end interface diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90 index 466b8b47b..7e83c34d8 100644 --- a/src/homogenization_damage.f90 +++ b/src/homogenization_damage.f90 @@ -151,20 +151,19 @@ end function homogenization_f_phi !-------------------------------------------------------------------------------------------------- !> @brief Set damage field. !-------------------------------------------------------------------------------------------------- -module subroutine homogenization_set_phi(phi,ce) +module subroutine homogenization_set_phi(phi) - integer, intent(in) :: ce - real(pREAL), intent(in) :: phi + real(pREAL), dimension(:), intent(in) :: phi - integer :: & - ho, & - en + integer :: ho, en, ce - ho = material_ID_homogenization(ce) - en = material_entry_homogenization(ce) - damagestate_h(ho)%state(1,en) = phi - current(ho)%phi(en) = phi + do ce=lbound(phi,1), ubound(phi,1) + ho = material_ID_homogenization(ce) + en = material_entry_homogenization(ce) + damagestate_h(ho)%state(1,en) = phi(ce) + current(ho)%phi(en) = phi(ce) + end do end subroutine homogenization_set_phi diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90 index 789ac994b..069c402eb 100644 --- a/src/homogenization_thermal.f90 +++ b/src/homogenization_thermal.f90 @@ -173,15 +173,20 @@ end function homogenization_f_T !-------------------------------------------------------------------------------------------------- !> @brief Set thermal field and its rate (T and dot_T). !-------------------------------------------------------------------------------------------------- -module subroutine homogenization_thermal_setField(T,dot_T, ce) +module subroutine homogenization_thermal_setField(T,dot_T) - integer, intent(in) :: ce - real(pREAL), intent(in) :: T, dot_T + real(pREAL), dimension(:), intent(in) :: T, dot_T + + integer :: ho, en, ce - current(material_ID_homogenization(ce))%T(material_entry_homogenization(ce)) = T - current(material_ID_homogenization(ce))%dot_T(material_entry_homogenization(ce)) = dot_T - call thermal_partition(ce) + do ce=max(lbound(T,1),lbound(dot_T,1)), min(ubound(T,1),ubound(dot_T,1)) + ho = material_ID_homogenization(ce) + en = material_entry_homogenization(ce) + current(ho)%T(en) = T(ce) + current(ho)%dot_T(en) = dot_T(ce) + call thermal_partition(ce) + end do end subroutine homogenization_thermal_setField diff --git a/src/phase.f90 b/src/phase.f90 index 573e71069..ab365dbd6 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -49,6 +49,29 @@ module phase type(tState), dimension(:), allocatable :: p !< tState for each active source mechanism in a phase end type + enum, bind(c); enumerator :: & + UNDEFINED, & + MECHANICAL_PLASTICITY_NONE, & + MECHANICAL_PLASTICITY_ISOTROPIC, & + MECHANICAL_PLASTICITY_PHENOPOWERLAW, & + MECHANICAL_PLASTICITY_KINEHARDENING, & + MECHANICAL_PLASTICITY_DISLOTWIN, & + MECHANICAL_PLASTICITY_DISLOTUNGSTEN, & + MECHANICAL_PLASTICITY_NONLOCAL, & + MECHANICAL_EIGEN_THERMALEXPANSION, & + DAMAGE_ISOBRITTLE, & + DAMAGE_ANISOBRITTLE, & + THERMAL_SOURCE_DISSIPATION, & + THERMAL_SOURCE_EXTERNALHEAT + end enum + + + integer(kind(UNDEFINED)), dimension(:), allocatable :: & + mechanical_plasticity_type, & !< plasticity of each phase + damage_type !< damage type of each phase + integer(kind(UNDEFINED)), dimension(:,:), allocatable :: & + thermal_source_type, & + mechanical_eigen_kinematics_type character(len=2), allocatable, dimension(:) :: phase_lattice real(pREAL), allocatable, dimension(:) :: phase_cOverA diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 729309b82..2c605c963 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -9,23 +9,14 @@ submodule(phase) damage l_c = 0.0_pREAL !< characteristic length end type tDamageParameters - enum, bind(c); enumerator :: & - DAMAGE_UNDEFINED_ID, & - DAMAGE_ISOBRITTLE_ID, & - DAMAGE_ANISOBRITTLE_ID - end enum - integer :: phase_damage_maxSizeDotState - - type :: tDataContainer + type :: tFieldQuantities real(pREAL), dimension(:), allocatable :: phi - end type tDataContainer + end type tFieldQuantities - integer(kind(DAMAGE_UNDEFINED_ID)), dimension(:), allocatable :: & - phase_damage !< active sources mechanisms of each phase - type(tDataContainer), dimension(:), allocatable :: current + type(tFieldQuantities), dimension(:), allocatable :: current type(tDamageParameters), dimension(:), allocatable :: param @@ -114,11 +105,11 @@ module subroutine damage_init() end do - allocate(phase_damage(phases%length), source = DAMAGE_UNDEFINED_ID) + allocate(damage_type(phases%length), source = UNDEFINED) if (damage_active) then - where(isobrittle_init() ) phase_damage = DAMAGE_ISOBRITTLE_ID - where(anisobrittle_init()) phase_damage = DAMAGE_ANISOBRITTLE_ID + where(isobrittle_init() ) damage_type = DAMAGE_ISOBRITTLE + where(anisobrittle_init()) damage_type = DAMAGE_ANISOBRITTLE end if phase_damage_maxSizeDotState = maxval(damageState%sizeDotState) @@ -159,8 +150,8 @@ module function phase_damage_C66(C66,ph,en) result(C66_degraded) real(pREAL), dimension(6,6) :: C66_degraded - damageType: select case (phase_damage(ph)) - case (DAMAGE_ISOBRITTLE_ID) damageType + damageType: select case (damage_type(ph)) + case (DAMAGE_ISOBRITTLE) damageType C66_degraded = C66 * damage_phi(ph,en)**2 case default damageType C66_degraded = C66 @@ -204,13 +195,14 @@ module function phase_f_phi(phi,co,ce) result(f) ph, & en + ph = material_ID_phase(co,ce) en = material_entry_phase(co,ce) - select case(phase_damage(ph)) - case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ANISOBRITTLE_ID) + select case(damage_type(ph)) + case(DAMAGE_ISOBRITTLE,DAMAGE_ANISOBRITTLE) f = 1.0_pREAL & - - 2.0_pREAL * phi*damageState(ph)%state(1,en) + - 2.0_pREAL * phi*damageState(ph)%state(1,en) ! ToDo: MD: seems to be phi**2 case default f = 0.0_pREAL end select @@ -318,8 +310,8 @@ module subroutine damage_restartWrite(groupHandle,ph) integer, intent(in) :: ph - select case(phase_damage(ph)) - case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ANISOBRITTLE_ID) + select case(damage_type(ph)) + case(DAMAGE_ISOBRITTLE,DAMAGE_ANISOBRITTLE) call HDF5_write(damageState(ph)%state,groupHandle,'omega_damage') end select @@ -332,8 +324,8 @@ module subroutine damage_restartRead(groupHandle,ph) integer, intent(in) :: ph - select case(phase_damage(ph)) - case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ANISOBRITTLE_ID) + select case(damage_type(ph)) + case(DAMAGE_ISOBRITTLE,DAMAGE_ANISOBRITTLE) call HDF5_read(damageState(ph)%state0,groupHandle,'omega_damage') end select @@ -350,15 +342,15 @@ module subroutine damage_result(group,ph) integer, intent(in) :: ph - if (phase_damage(ph) /= DAMAGE_UNDEFINED_ID) & + if (damage_type(ph) /= UNDEFINED) & call result_closeGroup(result_addGroup(group//'damage')) - sourceType: select case (phase_damage(ph)) + sourceType: select case (damage_type(ph)) - case (DAMAGE_ISOBRITTLE_ID) sourceType + case (DAMAGE_ISOBRITTLE) sourceType call isobrittle_result(ph,group//'damage/') - case (DAMAGE_ANISOBRITTLE_ID) sourceType + case (DAMAGE_ANISOBRITTLE) sourceType call anisobrittle_result(ph,group//'damage/') end select sourceType @@ -381,9 +373,9 @@ function phase_damage_collectDotState(ph,en) result(broken) if (damageState(ph)%sizeState > 0) then - sourceType: select case (phase_damage(ph)) + sourceType: select case (damage_type(ph)) - case (DAMAGE_ANISOBRITTLE_ID) sourceType + case (DAMAGE_ANISOBRITTLE) sourceType call anisobrittle_dotState(mechanical_S(ph,en), ph,en) ! ToDo: use M_d end select sourceType @@ -446,9 +438,9 @@ function phase_damage_deltaState(Fe, ph, en) result(broken) if (damageState(ph)%sizeState == 0) return - sourceType: select case (phase_damage(ph)) + sourceType: select case (damage_type(ph)) - case (DAMAGE_ISOBRITTLE_ID) sourceType + case (DAMAGE_ISOBRITTLE) sourceType call isobrittle_deltaState(phase_homogenizedC66(ph,en), Fe, ph,en) broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,en))) if (.not. broken) then diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 0f931517a..2ec8956aa 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -3,21 +3,6 @@ !---------------------------------------------------------------------------------------------------- submodule(phase) mechanical - - enum, bind(c); enumerator :: & - PLASTIC_UNDEFINED_ID, & - PLASTIC_NONE_ID, & - PLASTIC_ISOTROPIC_ID, & - PLASTIC_PHENOPOWERLAW_ID, & - PLASTIC_KINEHARDENING_ID, & - PLASTIC_DISLOTWIN_ID, & - PLASTIC_DISLOTUNGSTEN_ID, & - PLASTIC_NONLOCAL_ID, & - EIGEN_UNDEFINED_ID, & - EIGEN_CLEAVAGE_OPENING_ID, & - EIGEN_THERMAL_EXPANSION_ID - end enum - type(tTensorContainer), dimension(:), allocatable :: & ! current value phase_mechanical_Fe, & @@ -37,9 +22,6 @@ submodule(phase) mechanical phase_mechanical_S0 - integer(kind(PLASTIC_undefined_ID)), dimension(:), allocatable :: & - phase_plasticity !< plasticity of each phase - interface module subroutine eigen_init(phases) @@ -283,7 +265,7 @@ module subroutine mechanical_init(phases) call elastic_init(phases) allocate(plasticState(phases%length)) - allocate(phase_plasticity(phases%length),source = PLASTIC_UNDEFINED_ID) + allocate(mechanical_plasticity_type(phases%length),source = UNDEFINED) call plastic_init() do ph = 1,phases%length plasticState(ph)%state0 = plasticState(ph)%state @@ -327,24 +309,24 @@ module subroutine mechanical_result(group,ph) call results(group,ph) - select case(phase_plasticity(ph)) + select case(mechanical_plasticity_type(ph)) - case(PLASTIC_ISOTROPIC_ID) + case(MECHANICAL_PLASTICITY_ISOTROPIC) call plastic_isotropic_result(ph,group//'mechanical/') - case(PLASTIC_PHENOPOWERLAW_ID) + case(MECHANICAL_PLASTICITY_PHENOPOWERLAW) call plastic_phenopowerlaw_result(ph,group//'mechanical/') - case(PLASTIC_KINEHARDENING_ID) + case(MECHANICAL_PLASTICITY_KINEHARDENING) call plastic_kinehardening_result(ph,group//'mechanical/') - case(PLASTIC_DISLOTWIN_ID) + case(MECHANICAL_PLASTICITY_DISLOTWIN) call plastic_dislotwin_result(ph,group//'mechanical/') - case(PLASTIC_DISLOTUNGSTEN_ID) + case(MECHANICAL_PLASTICITY_DISLOTUNGSTEN) call plastic_dislotungsten_result(ph,group//'mechanical/') - case(PLASTIC_NONLOCAL_ID) + case(MECHANICAL_PLASTICITY_NONLOCAL) call plastic_nonlocal_result(ph,group//'mechanical/') end select diff --git a/src/phase_mechanical_eigen.f90 b/src/phase_mechanical_eigen.f90 index ec1bcfbbc..5bdb57f57 100644 --- a/src/phase_mechanical_eigen.f90 +++ b/src/phase_mechanical_eigen.f90 @@ -3,15 +3,7 @@ submodule(phase:mechanical) eigen integer, dimension(:), allocatable :: & Nmodels - integer(kind(EIGEN_UNDEFINED_ID)), dimension(:,:), allocatable :: & - model - integer(kind(EIGEN_UNDEFINED_ID)), dimension(:), allocatable :: & - model_damage - interface - module function damage_anisobrittle_init() result(myKinematics) - logical, dimension(:), allocatable :: myKinematics - end function damage_anisobrittle_init module function thermalexpansion_init(kinematics_length) result(myKinematics) integer, intent(in) :: kinematics_length @@ -60,17 +52,12 @@ module subroutine eigen_init(phases) Nmodels(ph) = kinematics%length end do - allocate(model(maxval(Nmodels),phases%length), source = EIGEN_undefined_ID) + allocate(mechanical_eigen_kinematics_type(maxval(Nmodels),phases%length), source = UNDEFINED) if (maxval(Nmodels) /= 0) then - where(thermalexpansion_init(maxval(Nmodels))) model = EIGEN_thermal_expansion_ID + where(thermalexpansion_init(maxval(Nmodels))) mechanical_eigen_kinematics_type = MECHANICAL_EIGEN_THERMALEXPANSION end if - allocate(model_damage(phases%length), source = EIGEN_UNDEFINED_ID) - - where(damage_anisobrittle_init()) model_damage = EIGEN_cleavage_opening_ID - - end subroutine eigen_init @@ -173,17 +160,9 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & dLi_dFi = 0.0_pREAL - plasticType: select case (phase_plasticity(ph)) - case (PLASTIC_isotropic_ID) plasticType - call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,ph,en) - Li = Li + my_Li - dLi_dS = dLi_dS + my_dLi_dS - active = .true. - end select plasticType - KinematicsLoop: do k = 1, Nmodels(ph) - kinematicsType: select case (model(k,ph)) - case (EIGEN_thermal_expansion_ID) kinematicsType + kinematicsType: select case (mechanical_eigen_kinematics_type(k,ph)) + case (MECHANICAL_EIGEN_THERMALEXPANSION) kinematicsType call thermalexpansion_LiAndItsTangent(my_Li, my_dLi_dS, ph,en) Li = Li + my_Li dLi_dS = dLi_dS + my_dLi_dS @@ -191,13 +170,21 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & end select kinematicsType end do KinematicsLoop - select case (model_damage(ph)) - case (EIGEN_cleavage_opening_ID) + plasticType: select case (mechanical_plasticity_type(ph)) + case (MECHANICAL_PLASTICITY_ISOTROPIC) plasticType + call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,ph,en) + Li = Li + my_Li + dLi_dS = dLi_dS + my_dLi_dS + active = .true. + end select plasticType + + damageType: select case (damage_type(ph)) + case (DAMAGE_ANISOBRITTLE) call damage_anisobrittle_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, en) Li = Li + my_Li dLi_dS = dLi_dS + my_dLi_dS active = .true. - end select + end select damageType if (.not. active) return diff --git a/src/phase_mechanical_eigen_cleavageopening.f90 b/src/phase_mechanical_eigen_cleavageopening.f90 deleted file mode 100644 index 780ed22b2..000000000 --- a/src/phase_mechanical_eigen_cleavageopening.f90 +++ /dev/null @@ -1,30 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH -!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine incorporating kinematics resulting from opening of cleavage planes -!> @details to be done -!-------------------------------------------------------------------------------------------------- -submodule(phase:eigen) cleavageopening - -contains - - -!-------------------------------------------------------------------------------------------------- -!> @brief module initialization -!> @details reads in material parameters, allocates arrays, and does sanity checks -!-------------------------------------------------------------------------------------------------- -module function damage_anisobrittle_init() result(myKinematics) - - logical, dimension(:), allocatable :: myKinematics - - - myKinematics = kinematics_active2('anisobrittle') - if (count(myKinematics) == 0) return - - print'(/,1x,a)', '<<<+- phase:mechanical:eigen:cleavageopening init -+>>>' - print'(/,a,i2)', ' # phases: ',count(myKinematics); flush(IO_STDOUT) - -end function damage_anisobrittle_init - - -end submodule cleavageopening diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index a2b5becff..89b5027ac 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -199,8 +199,8 @@ module function phase_homogenizedC66(ph,en) result(C) integer, intent(in) :: ph, en - plasticType: select case (phase_plasticity(ph)) - case (PLASTIC_DISLOTWIN_ID) plasticType + plasticType: select case (mechanical_plasticity_type(ph)) + case (MECHANICAL_PLASTICITY_DISLOTWIN) plasticType C = plastic_dislotwin_homogenizedC(ph,en) case default plasticType C = elastic_C66(ph,en) diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index 0c1959660..5d1a462e7 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -211,17 +211,17 @@ contains module subroutine plastic_init - print'(/,1x,a)', '<<<+- phase:mechanical:plastic init -+>>>' + print'(/,1x,a)', '<<<+- phase:mechanical:plasticity init -+>>>' - where(plastic_none_init()) phase_plasticity = PLASTIC_NONE_ID - where(plastic_isotropic_init()) phase_plasticity = PLASTIC_ISOTROPIC_ID - where(plastic_phenopowerlaw_init()) phase_plasticity = PLASTIC_PHENOPOWERLAW_ID - where(plastic_kinehardening_init()) phase_plasticity = PLASTIC_KINEHARDENING_ID - where(plastic_dislotwin_init()) phase_plasticity = PLASTIC_DISLOTWIN_ID - where(plastic_dislotungsten_init()) phase_plasticity = PLASTIC_DISLOTUNGSTEN_ID - where(plastic_nonlocal_init()) phase_plasticity = PLASTIC_NONLOCAL_ID + where(plastic_none_init()) mechanical_plasticity_type = MECHANICAL_PLASTICITY_NONE + where(plastic_isotropic_init()) mechanical_plasticity_type = MECHANICAL_PLASTICITY_ISOTROPIC + where(plastic_phenopowerlaw_init()) mechanical_plasticity_type = MECHANICAL_PLASTICITY_PHENOPOWERLAW + where(plastic_kinehardening_init()) mechanical_plasticity_type = MECHANICAL_PLASTICITY_KINEHARDENING + where(plastic_dislotwin_init()) mechanical_plasticity_type = MECHANICAL_PLASTICITY_DISLOTWIN + where(plastic_dislotungsten_init()) mechanical_plasticity_type = MECHANICAL_PLASTICITY_DISLOTUNGSTEN + where(plastic_nonlocal_init()) mechanical_plasticity_type = MECHANICAL_PLASTICITY_NONLOCAL - if (any(phase_plasticity == PLASTIC_undefined_ID)) call IO_error(201) + if (any(mechanical_plasticity_type == UNDEFINED)) call IO_error(201) end subroutine plastic_init @@ -251,7 +251,7 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & i, j - if (phase_plasticity(ph) == PLASTIC_NONE_ID) then + if (mechanical_plasticity_type(ph) == MECHANICAL_PLASTICITY_NONE) then Lp = 0.0_pREAL dLp_dFi = 0.0_pREAL dLp_dS = 0.0_pREAL @@ -259,24 +259,24 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & Mp = matmul(matmul(transpose(Fi),Fi),S) - plasticType: select case (phase_plasticity(ph)) + plasticType: select case (mechanical_plasticity_type(ph)) - case (PLASTIC_ISOTROPIC_ID) plasticType + case (MECHANICAL_PLASTICITY_ISOTROPIC) plasticType call isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) - case (PLASTIC_PHENOPOWERLAW_ID) plasticType + case (MECHANICAL_PLASTICITY_PHENOPOWERLAW) plasticType call phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) - case (PLASTIC_KINEHARDENING_ID) plasticType + case (MECHANICAL_PLASTICITY_KINEHARDENING) plasticType call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) - case (PLASTIC_NONLOCAL_ID) plasticType + case (MECHANICAL_PLASTICITY_NONLOCAL) plasticType call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) - case (PLASTIC_DISLOTWIN_ID) plasticType + case (MECHANICAL_PLASTICITY_DISLOTWIN) plasticType call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) - case (PLASTIC_DISLOTUNGSTEN_ID) plasticType + case (MECHANICAL_PLASTICITY_DISLOTUNGSTEN) plasticType call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) end select plasticType @@ -308,28 +308,28 @@ module function plastic_dotState(subdt,ph,en) result(dotState) dotState - if (phase_plasticity(ph) /= PLASTIC_NONE_ID) then + if (mechanical_plasticity_type(ph) /= MECHANICAL_PLASTICITY_NONE) then Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,en)),& phase_mechanical_Fi(ph)%data(1:3,1:3,en)),phase_mechanical_S(ph)%data(1:3,1:3,en)) - plasticType: select case (phase_plasticity(ph)) + plasticType: select case (mechanical_plasticity_type(ph)) - case (PLASTIC_ISOTROPIC_ID) plasticType + case (MECHANICAL_PLASTICITY_ISOTROPIC) plasticType dotState = isotropic_dotState(Mp,ph,en) - case (PLASTIC_PHENOPOWERLAW_ID) plasticType + case (MECHANICAL_PLASTICITY_PHENOPOWERLAW) plasticType dotState = phenopowerlaw_dotState(Mp,ph,en) - case (PLASTIC_KINEHARDENING_ID) plasticType + case (MECHANICAL_PLASTICITY_KINEHARDENING) plasticType dotState = plastic_kinehardening_dotState(Mp,ph,en) - case (PLASTIC_DISLOTWIN_ID) plasticType + case (MECHANICAL_PLASTICITY_DISLOTWIN) plasticType dotState = dislotwin_dotState(Mp,ph,en) - case (PLASTIC_DISLOTUNGSTEN_ID) plasticType + case (MECHANICAL_PLASTICITY_DISLOTUNGSTEN) plasticType dotState = dislotungsten_dotState(Mp,ph,en) - case (PLASTIC_NONLOCAL_ID) plasticType + case (MECHANICAL_PLASTICITY_NONLOCAL) plasticType call nonlocal_dotState(Mp,subdt,ph,en) dotState = plasticState(ph)%dotState(:,en) @@ -349,15 +349,15 @@ module subroutine plastic_dependentState(ph,en) en - plasticType: select case (phase_plasticity(ph)) + plasticType: select case (mechanical_plasticity_type(ph)) - case (PLASTIC_DISLOTWIN_ID) plasticType + case (MECHANICAL_PLASTICITY_DISLOTWIN) plasticType call dislotwin_dependentState(ph,en) - case (PLASTIC_DISLOTUNGSTEN_ID) plasticType + case (MECHANICAL_PLASTICITY_DISLOTUNGSTEN) plasticType call dislotungsten_dependentState(ph,en) - case (PLASTIC_NONLOCAL_ID) plasticType + case (MECHANICAL_PLASTICITY_NONLOCAL) plasticType call nonlocal_dependentState(ph,en) end select plasticType @@ -384,19 +384,19 @@ module function plastic_deltaState(ph, en) result(broken) broken = .false. - select case (phase_plasticity(ph)) - case (PLASTIC_NONLOCAL_ID,PLASTIC_KINEHARDENING_ID) + select case (mechanical_plasticity_type(ph)) + case (MECHANICAL_PLASTICITY_NONLOCAL,MECHANICAL_PLASTICITY_KINEHARDENING) Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,en)),& phase_mechanical_Fi(ph)%data(1:3,1:3,en)),& phase_mechanical_S(ph)%data(1:3,1:3,en)) - plasticType: select case (phase_plasticity(ph)) + plasticType: select case (mechanical_plasticity_type(ph)) - case (PLASTIC_KINEHARDENING_ID) plasticType + case (MECHANICAL_PLASTICITY_KINEHARDENING) plasticType call plastic_kinehardening_deltaState(Mp,ph,en) - case (PLASTIC_NONLOCAL_ID) plasticType + case (MECHANICAL_PLASTICITY_NONLOCAL) plasticType call plastic_nonlocal_deltaState(Mp,ph,en) end select plasticType diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index d436f1095..3e34256f1 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -1252,7 +1252,7 @@ function rhoDotFlux(timestep,ph,en) !* The entering flux from my neighbor will be distributed on my slip systems according to the !* compatibility if (neighbor_n > 0) then - if (phase_plasticity(np) == PLASTIC_NONLOCAL_ID .and. & + if (mechanical_plasticity_type(np) == MECHANICAL_PLASTICITY_NONLOCAL .and. & any(dependentState(ph)%compatibility(:,:,:,n,en) > 0.0_pREAL)) then forall (s = 1:ns, t = 1:4) @@ -1298,7 +1298,7 @@ function rhoDotFlux(timestep,ph,en) !* In case of reduced transmissivity, part of the leaving flux is stored as dead dislocation density. !* That means for an interface of zero transmissivity the leaving flux is fully converted to dead dislocations. if (opposite_n > 0) then - if (phase_plasticity(np) == PLASTIC_NONLOCAL_ID) then + if (mechanical_plasticity_type(np) == MECHANICAL_PLASTICITY_NONLOCAL) then normal_me2neighbor_defConf = math_det33(Favg) & * matmul(math_inv33(transpose(Favg)),geom(ph)%IPareaNormal(1:3,n,en)) ! normal of the interface in (average) deformed configuration (pointing en => neighbor) diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90 index 5ef4dc7fd..e21f173ad 100644 --- a/src/phase_thermal.f90 +++ b/src/phase_thermal.f90 @@ -15,19 +15,11 @@ submodule(phase) thermal type(tSourceState), allocatable, dimension(:) :: & thermalState - enum, bind(c); enumerator :: & - THERMAL_UNDEFINED_ID ,& - THERMAL_DISSIPATION_ID, & - THERMAL_EXTERNALHEAT_ID - end enum - - type :: tDataContainer ! ?? not very telling name. Better: "fieldQuantities" ?? + type :: tFieldQuantities real(pREAL), dimension(:), allocatable :: T, dot_T - end type tDataContainer - integer(kind(THERMAL_UNDEFINED_ID)), dimension(:,:), allocatable :: & - thermal_source + end type tFieldQuantities - type(tDataContainer), dimension(:), allocatable :: current ! ?? not very telling name. Better: "field" ?? MD: current(ho)%T(en) reads quite good + type(tFieldQuantities), dimension(:), allocatable :: current type(tThermalParameters), dimension(:), allocatable :: param @@ -36,36 +28,36 @@ submodule(phase) thermal interface - module function dissipation_init(source_length) result(mySources) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources - end function dissipation_init + module function source_dissipation_init(maxNsources) result(isMySource) + integer, intent(in) :: maxNsources + logical, dimension(:,:), allocatable :: isMySource + end function source_dissipation_init - module function externalheat_init(source_length) result(mySources) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources - end function externalheat_init + module function source_externalheat_init(maxNsources) result(isMySource) + integer, intent(in) :: maxNsources + logical, dimension(:,:), allocatable :: isMySource + end function source_externalheat_init - module subroutine externalheat_dotState(ph, en) + module subroutine source_externalheat_dotState(ph, en) integer, intent(in) :: & ph, & en - end subroutine externalheat_dotState + end subroutine source_externalheat_dotState - module function dissipation_f_T(ph,en) result(f_T) + module function source_dissipation_f_T(ph,en) result(f_T) integer, intent(in) :: & ph, & en real(pREAL) :: f_T - end function dissipation_f_T + end function source_dissipation_f_T - module function externalheat_f_T(ph,en) result(f_T) + module function source_externalheat_f_T(ph,en) result(f_T) integer, intent(in) :: & ph, & en real(pREAL) :: f_T - end function externalheat_f_T + end function source_externalheat_f_T end interface @@ -129,11 +121,11 @@ module subroutine thermal_init(phases) end do - allocate(thermal_source(maxval(thermal_Nsources),phases%length), source = THERMAL_UNDEFINED_ID) + allocate(thermal_source_type(maxval(thermal_Nsources),phases%length), source = UNDEFINED) if (maxval(thermal_Nsources) /= 0) then - where(dissipation_init (maxval(thermal_Nsources))) thermal_source = THERMAL_DISSIPATION_ID - where(externalheat_init(maxval(thermal_Nsources))) thermal_source = THERMAL_EXTERNALHEAT_ID + where(source_dissipation_init (maxval(thermal_Nsources))) thermal_source_type = THERMAL_SOURCE_DISSIPATION + where(source_externalheat_init(maxval(thermal_Nsources))) thermal_source_type = THERMAL_SOURCE_EXTERNALHEAT end if thermal_source_maxSizeDotState = 0 @@ -151,7 +143,7 @@ end subroutine thermal_init !---------------------------------------------------------------------------------------------- -!< @brief Calculate thermal source. +!< @brief Calculate thermal source (forcing term). !---------------------------------------------------------------------------------------------- module function phase_f_T(ph,en) result(f) @@ -165,13 +157,13 @@ module function phase_f_T(ph,en) result(f) f = 0.0_pREAL do so = 1, thermal_Nsources(ph) - select case(thermal_source(so,ph)) + select case(thermal_source_type(so,ph)) - case (THERMAL_DISSIPATION_ID) - f = f + dissipation_f_T(ph,en) + case (THERMAL_SOURCE_DISSIPATION) + f = f + source_dissipation_f_T(ph,en) - case (THERMAL_EXTERNALHEAT_ID) - f = f + externalheat_f_T(ph,en) + case (THERMAL_SOURCE_EXTERNALHEAT) + f = f + source_externalheat_f_T(ph,en) end select @@ -183,22 +175,22 @@ end function phase_f_T !-------------------------------------------------------------------------------------------------- !> @brief tbd. !-------------------------------------------------------------------------------------------------- -function phase_thermal_collectDotState(ph,en) result(broken) +function phase_thermal_collectDotState(ph,en) result(ok) integer, intent(in) :: ph, en - logical :: broken + logical :: ok integer :: i - broken = .false. + ok = .true. SourceLoop: do i = 1, thermal_Nsources(ph) - if (thermal_source(i,ph) == THERMAL_EXTERNALHEAT_ID) & - call externalheat_dotState(ph,en) + if (thermal_source_type(i,ph) == THERMAL_SOURCE_EXTERNALHEAT) & + call source_externalheat_dotState(ph,en) - broken = broken .or. any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,en))) + ok = ok .and. .not. any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,en))) end do SourceLoop @@ -241,34 +233,35 @@ module function phase_thermal_constitutive(Delta_t,ph,en) result(converged_) logical :: converged_ - converged_ = .not. integrateThermalState(Delta_t,ph,en) + converged_ = integrateThermalState(Delta_t,ph,en) end function phase_thermal_constitutive !-------------------------------------------------------------------------------------------------- -!> @brief integrate state with 1st order explicit Euler method +!> @brief Integrate state with 1st order explicit Euler method. !-------------------------------------------------------------------------------------------------- -function integrateThermalState(Delta_t, ph,en) result(broken) +function integrateThermalState(Delta_t, ph,en) result(converged) real(pREAL), intent(in) :: Delta_t integer, intent(in) :: ph, en - logical :: & - broken + logical :: converged integer :: & so, & sizeDotState - broken = phase_thermal_collectDotState(ph,en) - if (broken) return + converged = phase_thermal_collectDotState(ph,en) + if (converged) then - do so = 1, thermal_Nsources(ph) - sizeDotState = thermalState(ph)%p(so)%sizeDotState - thermalState(ph)%p(so)%state(1:sizeDotState,en) = thermalState(ph)%p(so)%state0(1:sizeDotState,en) & - + thermalState(ph)%p(so)%dotState(1:sizeDotState,en) * Delta_t - end do + do so = 1, thermal_Nsources(ph) + sizeDotState = thermalState(ph)%p(so)%sizeDotState + thermalState(ph)%p(so)%state(1:sizeDotState,en) = thermalState(ph)%p(so)%state0(1:sizeDotState,en) & + + thermalState(ph)%p(so)%dotState(1:sizeDotState,en) * Delta_t + end do + + end if end function integrateThermalState @@ -318,7 +311,7 @@ end subroutine thermal_forward !---------------------------------------------------------------------------------------------- -!< @brief Get temperature (for use by non-thermal physics) +!< @brief Get temperature (for use by non-thermal physics). !---------------------------------------------------------------------------------------------- pure module function thermal_T(ph,en) result(T) @@ -332,7 +325,7 @@ end function thermal_T !---------------------------------------------------------------------------------------------- -!< @brief Get rate of temperature (for use by non-thermal physics) +!< @brief Get rate of temperature (for use by non-thermal physics). !---------------------------------------------------------------------------------------------- module function thermal_dot_T(ph,en) result(dot_T) diff --git a/src/phase_thermal_dissipation.f90 b/src/phase_thermal_source_dissipation.f90 similarity index 76% rename from src/phase_thermal_dissipation.f90 rename to src/phase_thermal_source_dissipation.f90 index 573921670..760259286 100644 --- a/src/phase_thermal_dissipation.f90 +++ b/src/phase_thermal_source_dissipation.f90 @@ -5,7 +5,7 @@ !> @brief material subroutine for thermal source due to plastic dissipation !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule(phase:thermal) dissipation +submodule(phase:thermal) source_dissipation type :: tParameters !< container type for internal constitutive parameters real(pREAL) :: & @@ -22,10 +22,10 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function dissipation_init(source_length) result(mySources) +module function source_dissipation_init(maxNsources) result(isMySource) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources + integer, intent(in) :: maxNsources + logical, dimension(:,:), allocatable :: isMySource type(tDict), pointer :: & phases, & @@ -35,26 +35,29 @@ module function dissipation_init(source_length) result(mySources) class(tList), pointer :: & sources character(len=:), allocatable :: refs - integer :: so,Nmembers,ph + integer :: ph,Nmembers,so,Nsources - mySources = thermal_active('dissipation',source_length) - if (count(mySources) == 0) return + isMySource = thermal_active('dissipation',maxNsources) + if (count(isMySource) == 0) return - print'(/,1x,a)', '<<<+- phase:thermal:dissipation init -+>>>' - print'(/,a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- phase:thermal:source_dissipation init -+>>>' + print'(/,a,i2)', ' # phases: ',count(isMySource); flush(IO_STDOUT) phases => config_material%get_dict('phase') allocate(param(phases%length)) do ph = 1, phases%length + Nsources = count(isMySource(:,ph)) + if (Nsources == 0) cycle + if (Nsources > 1) call IO_error(600,ext_msg='dissipation') + Nmembers = count(material_ID_phase == ph) phase => phases%get_dict(ph) - if (count(mySources(:,ph)) == 0) cycle !ToDo: error if > 1 thermal => phase%get_dict('thermal') sources => thermal%get_list('source') do so = 1, sources%length - if (mySources(so,ph)) then + if (isMySource(so,ph)) then associate(prm => param(ph)) src => sources%get_dict(so) print'(1x,a,i0,a,i0)', 'phase ',ph,' source ',so @@ -62,22 +65,21 @@ module function dissipation_init(source_length) result(mySources) if (len(refs) > 0) print'(/,1x,a)', refs prm%kappa = src%get_asReal('kappa') - Nmembers = count(material_ID_phase == ph) call phase_allocateState(thermalState(ph)%p(so),Nmembers,0,0,0) - end associate + exit end if end do end do -end function dissipation_init +end function source_dissipation_init !-------------------------------------------------------------------------------------------------- !> @brief Ninstancess dissipation rate !-------------------------------------------------------------------------------------------------- -module function dissipation_f_T(ph,en) result(f_T) +module function source_dissipation_f_T(ph,en) result(f_T) integer, intent(in) :: ph, en real(pREAL) :: & @@ -91,6 +93,6 @@ module function dissipation_f_T(ph,en) result(f_T) f_T = prm%kappa*sum(abs(Mp*mechanical_L_p(ph,en))) end associate -end function dissipation_f_T +end function source_dissipation_f_T -end submodule dissipation +end submodule source_dissipation diff --git a/src/phase_thermal_externalheat.f90 b/src/phase_thermal_source_externalheat.f90 similarity index 63% rename from src/phase_thermal_externalheat.f90 rename to src/phase_thermal_source_externalheat.f90 index cdd037592..1a803549c 100644 --- a/src/phase_thermal_externalheat.f90 +++ b/src/phase_thermal_source_externalheat.f90 @@ -4,14 +4,14 @@ !> @author Philip Eisenlohr, Michigan State University !> @brief material subroutine for variable heat source !-------------------------------------------------------------------------------------------------- -submodule(phase:thermal) externalheat +submodule(phase:thermal) source_externalheat integer, dimension(:), allocatable :: & - source_thermal_externalheat_offset !< which source is my current thermal dissipation mechanism? + source_ID !< index in phase source list corresponding to this source type :: tParameters !< container type for internal constitutive parameters - type(tTable) :: f + type(tTable) :: f !< external heat power as (tabulated) function of time end type tParameters type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances) @@ -24,10 +24,10 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function externalheat_init(source_length) result(mySources) +module function source_externalheat_init(maxNsources) result(isMySource) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources + integer, intent(in) :: maxNsources + logical, dimension(:,:), allocatable :: isMySource type(tDict), pointer :: & phases, & @@ -37,28 +37,31 @@ module function externalheat_init(source_length) result(mySources) type(tList), pointer :: & sources character(len=:), allocatable :: refs - integer :: so,Nmembers,ph + integer :: ph,Nmembers,so,Nsources - mySources = thermal_active('externalheat',source_length) - if (count(mySources) == 0) return + isMySource = thermal_active('externalheat',maxNsources) + if (count(isMySource) == 0) return - print'(/,1x,a)', '<<<+- phase:thermal:externalheat init -+>>>' - print'(/,a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT) + print'(/,1x,a)', '<<<+- phase:thermal:source_externalheat init -+>>>' + print'(/,a,i2)', ' # phases: ',count(isMySource); flush(IO_STDOUT) phases => config_material%get_dict('phase') allocate(param(phases%length)) - allocate(source_thermal_externalheat_offset (phases%length), source=0) + allocate(source_ID(phases%length), source=0) do ph = 1, phases%length + Nsources = count(isMySource(:,ph)) + if (Nsources == 0) cycle + if (Nsources > 1) call IO_error(600,ext_msg='externalheat') + Nmembers = count(material_ID_phase == ph) phase => phases%get_dict(ph) - if (count(mySources(:,ph)) == 0) cycle thermal => phase%get_dict('thermal') sources => thermal%get_list('source') do so = 1, sources%length - if (mySources(so,ph)) then - source_thermal_externalheat_offset(ph) = so + if (isMySource(so,ph)) then + source_ID(ph) = so associate(prm => param(ph)) src => sources%get_dict(so) print'(1x,a,i0,a,i0)', 'phase ',ph,' source ',so @@ -66,41 +69,36 @@ module function externalheat_init(source_length) result(mySources) if (len(refs) > 0) print'(/,1x,a)', refs prm%f = table(src,'t','f') - - Nmembers = count(material_ID_phase == ph) call phase_allocateState(thermalState(ph)%p(so),Nmembers,1,1,0) end associate + exit end if end do end do -end function externalheat_init +end function source_externalheat_init !-------------------------------------------------------------------------------------------------- !> @brief rate of change of state !> @details state only contains current time to linearly interpolate given heat powers !-------------------------------------------------------------------------------------------------- -module subroutine externalheat_dotState(ph, en) +module subroutine source_externalheat_dotState(ph, en) integer, intent(in) :: & ph, & en - integer :: & - so - so = source_thermal_externalheat_offset(ph) + thermalState(ph)%p(source_ID(ph))%dotState(1,en) = 1.0_pREAL ! state is current time - thermalState(ph)%p(so)%dotState(1,en) = 1.0_pREAL ! state is current time - -end subroutine externalheat_dotState +end subroutine source_externalheat_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns local heat generation rate !-------------------------------------------------------------------------------------------------- -module function externalheat_f_T(ph,en) result(f_T) +module function source_externalheat_f_T(ph,en) result(f_T) integer, intent(in) :: & ph, & @@ -108,16 +106,11 @@ module function externalheat_f_T(ph,en) result(f_T) real(pREAL) :: & f_T - integer :: & - so - - - so = source_thermal_externalheat_offset(ph) associate(prm => param(ph)) - f_T = prm%f%at(thermalState(ph)%p(so)%state(1,en)) + f_T = prm%f%at(thermalState(ph)%p(source_ID(ph))%state(1,en)) end associate -end function externalheat_f_T +end function source_externalheat_f_T -end submodule externalheat +end submodule source_externalheat