From 5d88a782067714ad076dc8631d2d34c85da15366 Mon Sep 17 00:00:00 2001 From: Pratheek Shanthraj Date: Tue, 24 Jun 2014 22:59:16 +0000 Subject: [PATCH] added thermal and damage modules as examples of multi physics modules. only works with new state layout and still under testing. damage_none: does nothing damage_gradient: interacts with solver to solve gradient damage problems thermal_none: does nothing thermai_adiabatic: local heating only thermal_conduction: interacts with conduction solver to solve coupled heat transfer problems --- code/constitutive_damage.f90 | 270 ++++++++++++++++++++++++ code/constitutive_thermal.f90 | 266 ++++++++++++++++++++++++ code/damage_gradient.f90 | 377 ++++++++++++++++++++++++++++++++++ code/damage_none.f90 | 106 ++++++++++ code/thermal_adiabatic.f90 | 326 +++++++++++++++++++++++++++++ code/thermal_conduction.f90 | 325 +++++++++++++++++++++++++++++ code/thermal_none.f90 | 106 ++++++++++ 7 files changed, 1776 insertions(+) create mode 100644 code/constitutive_damage.f90 create mode 100644 code/constitutive_thermal.f90 create mode 100644 code/damage_gradient.f90 create mode 100644 code/damage_none.f90 create mode 100644 code/thermal_adiabatic.f90 create mode 100644 code/thermal_conduction.f90 create mode 100644 code/thermal_none.f90 diff --git a/code/constitutive_damage.f90 b/code/constitutive_damage.f90 new file mode 100644 index 000000000..d8a114b54 --- /dev/null +++ b/code/constitutive_damage.f90 @@ -0,0 +1,270 @@ +!-------------------------------------------------------------------------------------------------- +! $Id: constitutive_damage.f90 3205 2014-06-17 06:54:49Z MPIE\m.diehl $ +!-------------------------------------------------------------------------------------------------- +!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH +!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH +!> @brief damage internal microstructure state +!-------------------------------------------------------------------------------------------------- +module constitutive_damage + use prec, only: & + pInt, & + pReal + + implicit none + private + integer(pInt), public, dimension(:,:,:), allocatable :: & + constitutive_damage_sizePostResults !< size of postResults array per grain + integer(pInt), public, protected :: & + constitutive_damage_maxSizePostResults, & + constitutive_damage_maxSizeDotState + public :: & + constitutive_damage_init, & + constitutive_damage_microstructure, & + constitutive_damage_collectDotState, & + constitutive_damage_collectDeltaState, & + constitutive_damage_postResults + +contains + + +!-------------------------------------------------------------------------------------------------- +!> @brief allocates arrays pointing to array of the various constitutive modules +!-------------------------------------------------------------------------------------------------- +subroutine constitutive_damage_init + + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) + use IO, only: & + IO_open_file, & + IO_open_jobFile_stat, & + IO_write_jobFile, & + IO_timeStamp + use mesh, only: & + mesh_maxNips, & + mesh_NcpElems, & + mesh_element, & + FE_Nips, & + FE_geomtype + use material, only: & + material_phase, & + material_Nphase, & + material_localFileExt, & + material_configFile, & + phase_name, & + phase_damage, & + phase_damageInstance, & + phase_Noutput, & + homogenization_Ngrains, & + homogenization_maxNgrains, & + damageState, & + DAMAGE_none_ID, & + DAMAGE_NONE_label, & + DAMAGE_gradient_ID, & + DAMAGE_GRADIENT_label +use damage_none +use damage_gradient + + implicit none + integer(pInt), parameter :: FILEUNIT = 200_pInt + integer(pInt) :: & + g, & !< grain number + i, & !< integration point number + e, & !< element number + cMax, & !< maximum number of grains + iMax, & !< maximum number of integration points + eMax, & !< maximum number of elements + phase, & + s, & + p, & + instance,& + myNgrains + + integer(pInt), dimension(:,:), pointer :: thisSize + logical :: knownDamage + character(len=64), dimension(:,:), pointer :: thisOutput + character(len=32) :: outputName !< name of output, intermediate fix until HDF5 output is ready + +!-------------------------------------------------------------------------------------------------- +! parse plasticities from config file + if (.not. IO_open_jobFile_stat(FILEUNIT,material_localFileExt)) & ! no local material configuration present... + call IO_open_file(FILEUNIT,material_configFile) ! ... open material.config file + if (any(phase_damage == DAMAGE_none_ID)) call damage_none_init(FILEUNIT) + if (any(phase_damage == DAMAGE_gradient_ID)) call damage_gradient_init(FILEUNIT) + close(FILEUNIT) + + write(6,'(/,a)') ' <<<+- constitutive_damage init -+>>>' + write(6,'(a)') ' $Id: constitutive_damage.f90 3205 2014-06-17 06:54:49Z MPIE\m.diehl $' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() +#include "compilation_info.f90" + +!-------------------------------------------------------------------------------------------------- +! write description file for constitutive phase output + call IO_write_jobFile(FILEUNIT,'outputDamage') + do phase = 1_pInt,material_Nphase + instance = phase_damageInstance(phase) ! which instance of a plasticity is present phase + knownDamage = .true. + select case(phase_damage(phase)) ! split per constititution + case (DAMAGE_none_ID) + outputName = DAMAGE_NONE_label + thisOutput => null() + thisSize => null() + case (DAMAGE_gradient_ID) + outputName = DAMAGE_GRADIENT_label + thisOutput => damage_gradient_output + thisSize => damage_gradient_sizePostResult + case default + knownDamage = .false. + end select + write(FILEUNIT,'(/,a,/)') '['//trim(phase_name(phase))//']' + if (knownDamage) then + write(FILEUNIT,'(a)') '(damage)'//char(9)//trim(outputName) + if (phase_damage(phase) /= DAMAGE_none_ID) then + do e = 1_pInt,phase_Noutput(phase) + write(FILEUNIT,'(a,i4)') trim(thisOutput(e,instance))//char(9),thisSize(e,instance) + enddo + endif + endif + enddo + close(FILEUNIT) + +!-------------------------------------------------------------------------------------------------- +! allocation of states + cMax = homogenization_maxNgrains + iMax = mesh_maxNips + eMax = mesh_NcpElems + allocate(constitutive_damage_sizePostResults(cMax,iMax,eMax), source=0_pInt) + + ElemLoop:do e = 1_pInt,mesh_NcpElems ! loop over elements + myNgrains = homogenization_Ngrains(mesh_element(3,e)) + IPloop:do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e))) ! loop over IPs + GrainLoop:do g = 1_pInt,myNgrains ! loop over grains + phase = material_phase(g,i,e) + instance = phase_damageInstance(phase) + select case(phase_damage(phase)) + case (DAMAGE_gradient_ID) + constitutive_damage_sizePostResults(g,i,e) = damage_gradient_sizePostResults(instance) + + end select + enddo GrainLoop + enddo IPloop + enddo ElemLoop + + constitutive_damage_maxSizePostResults = maxval(constitutive_damage_sizePostResults) + constitutive_damage_maxSizeDotState = 0_pInt + do p = 1, size(damageState) + constitutive_damage_maxSizeDotState = max(constitutive_damage_maxSizeDotState, damageState(p)%sizeDotState) + enddo +end subroutine constitutive_damage_init + + +!-------------------------------------------------------------------------------------------------- +!> @brief calls microstructure function of the different constitutive models +!-------------------------------------------------------------------------------------------------- +subroutine constitutive_damage_microstructure(Tstar_v, Fe, ipc, ip, el) + use material, only: & + material_phase, & + phase_damage, & + DAMAGE_gradient_ID + use damage_gradient, only: & + damage_gradient_microstructure + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), intent(in), dimension(6) :: & + Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel) + real(pReal), intent(in), dimension(3,3) :: & + Fe + + select case (phase_damage(material_phase(ipc,ip,el))) + case (DAMAGE_gradient_ID) + call damage_gradient_microstructure(Tstar_v, Fe, ipc, ip, el) + + end select + +end subroutine constitutive_damage_microstructure + + +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the rate of change of microstructure +!-------------------------------------------------------------------------------------------------- +subroutine constitutive_damage_collectDotState(Tstar_v, Lp, ipc, ip, el) + use material, only: & + material_phase, & + phase_damage, & + DAMAGE_gradient_ID + use damage_gradient, only: & + damage_gradient_dotState + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), intent(in), dimension(6) :: & + Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel) + real(pReal), intent(in), dimension(3,3) :: & + Lp + + select case (phase_damage(material_phase(ipc,ip,el))) + case (DAMAGE_gradient_ID) + call damage_gradient_dotState(Tstar_v, Lp, ipc, ip, el) + + end select + +end subroutine constitutive_damage_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_damage_collectDeltaState(ipc, ip, el) + use material, only: & + material_phase, & + phase_damage + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + + select case (phase_damage(material_phase(ipc,ip,el))) + + end select + constitutive_damage_collectDeltaState = .true. + +end function constitutive_damage_collectDeltaState + + +!-------------------------------------------------------------------------------------------------- +!> @brief returns array of constitutive results +!-------------------------------------------------------------------------------------------------- +function constitutive_damage_postResults(ipc, ip, el) + use material, only: & + material_phase, & + phase_damage, & + DAMAGE_gradient_ID + use damage_gradient, only: & + damage_gradient_postResults + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), dimension(constitutive_damage_sizePostResults(ipc,ip,el)) :: & + constitutive_damage_postResults + + constitutive_damage_postResults = 0.0_pReal + + select case (phase_damage(material_phase(ipc,ip,el))) + case (DAMAGE_gradient_ID) + constitutive_damage_postResults = damage_gradient_postResults(ipc,ip,el) + end select + +end function constitutive_damage_postResults + + +end module constitutive_damage diff --git a/code/constitutive_thermal.f90 b/code/constitutive_thermal.f90 new file mode 100644 index 000000000..071fa2bc5 --- /dev/null +++ b/code/constitutive_thermal.f90 @@ -0,0 +1,266 @@ +!-------------------------------------------------------------------------------------------------- +! $Id: constitutive_thermal.f90 3205 2014-06-17 06:54:49Z MPIE\m.diehl $ +!-------------------------------------------------------------------------------------------------- +!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH +!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH +!> @brief thermal internal microstructure state +!-------------------------------------------------------------------------------------------------- +module constitutive_thermal + use prec, only: & + pInt, & + pReal + + implicit none + private + integer(pInt), public, dimension(:,:,:), allocatable :: & + constitutive_thermal_sizePostResults !< size of postResults array per grain + integer(pInt), public, protected :: & + constitutive_thermal_maxSizePostResults, & + constitutive_thermal_maxSizeDotState + public :: & + constitutive_thermal_init, & + constitutive_thermal_microstructure, & + constitutive_thermal_collectDotState, & + constitutive_thermal_collectDeltaState, & + constitutive_thermal_postResults + +contains + + +!-------------------------------------------------------------------------------------------------- +!> @brief allocates arrays pointing to array of the various constitutive modules +!-------------------------------------------------------------------------------------------------- +subroutine constitutive_thermal_init + + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) + use IO, only: & + IO_open_file, & + IO_open_jobFile_stat, & + IO_write_jobFile, & + IO_timeStamp + use mesh, only: & + mesh_maxNips, & + mesh_NcpElems, & + mesh_element, & + FE_Nips, & + FE_geomtype + use material, only: & + material_phase, & + material_Nphase, & + material_localFileExt, & + material_configFile, & + phase_name, & + phase_thermal, & + phase_thermalInstance, & + phase_Noutput, & + homogenization_Ngrains, & + homogenization_maxNgrains, & + thermalState, & + THERMAL_none_ID, & + THERMAL_NONE_label, & + THERMAL_conduction_ID, & + THERMAL_CONDUCTION_label + use thermal_none + use thermal_conduction + + implicit none + integer(pInt), parameter :: FILEUNIT = 200_pInt + integer(pInt) :: & + g, & !< grain number + i, & !< integration point number + e, & !< element number + cMax, & !< maximum number of grains + iMax, & !< maximum number of integration points + eMax, & !< maximum number of elements + phase, & + s, & + p, & + instance,& + myNgrains + + integer(pInt), dimension(:,:), pointer :: thisSize + logical :: knownThermal + character(len=64), dimension(:,:), pointer :: thisOutput + character(len=32) :: outputName !< name of output, intermediate fix until HDF5 output is ready + +!-------------------------------------------------------------------------------------------------- +! parse from config file + if (.not. IO_open_jobFile_stat(FILEUNIT,material_localFileExt)) & ! no local material configuration present... + call IO_open_file(FILEUNIT,material_configFile) ! ... open material.config file + if (any(phase_thermal == THERMAL_none_ID)) call thermal_none_init(FILEUNIT) + if (any(phase_thermal == THERMAL_conduction_ID)) call thermal_conduction_init(FILEUNIT) + close(FILEUNIT) + + write(6,'(/,a)') ' <<<+- constitutive_thermal init -+>>>' + write(6,'(a)') ' $Id: constitutive_thermal.f90 3205 2014-06-17 06:54:49Z MPIE\m.diehl $' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() +#include "compilation_info.f90" + +!-------------------------------------------------------------------------------------------------- +! write description file for constitutive phase output + call IO_write_jobFile(FILEUNIT,'outputThermal') + do phase = 1_pInt,material_Nphase + instance = phase_thermalInstance(phase) ! which instance is present phase + knownThermal = .true. + select case(phase_thermal(phase)) ! split per constititution + case (THERMAL_none_ID) + outputName = THERMAL_NONE_label + thisOutput => null() + thisSize => null() + case (THERMAL_conduction_ID) + outputName = THERMAL_CONDUCTION_label + thisOutput => thermal_conduction_output + thisSize => thermal_conduction_sizePostResult + case default + knownThermal = .false. + end select + write(FILEUNIT,'(/,a,/)') '['//trim(phase_name(phase))//']' + if (knownThermal) then + write(FILEUNIT,'(a)') '(thermal)'//char(9)//trim(outputName) + if (phase_thermal(phase) /= THERMAL_none_ID) then + do e = 1_pInt,phase_Noutput(phase) + write(FILEUNIT,'(a,i4)') trim(thisOutput(e,instance))//char(9),thisSize(e,instance) + enddo + endif + endif + enddo + close(FILEUNIT) + +!-------------------------------------------------------------------------------------------------- +! allocation of states + cMax = homogenization_maxNgrains + iMax = mesh_maxNips + eMax = mesh_NcpElems + allocate(constitutive_thermal_sizePostResults(cMax,iMax,eMax), source=0_pInt) + + ElemLoop:do e = 1_pInt,mesh_NcpElems ! loop over elements + myNgrains = homogenization_Ngrains(mesh_element(3,e)) + IPloop:do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e))) ! loop over IPs + GrainLoop:do g = 1_pInt,myNgrains ! loop over grains + phase = material_phase(g,i,e) + instance = phase_thermalInstance(phase) + select case(phase_thermal(phase)) + case (THERMAL_conduction_ID) + constitutive_thermal_sizePostResults(g,i,e) = thermal_conduction_sizePostResults(instance) + end select + enddo GrainLoop + enddo IPloop + enddo ElemLoop + + constitutive_thermal_maxSizePostResults = maxval(constitutive_thermal_sizePostResults) + constitutive_thermal_maxSizeDotState = 0_pInt + do p = 1, size(thermalState) + constitutive_thermal_maxSizeDotState = max(constitutive_thermal_maxSizeDotState, thermalState(p)%sizeDotState) + enddo +end subroutine constitutive_thermal_init + + +!-------------------------------------------------------------------------------------------------- +!> @brief calls microstructure function of the different constitutive models +!-------------------------------------------------------------------------------------------------- +subroutine constitutive_thermal_microstructure(Tstar_v, Lp, ipc, ip, el) + use material, only: & + material_phase, & + phase_thermal, & + THERMAL_conduction_ID + use thermal_conduction, only: & + thermal_conduction_microstructure + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), intent(in), dimension(6) :: & + Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel) + real(pReal), intent(in), dimension(3,3) :: & + Lp + + select case (phase_thermal(material_phase(ipc,ip,el))) + case (THERMAL_conduction_ID) + call thermal_conduction_microstructure(Tstar_v, Lp, ipc, ip, el) + end select + +end subroutine constitutive_thermal_microstructure + + +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the rate of change of microstructure +!-------------------------------------------------------------------------------------------------- +subroutine constitutive_thermal_collectDotState(Tstar_v, Lp, ipc, ip, el) + use material, only: & + material_phase, & + phase_thermal, & + THERMAL_adiabatic_ID +! use thermal_conduction, only: & +! thermal_adiabatic_microstructure + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), intent(in), dimension(6) :: & + Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel) + real(pReal), intent(in), dimension(3,3) :: & + Lp + + select case (phase_thermal(material_phase(ipc,ip,el))) + case (THERMAL_adiabatic_ID) +! call thermal_adiabatic_dotState(Tstar_v, Lp, ipc, ip, el) + end select + +end subroutine constitutive_thermal_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_thermal_collectDeltaState(ipc, ip, el) + use material, only: & + material_phase, & + phase_thermal + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + + select case (phase_thermal(material_phase(ipc,ip,el))) + + end select + +end function constitutive_thermal_collectDeltaState + + +!-------------------------------------------------------------------------------------------------- +!> @brief returns array of constitutive results +!-------------------------------------------------------------------------------------------------- +function constitutive_thermal_postResults(ipc, ip, el) + use material, only: & + material_phase, & + phase_thermal, & + THERMAL_conduction_ID + use thermal_conduction, only: & + thermal_conduction_postResults + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), dimension(constitutive_thermal_sizePostResults(ipc,ip,el)) :: & + constitutive_thermal_postResults + + constitutive_thermal_postResults = 0.0_pReal + + select case (phase_thermal(material_phase(ipc,ip,el))) + case (THERMAL_conduction_ID) + constitutive_thermal_postResults = thermal_conduction_postResults(ipc,ip,el) + end select + +end function constitutive_thermal_postResults + + +end module constitutive_thermal diff --git a/code/damage_gradient.f90 b/code/damage_gradient.f90 new file mode 100644 index 000000000..7f7221fc2 --- /dev/null +++ b/code/damage_gradient.f90 @@ -0,0 +1,377 @@ +!-------------------------------------------------------------------------------------------------- +! $Id: damage_gradient.f90 3210 2014-06-17 15:24:44Z MPIE\m.diehl $ +!-------------------------------------------------------------------------------------------------- +!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH +!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH +!> @brief material subroutine incoprorating dislocation and twinning physics +!> @details to be done +!-------------------------------------------------------------------------------------------------- +module damage_gradient + use prec, only: & + pReal, & + pInt + + implicit none + private + integer(pInt), dimension(:), allocatable, public, protected :: & + damage_gradient_sizeDotState, & !< number of dotStates + damage_gradient_sizeState, & !< total number of microstructural state variables + damage_gradient_sizePostResults !< cumulative size of post results + + integer(pInt), dimension(:,:), allocatable, target, public :: & + damage_gradient_sizePostResult !< size of each post result output + + character(len=64), dimension(:,:), allocatable, target, public :: & + damage_gradient_output !< name of each post result output + + integer(pInt), dimension(:), allocatable, private :: & + damage_gradient_Noutput !< number of outputs per instance of this damage + + real(pReal), dimension(:), allocatable, private :: & + damage_gradient_crack_mobility + + enum, bind(c) + enumerator :: undefined_ID, & + local_damage_ID, & + gradient_damage_ID + end enum + integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & + damage_gradient_outputID !< ID of each post result output + + + public :: & + damage_gradient_init, & + damage_gradient_stateInit, & + damage_gradient_aTolState, & + damage_gradient_microstructure, & + damage_gradient_dotState, & + damage_gradient_postResults + +contains + + +!-------------------------------------------------------------------------------------------------- +!> @brief module initialization +!> @details reads in material parameters, allocates arrays, and does sanity checks +!-------------------------------------------------------------------------------------------------- +subroutine damage_gradient_init(fileUnit) + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) + use debug, only: & + debug_level,& + debug_constitutive,& + debug_levelBasic + use mesh, only: & + mesh_maxNips, & + mesh_NcpElems + use IO, only: & + IO_read, & + IO_lc, & + IO_getTag, & + IO_isBlank, & + IO_stringPos, & + IO_stringValue, & + IO_floatValue, & + IO_intValue, & + IO_warning, & + IO_error, & + IO_timeStamp, & + IO_EOF + use material, only: & + homogenization_maxNgrains, & + phase_damage, & + phase_damageInstance, & + phase_Noutput, & + DAMAGE_GRADIENT_label, & + DAMAGE_gradient_ID, & + material_phase, & + damageState, & + MATERIAL_partPhase + use numerics,only: & + numerics_integrator + + implicit none + integer(pInt), intent(in) :: fileUnit + + integer(pInt), parameter :: MAXNCHUNKS = 7_pInt + integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions + integer(pInt) :: maxNinstance,mySize=0_pInt,phase,instance,o + integer(pInt) :: sizeState, sizeDotState + integer(pInt) :: NofMyPhase + character(len=65536) :: & + tag = '', & + line = '' + + write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_GRADIENT_label//' init -+>>>' + write(6,'(a)') ' $Id: damage_gradient.f90 3210 2014-06-17 15:24:44Z MPIE\m.diehl $' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() +#include "compilation_info.f90" + + maxNinstance = int(count(phase_damage == DAMAGE_gradient_ID),pInt) + if (maxNinstance == 0_pInt) return + + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & + write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance + + allocate(damage_gradient_sizeDotState(maxNinstance), source=0_pInt) + allocate(damage_gradient_sizeState(maxNinstance), source=0_pInt) + allocate(damage_gradient_sizePostResults(maxNinstance), source=0_pInt) + allocate(damage_gradient_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt) + allocate(damage_gradient_output(maxval(phase_Noutput),maxNinstance)) + damage_gradient_output = '' + allocate(damage_gradient_outputID(maxval(phase_Noutput),maxNinstance), source=undefined_ID) + allocate(damage_gradient_Noutput(maxNinstance), source=0_pInt) + allocate(damage_gradient_crack_mobility(maxNinstance), source=0.0_pReal) + + rewind(fileUnit) + phase = 0_pInt + do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to + line = IO_read(fileUnit) + enddo + + parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part + line = IO_read(fileUnit) + if (IO_isBlank(line)) cycle ! skip empty lines + if (IO_getTag(line,'<','>') /= '') then ! stop at next part + line = IO_read(fileUnit, .true.) ! reset IO_read + exit + endif + if (IO_getTag(line,'[',']') /= '') then ! next phase section + phase = phase + 1_pInt ! advance phase section counter + cycle ! skip to next line + endif + if (phase > 0_pInt ) then; if (phase_damage(phase) == DAMAGE_gradient_ID) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran + instance = phase_damageInstance(phase) ! which instance of my damage is present phase + positions = IO_stringPos(line,MAXNCHUNKS) + tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key + select case(tag) + case ('(output)') + select case(IO_lc(IO_stringValue(line,positions,2_pInt))) + case ('local_damage') + damage_gradient_outputID(damage_gradient_Noutput(instance),instance) = local_damage_ID + damage_gradient_Noutput(instance) = damage_gradient_Noutput(instance) + 1_pInt + damage_gradient_output(damage_gradient_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,positions,2_pInt)) + case ('gradient_damage') + damage_gradient_outputID(damage_gradient_Noutput(instance),instance) = gradient_damage_ID + damage_gradient_Noutput(instance) = damage_gradient_Noutput(instance) + 1_pInt + damage_gradient_output(damage_gradient_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,positions,2_pInt)) + end select + + case ('crack_mobility') + damage_gradient_crack_mobility(instance) = IO_floatValue(line,positions,2_pInt) + end select + endif; endif + enddo parsingFile + + initializeInstances: do phase = 1_pInt, size(phase_damage) + if (phase_damage(phase) == DAMAGE_gradient_ID) then + NofMyPhase=count(material_phase==phase) + instance = phase_damageInstance(phase) + damage_gradient_sizeDotState(instance) = 1_pInt + damage_gradient_sizeState(instance) = 3_pInt + +!-------------------------------------------------------------------------------------------------- +! Determine size of postResults array + outputsLoop: do o = 1_pInt,damage_gradient_Noutput(instance) + select case(damage_gradient_outputID(o,instance)) + case(local_damage_ID, & + gradient_damage_ID & + ) + mySize = 1_pInt + end select + + if (mySize > 0_pInt) then ! any meaningful output found + damage_gradient_sizePostResult(o,instance) = mySize + damage_gradient_sizePostResults(instance) = damage_gradient_sizePostResults(instance) + mySize + endif + enddo outputsLoop +! Determine size of state array + sizeDotState = damage_gradient_sizeDotState(instance) + sizeState = damage_gradient_sizeState (instance) + + damageState(phase)%sizeState = sizeState + damageState(phase)%sizeDotState = sizeDotState + allocate(damageState(phase)%aTolState (sizeState), source=0.0_pReal) + allocate(damageState(phase)%state0 (sizeState,NofMyPhase), source=0.0_pReal) + allocate(damageState(phase)%partionedState0 (sizeState,NofMyPhase), source=0.0_pReal) + allocate(damageState(phase)%subState0 (sizeState,NofMyPhase), source=0.0_pReal) + allocate(damageState(phase)%state (sizeState,NofMyPhase), source=0.0_pReal) + allocate(damageState(phase)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) + + allocate(damageState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) + allocate(damageState(phase)%deltaState (sizeDotState,NofMyPhase), source=0.0_pReal) + allocate(damageState(phase)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) + if (any(numerics_integrator == 1_pInt)) then + allocate(damageState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) + allocate(damageState(phase)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal) + endif + if (any(numerics_integrator == 4_pInt)) & + allocate(damageState(phase)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal) + if (any(numerics_integrator == 5_pInt)) & + allocate(damageState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal) + + call damage_gradient_stateInit(phase,instance) + call damage_gradient_aTolState(phase,instance) + endif + + enddo initializeInstances +end subroutine damage_gradient_init + +!-------------------------------------------------------------------------------------------------- +!> @brief sets the relevant NEW state values for a given instance of this damage +!-------------------------------------------------------------------------------------------------- +subroutine damage_gradient_stateInit(phase,instance) + use material, only: & + damageState + + implicit none + integer(pInt), intent(in) :: instance !< number specifying the instance of the damage + integer(pInt), intent(in) :: phase !< number specifying the phase of the damage + + real(pReal), dimension(damageState(phase)%sizeState) :: tempState + + tempState(1) = 0.0_pReal + tempState(2:3) = 1.0_pReal + damageState(phase)%state = spread(tempState,2,size(damageState(phase)%state(1,:))) + damageState(phase)%state0 = damageState(phase)%state + damageState(phase)%partionedState0 = damageState(phase)%state +end subroutine damage_gradient_stateInit + +!-------------------------------------------------------------------------------------------------- +!> @brief sets the relevant state values for a given instance of this damage +!-------------------------------------------------------------------------------------------------- +subroutine damage_gradient_aTolState(phase,instance) + use material, only: & + damageState + + implicit none + integer(pInt), intent(in) :: & + phase, & + instance ! number specifying the current instance of the damage + real(pReal), dimension(damageState(phase)%sizeState) :: tempTol + + tempTol = 0.0_pReal + damageState(phase)%aTolState = tempTol +end subroutine damage_gradient_aTolState + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates derived quantities from state +!-------------------------------------------------------------------------------------------------- +subroutine damage_gradient_microstructure(Tstar_v, Fe, ipc, ip, el) + use material, only: & + mappingConstitutive, & + phase_damageInstance, & + damageState + use mesh, only: & + charLength + use math, only: & + math_Mandel6to33, & + math_mul33x33, & + math_transpose33, & + math_I3 + use lattice, only: & + lattice_surfaceEnergy33 + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(6) :: & + Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel) + real(pReal), intent(in), dimension(3,3) :: & + Fe + integer(pInt) :: & + instance, phase, constituent + real(pReal) :: & + damage + + phase = mappingConstitutive(2,ipc,ip,el) + constituent = mappingConstitutive(1,ipc,ip,el) + instance = phase_damageInstance(phase) + damage = damageState(phase)%state(3,constituent)*damageState(phase)%state(3,constituent) + + damageState(phase)%state(2,constituent) = & + min(1.0_pReal, & + 2.0_pReal*charLength*maxval(lattice_surfaceEnergy33(1:3,1:3,phase))/ & + (0.125_pReal*sum(math_Mandel6to33(Tstar_v/damage)*(math_mul33x33(math_transpose33(Fe),Fe)-math_I3)) + & + 0.5_pReal*damageState(phase)%state(1,constituent)) & + ) + +end subroutine damage_gradient_microstructure + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates derived quantities from state +!-------------------------------------------------------------------------------------------------- +subroutine damage_gradient_dotState(Tstar_v, Lp, ipc, ip, el) + use material, only: & + mappingConstitutive, & + phase_damageInstance, & + damageState + use math, only: & + math_Mandel6to33 + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(6) :: & + Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel) + real(pReal), intent(in), dimension(3,3) :: & + Lp + integer(pInt) :: & + instance, phase, constituent + + phase = mappingConstitutive(2,ipc,ip,el) + constituent = mappingConstitutive(1,ipc,ip,el) + instance = phase_damageInstance(phase) + + damageState(phase)%dotState(1,constituent) = & + sum(abs(math_Mandel6to33(Tstar_v)*Lp)) + +end subroutine damage_gradient_dotState + + +!-------------------------------------------------------------------------------------------------- +!> @brief return array of constitutive results +!-------------------------------------------------------------------------------------------------- +function damage_gradient_postResults(ipc,ip,el) + use material, only: & + mappingConstitutive, & + phase_damageInstance,& + damageState + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), dimension(damage_gradient_sizePostResults(phase_damageInstance(mappingConstitutive(2,ipc,ip,el)))) :: & + damage_gradient_postResults + + integer(pInt) :: & + instance, phase, constituent, o, c + + phase = mappingConstitutive(2,ipc,ip,el) + constituent = mappingConstitutive(1,ipc,ip,el) + instance = phase_damageInstance(phase) + + c = 0_pInt + damage_gradient_postResults = 0.0_pReal + + do o = 1_pInt,damage_gradient_Noutput(instance) + select case(damage_gradient_outputID(o,instance)) + + case (local_damage_ID) + damage_gradient_postResults(c+1_pInt) = damageState(phase)%state(2,constituent) + c = c + 1 + case (gradient_damage_ID) + damage_gradient_postResults(c+1_pInt) = damageState(phase)%state(3,constituent) + c = c + 1 + end select + enddo +end function damage_gradient_postResults + +end module damage_gradient diff --git a/code/damage_none.f90 b/code/damage_none.f90 new file mode 100644 index 000000000..74e2340ad --- /dev/null +++ b/code/damage_none.f90 @@ -0,0 +1,106 @@ +!-------------------------------------------------------------------------------------------------- +! $Id: damage_none.f90 3148 2014-05-27 14:46:03Z MPIE\m.diehl $ +!-------------------------------------------------------------------------------------------------- +!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH +!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH +!> @brief material subroutine for purely elastic material +!-------------------------------------------------------------------------------------------------- +module damage_none + use prec, only: & + pInt + + implicit none + private + integer(pInt), dimension(:), allocatable, public, protected :: & + damage_none_sizeDotState, & + damage_none_sizeState, & + damage_none_sizePostResults + + integer(pInt), dimension(:,:), allocatable, target, public :: & + damage_none_sizePostResult !< size of each post result output + + public :: & + damage_none_init + +contains + + +!-------------------------------------------------------------------------------------------------- +!> @brief module initialization +!> @details reads in material parameters, allocates arrays, and does sanity checks +!-------------------------------------------------------------------------------------------------- +subroutine damage_none_init(fileUnit) + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) + use debug, only: & + debug_level, & + debug_constitutive, & + debug_levelBasic + use IO, only: & + IO_timeStamp + use numerics, only: & + numerics_integrator + use material, only: & + phase_damage, & + phase_Noutput, & + DAMAGE_NONE_label, & + material_phase, & + damageState, & + DAMAGE_NONE_ID, & + MATERIAL_partPhase + + implicit none + + integer(pInt), intent(in) :: fileUnit + integer(pInt) :: & + maxNinstance, & + phase, & + NofMyPhase, & + sizeState, & + sizeDotState + + write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_NONE_label//' init -+>>>' + write(6,'(a)') ' $Id: damage_none.f90 3148 2014-05-27 14:46:03Z MPIE\m.diehl $' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() +#include "compilation_info.f90" + + maxNinstance = int(count(phase_damage == DAMAGE_NONE_ID),pInt) + if (maxNinstance == 0_pInt) return + + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & + write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance + +#ifdef NEWSTATE + initializeInstances: do phase = 1_pInt, size(phase_damage) + NofMyPhase=count(material_phase==phase) + if (phase_damage(phase) == DAMAGE_none_ID .and. NofMyPhase/=0) then + sizeState = 0_pInt + damageState(phase)%sizeState = sizeState + sizeDotState = sizeState + damageState(phase)%sizeDotState = sizeDotState + allocate(damageState(phase)%state0 (sizeState,NofMyPhase)) + allocate(damageState(phase)%partionedState0(sizeState,NofMyPhase)) + allocate(damageState(phase)%subState0 (sizeState,NofMyPhase)) + allocate(damageState(phase)%state (sizeState,NofMyPhase)) + allocate(damageState(phase)%state_backup (sizeState,NofMyPhase)) + allocate(damageState(phase)%aTolState (NofMyPhase)) + allocate(damageState(phase)%dotState (sizeDotState,NofMyPhase)) + allocate(damageState(phase)%dotState_backup(sizeDotState,NofMyPhase)) + if (any(numerics_integrator == 1_pInt)) then + allocate(damageState(phase)%previousDotState (sizeDotState,NofMyPhase)) + allocate(damageState(phase)%previousDotState2 (sizeDotState,NofMyPhase)) + endif + if (any(numerics_integrator == 4_pInt)) & + allocate(damageState(phase)%RK4dotState (sizeDotState,NofMyPhase)) + if (any(numerics_integrator == 5_pInt)) & + allocate(damageState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase)) + endif + enddo initializeInstances +#else + allocate(damage_none_sizeDotState(maxNinstance), source=1_pInt) + allocate(damage_none_sizeState(maxNinstance), source=1_pInt) +#endif + allocate(damage_none_sizePostResults(maxNinstance), source=0_pInt) + +end subroutine damage_none_init + +end module damage_none diff --git a/code/thermal_adiabatic.f90 b/code/thermal_adiabatic.f90 new file mode 100644 index 000000000..0a98cba97 --- /dev/null +++ b/code/thermal_adiabatic.f90 @@ -0,0 +1,326 @@ +!-------------------------------------------------------------------------------------------------- +! $Id: thermal_adiabatic.f90 3210 2014-06-17 15:24:44Z MPIE\m.diehl $ +!-------------------------------------------------------------------------------------------------- +!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH +!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH +!> @brief material subroutine incoprorating dislocation and twinning physics +!> @details to be done +!-------------------------------------------------------------------------------------------------- +module thermal_adiabatic + use prec, only: & + pReal, & + pInt + + implicit none + private + integer(pInt), dimension(:), allocatable, public, protected :: & + thermal_adiabatic_sizeDotState, & !< number of dotStates + thermal_adiabatic_sizeState, & !< total number of microstructural state variables + thermal_adiabatic_sizePostResults !< cumulative size of post results + + integer(pInt), dimension(:,:), allocatable, target, public :: & + thermal_adiabatic_sizePostResult !< size of each post result output + + character(len=64), dimension(:,:), allocatable, target, public :: & + thermal_adiabatic_output !< name of each post result output + + integer(pInt), dimension(:), allocatable, private :: & + thermal_adiabatic_Noutput !< number of outputs per instance of this damage + + real(pReal), dimension(:), allocatable, private :: & + thermal_adiabatic_specific_heat, & + thermal_adiabatic_density + + enum, bind(c) + enumerator :: undefined_ID, & + temperature_ID + end enum + integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & + thermal_adiabatic_outputID !< ID of each post result output + + + public :: & + thermal_adiabatic_init, & + thermal_adiabatic_stateInit, & + thermal_adiabatic_aTolState, & + thermal_adiabatic_dotState, & + thermal_adiabatic_postResults + +contains + + +!-------------------------------------------------------------------------------------------------- +!> @brief module initialization +!> @details reads in material parameters, allocates arrays, and does sanity checks +!-------------------------------------------------------------------------------------------------- +subroutine thermal_adiabatic_init(fileUnit) + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) + use debug, only: & + debug_level,& + debug_constitutive,& + debug_levelBasic + use mesh, only: & + mesh_maxNips, & + mesh_NcpElems + use IO, only: & + IO_read, & + IO_lc, & + IO_getTag, & + IO_isBlank, & + IO_stringPos, & + IO_stringValue, & + IO_floatValue, & + IO_intValue, & + IO_warning, & + IO_error, & + IO_timeStamp, & + IO_EOF + use material, only: & + homogenization_maxNgrains, & + phase_thermal, & + phase_thermalInstance, & + phase_Noutput, & + THERMAL_ADIABATIC_label, & + THERMAL_adiabatic_ID, & + material_phase, & + thermalState, & + MATERIAL_partPhase + use numerics,only: & + numerics_integrator + + implicit none + integer(pInt), intent(in) :: fileUnit + + integer(pInt), parameter :: MAXNCHUNKS = 7_pInt + integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions + integer(pInt) :: maxNinstance,mySize=0_pInt,phase,instance,o + integer(pInt) :: sizeState, sizeDotState + integer(pInt) :: NofMyPhase + character(len=65536) :: & + tag = '', & + line = '' + + write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_ADIABATIC_label//' init -+>>>' + write(6,'(a)') ' $Id: thermal_adiabatic.f90 3210 2014-06-17 15:24:44Z MPIE\m.diehl $' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() +#include "compilation_info.f90" + + maxNinstance = int(count(phase_thermal == THERMAL_adiabatic_ID),pInt) + if (maxNinstance == 0_pInt) return + + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & + write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance + + allocate(thermal_adiabatic_sizeDotState(maxNinstance), source=0_pInt) + allocate(thermal_adiabatic_sizeState(maxNinstance), source=0_pInt) + allocate(thermal_adiabatic_sizePostResults(maxNinstance), source=0_pInt) + allocate(thermal_adiabatic_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt) + allocate(thermal_adiabatic_output(maxval(phase_Noutput),maxNinstance)) + thermal_adiabatic_output = '' + allocate(thermal_adiabatic_outputID(maxval(phase_Noutput),maxNinstance), source=undefined_ID) + allocate(thermal_adiabatic_Noutput(maxNinstance), source=0_pInt) + allocate(thermal_adiabatic_specific_heat(maxNinstance), source=0.0_pReal) + allocate(thermal_adiabatic_density(maxNinstance), source=0.0_pReal) + + rewind(fileUnit) + phase = 0_pInt + do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to + line = IO_read(fileUnit) + enddo + + parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part + line = IO_read(fileUnit) + if (IO_isBlank(line)) cycle ! skip empty lines + if (IO_getTag(line,'<','>') /= '') then ! stop at next part + line = IO_read(fileUnit, .true.) ! reset IO_read + exit + endif + if (IO_getTag(line,'[',']') /= '') then ! next phase section + phase = phase + 1_pInt ! advance phase section counter + cycle ! skip to next line + endif + if (phase > 0_pInt ) then; if (phase_thermal(phase) == THERMAL_adiabatic_ID) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran + instance = phase_thermalInstance(phase) ! which instance of my thermal is present phase + positions = IO_stringPos(line,MAXNCHUNKS) + tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key + select case(tag) + case ('(output)') + select case(IO_lc(IO_stringValue(line,positions,2_pInt))) + case ('temperature') + thermal_adiabatic_outputID(thermal_adiabatic_Noutput(instance),instance) = temperature_ID + thermal_adiabatic_Noutput(instance) = thermal_adiabatic_Noutput(instance) + 1_pInt + thermal_adiabatic_output(thermal_adiabatic_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,positions,2_pInt)) + end select + + case ('specific_heat') + thermal_adiabatic_specific_heat(instance) = IO_floatValue(line,positions,2_pInt) + case ('density') + thermal_adiabatic_density(instance) = IO_floatValue(line,positions,2_pInt) + end select + endif; endif + enddo parsingFile + + initializeInstances: do phase = 1_pInt, size(phase_thermal) + if (phase_thermal(phase) == THERMAL_adiabatic_ID) then + NofMyPhase=count(material_phase==phase) + instance = phase_thermalInstance(phase) + thermal_adiabatic_sizeDotState(instance) = 1_pInt + thermal_adiabatic_sizeState(instance) = 1_pInt + +!-------------------------------------------------------------------------------------------------- +! Determine size of postResults array + outputsLoop: do o = 1_pInt,thermal_adiabatic_Noutput(instance) + select case(thermal_adiabatic_outputID(o,instance)) + case(temperature_ID) + mySize = 1_pInt + end select + + if (mySize > 0_pInt) then ! any meaningful output found + thermal_adiabatic_sizePostResult(o,instance) = mySize + thermal_adiabatic_sizePostResults(instance) = thermal_adiabatic_sizePostResults(instance) + mySize + endif + enddo outputsLoop +! Determine size of state array + sizeDotState = thermal_adiabatic_sizeDotState(instance) + sizeState = thermal_adiabatic_sizeState (instance) + + thermalState(phase)%sizeState = sizeState + thermalState(phase)%sizeDotState = sizeDotState + allocate(thermalState(phase)%aTolState (sizeState), source=0.0_pReal) + allocate(thermalState(phase)%state0 (sizeState,NofMyPhase), source=0.0_pReal) + allocate(thermalState(phase)%partionedState0 (sizeState,NofMyPhase), source=0.0_pReal) + allocate(thermalState(phase)%subState0 (sizeState,NofMyPhase), source=0.0_pReal) + allocate(thermalState(phase)%state (sizeState,NofMyPhase), source=0.0_pReal) + allocate(thermalState(phase)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) + + allocate(thermalState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) + allocate(thermalState(phase)%deltaState (sizeDotState,NofMyPhase), source=0.0_pReal) + allocate(thermalState(phase)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) + if (any(numerics_integrator == 1_pInt)) then + allocate(thermalState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) + allocate(thermalState(phase)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal) + endif + if (any(numerics_integrator == 4_pInt)) & + allocate(thermalState(phase)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal) + if (any(numerics_integrator == 5_pInt)) & + allocate(thermalState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal) + + call thermal_adiabatic_stateInit(phase,instance) + call thermal_adiabatic_aTolState(phase,instance) + endif + + enddo initializeInstances +end subroutine thermal_adiabatic_init + +!-------------------------------------------------------------------------------------------------- +!> @brief sets the relevant NEW state values for a given instance of this thermal +!-------------------------------------------------------------------------------------------------- +subroutine thermal_adiabatic_stateInit(phase,instance) + use material, only: & + thermalState + use lattice, only: & + lattice_referenceTemperature + + implicit none + integer(pInt), intent(in) :: instance !< number specifying the instance of the thermal + integer(pInt), intent(in) :: phase !< number specifying the phase of the thermal + + real(pReal), dimension(thermalState(phase)%sizeState) :: tempState + + tempState(1) = lattice_referenceTemperature(phase) + thermalState(phase)%state = spread(tempState,2,size(thermalState(phase)%state(1,:))) + thermalState(phase)%state0 = thermalState(phase)%state + thermalState(phase)%partionedState0 = thermalState(phase)%state +end subroutine thermal_adiabatic_stateInit + +!-------------------------------------------------------------------------------------------------- +!> @brief sets the relevant state values for a given instance of this thermal +!-------------------------------------------------------------------------------------------------- +subroutine thermal_adiabatic_aTolState(phase,instance) + use material, only: & + thermalState + + implicit none + integer(pInt), intent(in) :: & + phase, & + instance ! number specifying the current instance of the thermal + real(pReal), dimension(thermalState(phase)%sizeState) :: tempTol + + tempTol = 0.0_pReal + thermalState(phase)%aTolState = tempTol +end subroutine thermal_adiabatic_aTolState + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates derived quantities from state +!-------------------------------------------------------------------------------------------------- +subroutine thermal_adiabatic_dotState(Tstar_v, Lp, ipc, ip, el) + use material, only: & + mappingConstitutive, & + phase_thermalInstance, & + thermalState + use math, only: & + math_Mandel6to33 + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(6) :: & + Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel) + real(pReal), intent(in), dimension(3,3) :: & + Lp + integer(pInt) :: & + instance, phase, constituent + + phase = mappingConstitutive(2,ipc,ip,el) + constituent = mappingConstitutive(1,ipc,ip,el) + instance = phase_thermalInstance(phase) + + thermalState(phase)%dotState(1,constituent) = & + 0.95_pReal & + * sum(abs(math_Mandel6to33(Tstar_v)*Lp)) & + / (thermal_adiabatic_density(phase)*thermal_adiabatic_specific_heat(phase)) + +end subroutine thermal_adiabatic_dotState + + +!-------------------------------------------------------------------------------------------------- +!> @brief return array of constitutive results +!-------------------------------------------------------------------------------------------------- +function thermal_adiabatic_postResults(ipc,ip,el) + use material, only: & + mappingConstitutive, & + phase_thermalInstance, & + thermalState + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), dimension(thermal_adiabatic_sizePostResults(phase_thermalInstance(mappingConstitutive(2,ipc,ip,el)))) :: & + thermal_adiabatic_postResults + + integer(pInt) :: & + instance, phase, constituent, o, c + + phase = mappingConstitutive(2,ipc,ip,el) + constituent = mappingConstitutive(1,ipc,ip,el) + instance = phase_thermalInstance(phase) + + c = 0_pInt + thermal_adiabatic_postResults = 0.0_pReal + + do o = 1_pInt,thermal_adiabatic_Noutput(instance) + select case(thermal_adiabatic_outputID(o,instance)) + + case (temperature_ID) + thermal_adiabatic_postResults(c+1_pInt) = thermalState(phase)%state(1,constituent) + c = c + 1 + end select + enddo +end function thermal_adiabatic_postResults + +end module thermal_adiabatic diff --git a/code/thermal_conduction.f90 b/code/thermal_conduction.f90 new file mode 100644 index 000000000..b11ebd51d --- /dev/null +++ b/code/thermal_conduction.f90 @@ -0,0 +1,325 @@ +!-------------------------------------------------------------------------------------------------- +! $Id: thermal_conduction.f90 3210 2014-06-17 15:24:44Z MPIE\m.diehl $ +!-------------------------------------------------------------------------------------------------- +!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH +!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH +!> @brief material subroutine incoprorating dislocation and twinning physics +!> @details to be done +!-------------------------------------------------------------------------------------------------- +module thermal_conduction + use prec, only: & + pReal, & + pInt + + implicit none + private + integer(pInt), dimension(:), allocatable, public, protected :: & + thermal_conduction_sizeDotState, & !< number of dotStates + thermal_conduction_sizeState, & !< total number of microstructural state variables + thermal_conduction_sizePostResults !< cumulative size of post results + + integer(pInt), dimension(:,:), allocatable, target, public :: & + thermal_conduction_sizePostResult !< size of each post result output + + character(len=64), dimension(:,:), allocatable, target, public :: & + thermal_conduction_output !< name of each post result output + + integer(pInt), dimension(:), allocatable, private :: & + thermal_conduction_Noutput !< number of outputs per instance of this damage + + real(pReal), dimension(:), allocatable, private :: & + thermal_conduction_specific_heat, & + thermal_conduction_density + + enum, bind(c) + enumerator :: undefined_ID, & + temperature_ID + end enum + integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & + thermal_conduction_outputID !< ID of each post result output + + + public :: & + thermal_conduction_init, & + thermal_conduction_stateInit, & + thermal_conduction_aTolState, & + thermal_conduction_microstructure, & + thermal_conduction_postResults + +contains + + +!-------------------------------------------------------------------------------------------------- +!> @brief module initialization +!> @details reads in material parameters, allocates arrays, and does sanity checks +!-------------------------------------------------------------------------------------------------- +subroutine thermal_conduction_init(fileUnit) + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) + use debug, only: & + debug_level,& + debug_constitutive,& + debug_levelBasic + use mesh, only: & + mesh_maxNips, & + mesh_NcpElems + use IO, only: & + IO_read, & + IO_lc, & + IO_getTag, & + IO_isBlank, & + IO_stringPos, & + IO_stringValue, & + IO_floatValue, & + IO_intValue, & + IO_warning, & + IO_error, & + IO_timeStamp, & + IO_EOF + use material, only: & + homogenization_maxNgrains, & + phase_thermal, & + phase_thermalInstance, & + phase_Noutput, & + THERMAL_CONDUCTION_label, & + THERMAL_conduction_ID, & + material_phase, & + thermalState, & + MATERIAL_partPhase + use numerics,only: & + numerics_integrator + + implicit none + integer(pInt), intent(in) :: fileUnit + + integer(pInt), parameter :: MAXNCHUNKS = 7_pInt + integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions + integer(pInt) :: maxNinstance,mySize=0_pInt,phase,instance,o + integer(pInt) :: sizeState, sizeDotState + integer(pInt) :: NofMyPhase + character(len=65536) :: & + tag = '', & + line = '' + + write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_CONDUCTION_label//' init -+>>>' + write(6,'(a)') ' $Id: thermal_conduction.f90 3210 2014-06-17 15:24:44Z MPIE\m.diehl $' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() +#include "compilation_info.f90" + + maxNinstance = int(count(phase_thermal == THERMAL_conduction_ID),pInt) + if (maxNinstance == 0_pInt) return + + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & + write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance + + allocate(thermal_conduction_sizeDotState(maxNinstance), source=0_pInt) + allocate(thermal_conduction_sizeState(maxNinstance), source=0_pInt) + allocate(thermal_conduction_sizePostResults(maxNinstance), source=0_pInt) + allocate(thermal_conduction_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt) + allocate(thermal_conduction_output(maxval(phase_Noutput),maxNinstance)) + thermal_conduction_output = '' + allocate(thermal_conduction_outputID(maxval(phase_Noutput),maxNinstance), source=undefined_ID) + allocate(thermal_conduction_Noutput(maxNinstance), source=0_pInt) + allocate(thermal_conduction_specific_heat(maxNinstance), source=0.0_pReal) + allocate(thermal_conduction_density(maxNinstance), source=0.0_pReal) + + rewind(fileUnit) + phase = 0_pInt + do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to + line = IO_read(fileUnit) + enddo + + parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part + line = IO_read(fileUnit) + if (IO_isBlank(line)) cycle ! skip empty lines + if (IO_getTag(line,'<','>') /= '') then ! stop at next part + line = IO_read(fileUnit, .true.) ! reset IO_read + exit + endif + if (IO_getTag(line,'[',']') /= '') then ! next phase section + phase = phase + 1_pInt ! advance phase section counter + cycle ! skip to next line + endif + if (phase > 0_pInt ) then; if (phase_thermal(phase) == THERMAL_conduction_ID) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran + instance = phase_thermalInstance(phase) ! which instance of my thermal is present phase + positions = IO_stringPos(line,MAXNCHUNKS) + tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key + select case(tag) + case ('(output)') + select case(IO_lc(IO_stringValue(line,positions,2_pInt))) + case ('temperature') + thermal_conduction_outputID(thermal_conduction_Noutput(instance),instance) = temperature_ID + thermal_conduction_Noutput(instance) = thermal_conduction_Noutput(instance) + 1_pInt + thermal_conduction_output(thermal_conduction_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,positions,2_pInt)) + end select + + case ('specific_heat') + thermal_conduction_specific_heat(instance) = IO_floatValue(line,positions,2_pInt) + case ('density') + thermal_conduction_density(instance) = IO_floatValue(line,positions,2_pInt) + end select + endif; endif + enddo parsingFile + + initializeInstances: do phase = 1_pInt, size(phase_thermal) + if (phase_thermal(phase) == THERMAL_conduction_ID) then + NofMyPhase=count(material_phase==phase) + instance = phase_thermalInstance(phase) + thermal_conduction_sizeDotState(instance) = 0_pInt + thermal_conduction_sizeState(instance) = 2_pInt + +!-------------------------------------------------------------------------------------------------- +! Determine size of postResults array + outputsLoop: do o = 1_pInt,thermal_conduction_Noutput(instance) + select case(thermal_conduction_outputID(o,instance)) + case(temperature_ID) + mySize = 1_pInt + end select + + if (mySize > 0_pInt) then ! any meaningful output found + thermal_conduction_sizePostResult(o,instance) = mySize + thermal_conduction_sizePostResults(instance) = thermal_conduction_sizePostResults(instance) + mySize + endif + enddo outputsLoop +! Determine size of state array + sizeDotState = thermal_conduction_sizeDotState(instance) + sizeState = thermal_conduction_sizeState (instance) + + thermalState(phase)%sizeState = sizeState + thermalState(phase)%sizeDotState = sizeDotState + allocate(thermalState(phase)%aTolState (sizeState), source=0.0_pReal) + allocate(thermalState(phase)%state0 (sizeState,NofMyPhase), source=0.0_pReal) + allocate(thermalState(phase)%partionedState0 (sizeState,NofMyPhase), source=0.0_pReal) + allocate(thermalState(phase)%subState0 (sizeState,NofMyPhase), source=0.0_pReal) + allocate(thermalState(phase)%state (sizeState,NofMyPhase), source=0.0_pReal) + allocate(thermalState(phase)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) + + allocate(thermalState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) + allocate(thermalState(phase)%deltaState (sizeDotState,NofMyPhase), source=0.0_pReal) + allocate(thermalState(phase)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) + if (any(numerics_integrator == 1_pInt)) then + allocate(thermalState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) + allocate(thermalState(phase)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal) + endif + if (any(numerics_integrator == 4_pInt)) & + allocate(thermalState(phase)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal) + if (any(numerics_integrator == 5_pInt)) & + allocate(thermalState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal) + + call thermal_conduction_stateInit(phase,instance) + call thermal_conduction_aTolState(phase,instance) + endif + + enddo initializeInstances +end subroutine thermal_conduction_init + +!-------------------------------------------------------------------------------------------------- +!> @brief sets the relevant NEW state values for a given instance of this thermal +!-------------------------------------------------------------------------------------------------- +subroutine thermal_conduction_stateInit(phase,instance) + use material, only: & + thermalState + use lattice, only: & + lattice_referenceTemperature + + implicit none + integer(pInt), intent(in) :: instance !< number specifying the instance of the thermal + integer(pInt), intent(in) :: phase !< number specifying the phase of the thermal + + real(pReal), dimension(thermalState(phase)%sizeState) :: tempState + + tempState(1) = 0.0_pReal + tempState(2) = lattice_referenceTemperature(phase) + thermalState(phase)%state = spread(tempState,2,size(thermalState(phase)%state(1,:))) + thermalState(phase)%state0 = thermalState(phase)%state + thermalState(phase)%partionedState0 = thermalState(phase)%state +end subroutine thermal_conduction_stateInit + +!-------------------------------------------------------------------------------------------------- +!> @brief sets the relevant state values for a given instance of this thermal +!-------------------------------------------------------------------------------------------------- +subroutine thermal_conduction_aTolState(phase,instance) + use material, only: & + thermalState + + implicit none + integer(pInt), intent(in) :: & + phase, & + instance ! number specifying the current instance of the thermal + real(pReal), dimension(thermalState(phase)%sizeState) :: tempTol + + tempTol = 0.0_pReal + thermalState(phase)%aTolState = tempTol +end subroutine thermal_conduction_aTolState + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates derived quantities from state +!-------------------------------------------------------------------------------------------------- +subroutine thermal_conduction_microstructure(Tstar_v, Lp, ipc, ip, el) + use material, only: & + mappingConstitutive, & + phase_thermalInstance, & + thermalState + use math, only: & + math_Mandel6to33 + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(6) :: & + Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel) + real(pReal), intent(in), dimension(3,3) :: & + Lp + integer(pInt) :: & + instance, phase, constituent + + phase = mappingConstitutive(2,ipc,ip,el) + constituent = mappingConstitutive(1,ipc,ip,el) + instance = phase_thermalInstance(phase) + + thermalState(phase)%state(1,constituent) = & + sum(abs(math_Mandel6to33(Tstar_v)*Lp)) + +end subroutine thermal_conduction_microstructure + + +!-------------------------------------------------------------------------------------------------- +!> @brief return array of constitutive results +!-------------------------------------------------------------------------------------------------- +function thermal_conduction_postResults(ipc,ip,el) + use material, only: & + mappingConstitutive, & + phase_thermalInstance, & + thermalState + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), dimension(thermal_conduction_sizePostResults(phase_thermalInstance(mappingConstitutive(2,ipc,ip,el)))) :: & + thermal_conduction_postResults + + integer(pInt) :: & + instance, phase, constituent, o, c + + phase = mappingConstitutive(2,ipc,ip,el) + constituent = mappingConstitutive(1,ipc,ip,el) + instance = phase_thermalInstance(phase) + + c = 0_pInt + thermal_conduction_postResults = 0.0_pReal + + do o = 1_pInt,thermal_conduction_Noutput(instance) + select case(thermal_conduction_outputID(o,instance)) + + case (temperature_ID) + thermal_conduction_postResults(c+1_pInt) = thermalState(phase)%state(2,constituent) + c = c + 1 + end select + enddo +end function thermal_conduction_postResults + +end module thermal_conduction diff --git a/code/thermal_none.f90 b/code/thermal_none.f90 new file mode 100644 index 000000000..f757b9880 --- /dev/null +++ b/code/thermal_none.f90 @@ -0,0 +1,106 @@ +!-------------------------------------------------------------------------------------------------- +! $Id: thermal_none.f90 3148 2014-05-27 14:46:03Z MPIE\m.diehl $ +!-------------------------------------------------------------------------------------------------- +!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH +!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH +!> @brief material subroutine for purely elastic material +!-------------------------------------------------------------------------------------------------- +module thermal_none + use prec, only: & + pInt + + implicit none + private + integer(pInt), dimension(:), allocatable, public, protected :: & + thermal_none_sizeDotState, & + thermal_none_sizeState, & + thermal_none_sizePostResults + + integer(pInt), dimension(:,:), allocatable, target, public :: & + thermal_none_sizePostResult !< size of each post result output + + public :: & + thermal_none_init + +contains + + +!-------------------------------------------------------------------------------------------------- +!> @brief module initialization +!> @details reads in material parameters, allocates arrays, and does sanity checks +!-------------------------------------------------------------------------------------------------- +subroutine thermal_none_init(fileUnit) + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) + use debug, only: & + debug_level, & + debug_constitutive, & + debug_levelBasic + use IO, only: & + IO_timeStamp + use numerics, only: & + numerics_integrator + use material, only: & + phase_thermal, & + phase_Noutput, & + THERMAL_NONE_label, & + material_phase, & + thermalState, & + THERMAL_NONE_ID, & + MATERIAL_partPhase + + implicit none + + integer(pInt), intent(in) :: fileUnit + integer(pInt) :: & + maxNinstance, & + phase, & + NofMyPhase, & + sizeState, & + sizeDotState + + write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_NONE_label//' init -+>>>' + write(6,'(a)') ' $Id: thermal_none.f90 3148 2014-05-27 14:46:03Z MPIE\m.diehl $' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() +#include "compilation_info.f90" + + maxNinstance = int(count(phase_thermal == THERMAL_NONE_ID),pInt) + if (maxNinstance == 0_pInt) return + + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & + write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance + +#ifdef NEWSTATE + initializeInstances: do phase = 1_pInt, size(phase_thermal) + NofMyPhase=count(material_phase==phase) + if (phase_thermal(phase) == THERMAL_none_ID .and. NofMyPhase/=0) then + sizeState = 0_pInt + thermalState(phase)%sizeState = sizeState + sizeDotState = sizeState + thermalState(phase)%sizeDotState = sizeDotState + allocate(thermalState(phase)%state0 (sizeState,NofMyPhase)) + allocate(thermalState(phase)%partionedState0(sizeState,NofMyPhase)) + allocate(thermalState(phase)%subState0 (sizeState,NofMyPhase)) + allocate(thermalState(phase)%state (sizeState,NofMyPhase)) + allocate(thermalState(phase)%state_backup (sizeState,NofMyPhase)) + allocate(thermalState(phase)%aTolState (NofMyPhase)) + allocate(thermalState(phase)%dotState (sizeDotState,NofMyPhase)) + allocate(thermalState(phase)%dotState_backup(sizeDotState,NofMyPhase)) + if (any(numerics_integrator == 1_pInt)) then + allocate(thermalState(phase)%previousDotState (sizeDotState,NofMyPhase)) + allocate(thermalState(phase)%previousDotState2 (sizeDotState,NofMyPhase)) + endif + if (any(numerics_integrator == 4_pInt)) & + allocate(thermalState(phase)%RK4dotState (sizeDotState,NofMyPhase)) + if (any(numerics_integrator == 5_pInt)) & + allocate(thermalState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase)) + endif + enddo initializeInstances +#else + allocate(thermal_none_sizeDotState(maxNinstance), source=1_pInt) + allocate(thermal_none_sizeState(maxNinstance), source=1_pInt) +#endif + allocate(thermal_none_sizePostResults(maxNinstance), source=0_pInt) + +end subroutine thermal_none_init + +end module thermal_none