diff --git a/cmake/Compiler-Intel.cmake b/cmake/Compiler-Intel.cmake index 998f60326..60ed46cbc 100644 --- a/cmake/Compiler-Intel.cmake +++ b/cmake/Compiler-Intel.cmake @@ -32,6 +32,8 @@ # disables warnings ... set (COMPILE_FLAGS "${COMPILE_FLAGS} 5268") # ... the text exceeds right hand column allowed on the line (we have only comments there) + set (COMPILE_FLAGS "${COMPILE_FLAGS},7624") + # ... about deprecated forall (has nice syntax and most likely a performance advantage) set (COMPILE_FLAGS "${COMPILE_FLAGS} -warn") # enables warnings ... diff --git a/src/DAMASK_interface.f90 b/src/DAMASK_interface.f90 index b76813fe6..cb13bfaea 100644 --- a/src/DAMASK_interface.f90 +++ b/src/DAMASK_interface.f90 @@ -14,7 +14,11 @@ #define PETSC_MAJOR 3 #define PETSC_MINOR_MIN 10 #define PETSC_MINOR_MAX 11 + module DAMASK_interface + use, intrinsic :: iso_fortran_env + use PETScSys + use prec use system_routines @@ -50,9 +54,6 @@ contains !! information on computation to screen !-------------------------------------------------------------------------------------------------- subroutine DAMASK_interface_init - use, intrinsic :: iso_fortran_env - use PETScSys - #include #if defined(__GFORTRAN__) && __GNUC__ Modeling and Simulations in Materials Science and Engineering 22, 075013 (2014). !-------------------------------------------------------------------------- module Lambert + use prec use math implicit none private - real(pReal), parameter, private :: & + + real(pReal), parameter :: & SPI = sqrt(PI), & PREF = sqrt(6.0_pReal/PI), & A = PI**(5.0_pReal/6.0_pReal)/6.0_pReal**(1.0_pReal/6.0_pReal), & @@ -55,10 +57,8 @@ module Lambert PREK = R1 * 2.0_pReal**(1.0_pReal/4.0_pReal)/BETA public :: & - LambertCubeToBall, & - LambertBallToCube - private :: & - GetPyramidOrder + Lambert_CubeToBall, & + Lambert_BallToCube contains @@ -68,7 +68,7 @@ contains !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief map from 3D cubic grid to 3D ball !-------------------------------------------------------------------------- -function LambertCubeToBall(cube) result(ball) +function Lambert_CubeToBall(cube) result(ball) real(pReal), intent(in), dimension(3) :: cube real(pReal), dimension(3) :: ball, LamXYZ, XYZ @@ -116,7 +116,7 @@ function LambertCubeToBall(cube) result(ball) endif center -end function LambertCubeToBall +end function Lambert_CubeToBall !-------------------------------------------------------------------------- @@ -124,7 +124,7 @@ end function LambertCubeToBall !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief map from 3D ball to 3D cubic grid !-------------------------------------------------------------------------- -pure function LambertBallToCube(xyz) result(cube) +pure function Lambert_BallToCube(xyz) result(cube) real(pReal), intent(in), dimension(3) :: xyz real(pReal), dimension(3) :: cube, xyz1, xyz3 @@ -170,7 +170,7 @@ pure function LambertBallToCube(xyz) result(cube) endif center -end function LambertBallToCube +end function Lambert_BallToCube !-------------------------------------------------------------------------- diff --git a/src/config.f90 b/src/config.f90 index 6bc9e9c0b..cd67c4641 100644 --- a/src/config.f90 +++ b/src/config.f90 @@ -7,10 +7,14 @@ !-------------------------------------------------------------------------------------------------- module config use prec + use DAMASK_interface + use IO + use debug use list implicit none private + type(tPartitionedStringList), public, protected, allocatable, dimension(:) :: & config_phase, & config_microstructure, & @@ -18,10 +22,11 @@ module config config_texture, & config_crystallite - type(tPartitionedStringList), public, protected :: & + type(tPartitionedStringList), public, protected :: & config_numerics, & config_debug + !ToDo: bad names (how should one know that those variables are defined in config?) character(len=64), dimension(:), allocatable, public, protected :: & phase_name, & !< name of each phase homogenization_name, & !< name of each homogenization @@ -45,19 +50,9 @@ contains !> @brief reads material.config and stores its content per part !-------------------------------------------------------------------------------------------------- subroutine config_init - use DAMASK_interface, only: & - getSolverJobName - use IO, only: & - IO_read_ASCII, & - IO_error, & - IO_lc, & - IO_getTag - use debug, only: & - debug_level, & - debug_material, & - debug_levelBasic - integer :: myDebug,i + integer :: i + logical :: verbose character(len=pStringLen) :: & line, & @@ -67,7 +62,7 @@ subroutine config_init write(6,'(/,a)') ' <<<+- config init -+>>>' - myDebug = debug_level(debug_material) + verbose = iand(debug_level(debug_material),debug_levelBasic) /= 0 inquire(file=trim(getSolverJobName())//'.materialConfig',exist=fileExists) if(fileExists) then @@ -87,23 +82,23 @@ subroutine config_init case (trim('phase')) call parse_materialConfig(phase_name,config_phase,line,fileContent(i+1:)) - if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Phase parsed'; flush(6) + if (verbose) write(6,'(a)') ' Phase parsed'; flush(6) case (trim('microstructure')) call parse_materialConfig(microstructure_name,config_microstructure,line,fileContent(i+1:)) - if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Microstructure parsed'; flush(6) + if (verbose) write(6,'(a)') ' Microstructure parsed'; flush(6) case (trim('crystallite')) call parse_materialConfig(crystallite_name,config_crystallite,line,fileContent(i+1:)) - if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Crystallite parsed'; flush(6) + if (verbose) write(6,'(a)') ' Crystallite parsed'; flush(6) case (trim('homogenization')) call parse_materialConfig(homogenization_name,config_homogenization,line,fileContent(i+1:)) - if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Homogenization parsed'; flush(6) + if (verbose) write(6,'(a)') ' Homogenization parsed'; flush(6) case (trim('texture')) call parse_materialConfig(texture_name,config_texture,line,fileContent(i+1:)) - if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Texture parsed'; flush(6) + if (verbose) write(6,'(a)') ' Texture parsed'; flush(6) end select @@ -141,8 +136,6 @@ contains !! Recursion is triggered by "{path/to/inputfile}" in a line !-------------------------------------------------------------------------------------------------- recursive function read_materialConfig(fileName,cnt) result(fileContent) - use IO, only: & - IO_warning character(len=*), intent(in) :: fileName integer, intent(in), optional :: cnt !< recursion counter @@ -226,9 +219,6 @@ end function read_materialConfig subroutine parse_materialConfig(sectionNames,part,line, & fileContent) - use IO, only: & - IO_intOut - character(len=64), allocatable, dimension(:), intent(out) :: sectionNames type(tPartitionedStringList), allocatable, dimension(:), intent(inout) :: part character(len=pStringLen), intent(inout) :: line @@ -298,8 +288,6 @@ end subroutine config_init !> @brief deallocates the linked lists that store the content of the configuration files !-------------------------------------------------------------------------------------------------- subroutine config_deallocate(what) - use IO, only: & - IO_error character(len=*), intent(in) :: what diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 3116345b6..31c859e30 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -9,36 +9,43 @@ !-------------------------------------------------------------------------------------------------- module crystallite - use prec, only: & - pReal, & - pStringLen - use rotations, only: & - rotation - use FEsolving, only: & - FEsolving_execElem, & - FEsolving_execIP - use material, only: & - homogenization_Ngrains + use prec + use IO + use config + use debug + use numerics + use rotations + use math + use mesh + use FEsolving + use material + use constitutive + use lattice use future + use plastic_nonlocal +#if defined(PETSc) || defined(DAMASK_HDF5) + use HDF5_utilities + use results +#endif implicit none private - character(len=64), dimension(:,:), allocatable, private :: & + character(len=64), dimension(:,:), allocatable :: & crystallite_output !< name of each post result output integer, public, protected :: & crystallite_maxSizePostResults !< description not available integer, dimension(:), allocatable, public, protected :: & crystallite_sizePostResults !< description not available - integer, dimension(:,:), allocatable, private :: & + integer, dimension(:,:), allocatable :: & crystallite_sizePostResult !< description not available real(pReal), dimension(:,:,:), allocatable, public :: & crystallite_dt !< requested time increment of each grain - real(pReal), dimension(:,:,:), allocatable, private :: & + real(pReal), dimension(:,:,:), allocatable :: & crystallite_subdt, & !< substepped time increment of each grain crystallite_subFrac, & !< already calculated fraction of increment crystallite_subStep !< size of next integration step - type(rotation), dimension(:,:,:), allocatable, private :: & + type(rotation), dimension(:,:,:), allocatable :: & crystallite_orientation, & !< orientation crystallite_orientation0 !< initial orientation real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: & @@ -63,7 +70,7 @@ module crystallite crystallite_Li, & !< current intermediate velocitiy grad (end of converged time step) crystallite_Li0, & !< intermediate velocitiy grad at start of FE inc crystallite_partionedLi0 !< intermediate velocity grad at start of homog inc - real(pReal), dimension(:,:,:,:,:), allocatable, private :: & + real(pReal), dimension(:,:,:,:,:), allocatable :: & crystallite_subS0, & !< 2nd Piola-Kirchhoff stress vector at start of crystallite inc crystallite_invFp, & !< inverse of current plastic def grad (end of converged time step) crystallite_subFp0,& !< plastic def grad at start of crystallite inc @@ -77,7 +84,7 @@ module crystallite crystallite_dPdF !< current individual dPdF per grain (end of converged time step) logical, dimension(:,:,:), allocatable, public :: & crystallite_requested !< used by upper level (homogenization) to request crystallite calculation - logical, dimension(:,:,:), allocatable, private :: & + logical, dimension(:,:,:), allocatable :: & crystallite_converged, & !< convergence flag crystallite_todo, & !< flag to indicate need for further computation crystallite_localPlasticity !< indicates this grain to have purely local constitutive law @@ -101,16 +108,16 @@ module crystallite neighboringip_ID, & neighboringelement_ID end enum - integer(kind(undefined_ID)),dimension(:,:), allocatable, private :: & + integer(kind(undefined_ID)),dimension(:,:), allocatable :: & crystallite_outputID !< ID of each post result output - type, private :: tOutput !< new requested output (per phase) + type :: tOutput !< new requested output (per phase) character(len=65536), allocatable, dimension(:) :: & label end type tOutput - type(tOutput), allocatable, dimension(:), private :: output_constituent + type(tOutput), allocatable, dimension(:) :: output_constituent - type, private :: tNumerics + type :: tNumerics integer :: & iJacoLpresiduum, & !< frequency of Jacobian update of residuum in Lp nState, & !< state loop limit @@ -138,15 +145,6 @@ module crystallite crystallite_push33ToRef, & crystallite_postResults, & crystallite_results - private :: & - integrateStress, & - integrateState, & - integrateStateFPI, & - integrateStateEuler, & - integrateStateAdaptiveEuler, & - integrateStateRK4, & - integrateStateRKCK45, & - stateJump contains @@ -155,39 +153,6 @@ contains !> @brief allocates and initialize per grain variables !-------------------------------------------------------------------------------------------------- subroutine crystallite_init -#ifdef DEBUG - use debug, only: & - debug_info, & - debug_reset, & - debug_level, & - debug_crystallite, & - debug_levelBasic -#endif - use numerics, only: & - numerics_integrator, & - worldrank, & - usePingPong - use math, only: & - math_I3, & - math_EulerToR, & - math_inv33 - use mesh, only: & - theMesh, & - mesh_element - use IO, only: & - IO_stringValue, & - IO_write_jobFile, & - IO_error - use material - use config, only: & - config_deallocate, & - config_crystallite, & - config_numerics, & - config_phase, & - crystallite_name - use constitutive, only: & - constitutive_initialFi, & - constitutive_microstructure ! derived (shortcut) quantities of given state integer, parameter :: FILEUNIT=434 logical, dimension(:,:), allocatable :: devNull @@ -478,34 +443,6 @@ end subroutine crystallite_init !> @brief calculate stress (P) !-------------------------------------------------------------------------------------------------- function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) - use prec, only: & - tol_math_check, & - dNeq0 -#ifdef DEBUG - use debug, only: & - debug_level, & - debug_crystallite, & - debug_levelBasic, & - debug_levelExtensive, & - debug_levelSelective, & - debug_e, & - debug_i, & - debug_g -#endif - use IO, only: & - IO_warning, & - IO_error - use math, only: & - math_inv33 - use mesh, only: & - theMesh, & - mesh_element - use material, only: & - homogenization_Ngrains, & - plasticState, & - sourceState, & - phase_Nsources, & - phaseAt, phasememberAt logical, dimension(theMesh%elem%nIPs,theMesh%Nelems) :: crystallite_stress real(pReal), intent(in), optional :: & @@ -746,30 +683,6 @@ end function crystallite_stress !> @brief calculate tangent (dPdF) !-------------------------------------------------------------------------------------------------- subroutine crystallite_stressTangent - use prec, only: & - tol_math_check, & - dNeq0 - use IO, only: & - IO_warning, & - IO_error - use math, only: & - math_inv33, & - math_identity2nd, & - math_3333to99, & - math_99to3333, & - math_I3, & - math_mul3333xx3333, & - math_mul33xx33, & - math_invert2, & - math_det33 - use mesh, only: & - mesh_element - use material, only: & - homogenization_Ngrains - use constitutive, only: & - constitutive_SandItsTangents, & - constitutive_LpAndItsTangents, & - constitutive_LiAndItsTangents integer :: & c, & !< counter in integration point component loop @@ -910,19 +823,6 @@ end subroutine crystallite_stressTangent !> @brief calculates orientations !-------------------------------------------------------------------------------------------------- subroutine crystallite_orientations - use math, only: & - math_rotationalPart33, & - math_RtoQ - use material, only: & - plasticState, & - material_phase, & - homogenization_Ngrains - use mesh, only: & - mesh_element - use lattice, only: & - lattice_qDisorientation - use plastic_nonlocal, only: & - plastic_nonlocal_updateCompatibility integer & c, & !< counter in integration point component loop @@ -979,28 +879,6 @@ end function crystallite_push33ToRef !> @brief return results of particular grain !-------------------------------------------------------------------------------------------------- function crystallite_postResults(ipc, ip, el) - use math, only: & - math_det33, & - math_I3, & - inDeg - use mesh, only: & - theMesh, & - mesh_element, & - mesh_ipVolume, & - mesh_ipNeighborhood - use material, only: & - plasticState, & - sourceState, & - microstructure_crystallite, & - crystallite_Noutput, & - material_phase, & - material_texture, & - homogenization_Ngrains - use constitutive, only: & - constitutive_homogenizedC, & - constitutive_postResults - use rotations, only: & - rotation integer, intent(in):: & el, & !< element index @@ -1118,16 +996,9 @@ end function crystallite_postResults !-------------------------------------------------------------------------------------------------- subroutine crystallite_results #if defined(PETSc) || defined(DAMASK_HDF5) - use lattice - use results - use HDF5_utilities - use rotations use config, only: & config_name_phase => phase_name ! anticipate logical name - use material, only: & - material_phase_plasticity_type => phase_plasticity - integer :: p,o real(pReal), allocatable, dimension(:,:,:) :: selected_tensors type(rotation), allocatable, dimension(:) :: selected_rotations @@ -1267,33 +1138,6 @@ end subroutine crystallite_results !> intermediate acceleration of the Newton-Raphson correction !-------------------------------------------------------------------------------------------------- logical function integrateStress(ipc,ip,el,timeFraction) - use, intrinsic :: & - IEEE_arithmetic - use prec, only: tol_math_check, & - dEq0 -#ifdef DEBUG - use debug, only: debug_level, & - debug_e, & - debug_i, & - debug_g, & - debug_crystallite, & - debug_levelBasic, & - debug_levelExtensive, & - debug_levelSelective -#endif - - use constitutive, only: constitutive_LpAndItsTangents, & - constitutive_LiAndItsTangents, & - constitutive_SandItsTangents - use math, only: math_mul33xx33, & - math_mul3333xx3333, & - math_inv33, & - math_det33, & - math_I3, & - math_identity2nd, & - math_3333to99, & - math_33to9, & - math_9to33 integer, intent(in):: el, & ! element index ip, & ! integration point index @@ -1693,27 +1537,6 @@ end function integrateStress !> using Fixed Point Iteration to adapt the stepsize !-------------------------------------------------------------------------------------------------- subroutine integrateStateFPI -#ifdef DEBUG - use debug, only: debug_level, & - debug_e, & - debug_i, & - debug_g, & - debug_crystallite, & - debug_levelBasic, & - debug_levelExtensive, & - debug_levelSelective -#endif - use mesh, only: & - mesh_element - use material, only: & - plasticState, & - sourceState, & - phaseAt, phasememberAt, & - phase_Nsources, & - homogenization_Ngrains - use constitutive, only: & - constitutive_plasticity_maxSizeDotState, & - constitutive_source_maxSizeDotState integer :: & NiterationState, & !< number of iterations in state loop @@ -1901,8 +1724,6 @@ end subroutine integrateStateFPI !> @brief integrate state with 1st order explicit Euler method !-------------------------------------------------------------------------------------------------- subroutine integrateStateEuler - use material, only: & - plasticState call update_dotState(1.0_pReal) call update_state(1.0_pReal) @@ -1919,19 +1740,6 @@ end subroutine integrateStateEuler !> @brief integrate stress, state with 1st order Euler method with adaptive step size !-------------------------------------------------------------------------------------------------- subroutine integrateStateAdaptiveEuler - use mesh, only: & - theMesh, & - mesh_element - use material, only: & - homogenization_Ngrains, & - plasticState, & - sourceState, & - phaseAt, phasememberAt, & - phase_Nsources, & - homogenization_maxNgrains - use constitutive, only: & - constitutive_plasticity_maxSizeDotState, & - constitutive_source_maxSizeDotState integer :: & e, & ! element index in element loop @@ -2025,14 +1833,6 @@ end subroutine integrateStateAdaptiveEuler ! ToDo: This is totally BROKEN: RK4dotState is never used!!! !-------------------------------------------------------------------------------------------------- subroutine integrateStateRK4 - use mesh, only: & - mesh_element - use material, only: & - homogenization_Ngrains, & - plasticState, & - sourceState, & - phase_Nsources, & - phaseAt, phasememberAt real(pReal), dimension(4), parameter :: & TIMESTEPFRACTION = [0.5_pReal, 0.5_pReal, 1.0_pReal, 1.0_pReal] ! factor giving the fraction of the original timestep used for Runge Kutta Integration @@ -2092,19 +1892,6 @@ end subroutine integrateStateRK4 !> adaptive step size (use 5th order solution to advance = "local extrapolation") !-------------------------------------------------------------------------------------------------- subroutine integrateStateRKCK45 - use mesh, only: & - mesh_element, & - theMesh - use material, only: & - homogenization_Ngrains, & - plasticState, & - sourceState, & - phase_Nsources, & - phaseAt, phasememberAt, & - homogenization_maxNgrains - use constitutive, only: & - constitutive_plasticity_maxSizeDotState, & - constitutive_source_maxSizeDotState real(pReal), dimension(5,5), parameter :: & A = reshape([& @@ -2287,8 +2074,6 @@ end subroutine nonlocalConvergenceCheck !> @details: For explicitEuler, RK4 and RKCK45, adaptive Euler and FPI have their on criteria !-------------------------------------------------------------------------------------------------- subroutine setConvergenceFlag - use mesh, only: & - mesh_element integer :: & e, & !< element index in element loop @@ -2327,8 +2112,6 @@ end subroutine setConvergenceFlag !> @brief Standard forwarding of state as state = state0 + dotState * (delta t) !-------------------------------------------------------------------------------------------------- subroutine update_stress(timeFraction) - use mesh, only: & - mesh_element real(pReal), intent(in) :: & timeFraction @@ -2360,8 +2143,6 @@ end subroutine update_stress !> @brief tbd !-------------------------------------------------------------------------------------------------- subroutine update_dependentState - use mesh, only: & - mesh_element use constitutive, only: & constitutive_dependentState => constitutive_microstructure @@ -2387,13 +2168,6 @@ end subroutine update_dependentState !> @brief Standard forwarding of state as state = state0 + dotState * (delta t) !-------------------------------------------------------------------------------------------------- subroutine update_state(timeFraction) - use material, only: & - plasticState, & - sourceState, & - phase_Nsources, & - phaseAt, phasememberAt - use mesh, only: & - mesh_element real(pReal), intent(in) :: & timeFraction @@ -2435,17 +2209,6 @@ end subroutine update_state !> if NaN occurs, crystallite_todo is set to FALSE. Any NaN in a nonlocal propagates to all others !-------------------------------------------------------------------------------------------------- subroutine update_dotState(timeFraction) - use, intrinsic :: & - IEEE_arithmetic - use material, only: & - plasticState, & - sourceState, & - phaseAt, phasememberAt, & - phase_Nsources - use mesh, only: & - mesh_element - use constitutive, only: & - constitutive_collectDotState real(pReal), intent(in) :: & timeFraction @@ -2492,19 +2255,7 @@ end subroutine update_DotState subroutine update_deltaState - use, intrinsic :: & - IEEE_arithmetic - use prec, only: & - dNeq0 - use mesh, only: & - mesh_element - use material, only: & - plasticState, & - sourceState, & - phase_Nsources, & - phaseAt, phasememberAt - use constitutive, only: & - constitutive_collectDeltaState + integer :: & e, & !< element index in element loop i, & !< integration point index in ip loop @@ -2569,29 +2320,6 @@ end subroutine update_deltaState !> returns true, if state jump was successfull or not needed. false indicates NaN in delta state !-------------------------------------------------------------------------------------------------- logical function stateJump(ipc,ip,el) - use, intrinsic :: & - IEEE_arithmetic - use prec, only: & - dNeq0 -#ifdef DEBUG - use debug, only: & - debug_e, & - debug_i, & - debug_g, & - debug_level, & - debug_crystallite, & - debug_levelExtensive, & - debug_levelSelective -#endif - use material, only: & - plasticState, & - sourceState, & - phase_Nsources, & - phaseAt, phasememberAt - use mesh, only: & - mesh_element - use constitutive, only: & - constitutive_collectDeltaState integer, intent(in):: & el, & ! element index diff --git a/src/damage_local.f90 b/src/damage_local.f90 index 1ec42f863..bd71ae95b 100644 --- a/src/damage_local.f90 +++ b/src/damage_local.f90 @@ -4,9 +4,13 @@ !-------------------------------------------------------------------------------------------------- module damage_local use prec + use material + use numerics + use config implicit none private + integer, dimension(:,:), allocatable, target, public :: & damage_local_sizePostResult !< size of each post result output @@ -20,23 +24,22 @@ module damage_local enumerator :: undefined_ID, & damage_ID end enum - integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & + integer(kind(undefined_ID)), dimension(:,:), allocatable :: & damage_local_outputID !< ID of each post result output - type, private :: tParameters + type :: tParameters integer(kind(undefined_ID)), dimension(:), allocatable :: & outputID end type tParameters - type(tparameters), dimension(:), allocatable, private :: & + type(tparameters), dimension(:), allocatable :: & param public :: & damage_local_init, & damage_local_updateState, & damage_local_postResults - private :: & - damage_local_getSourceAndItsTangent + contains @@ -45,23 +48,8 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine damage_local_init - use material, only: & - damage_type, & - damage_typeInstance, & - homogenization_Noutput, & - DAMAGE_local_label, & - DAMAGE_local_ID, & - material_homogenizationAt, & - mappingHomogenization, & - damageState, & - damageMapping, & - damage, & - damage_initialPhi - use config, only: & - config_homogenization - - integer :: maxNinstance,homog,instance,o,i + integer :: maxNinstance,homog,instance,i integer :: sizeState integer :: NofMyHomog, h integer(kind(undefined_ID)) :: & @@ -72,7 +60,7 @@ subroutine damage_local_init write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_local_label//' init -+>>>' - maxNinstance = int(count(damage_type == DAMAGE_local_ID),pInt) + maxNinstance = count(damage_type == DAMAGE_local_ID) if (maxNinstance == 0) return allocate(damage_local_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0) @@ -135,14 +123,6 @@ end subroutine damage_local_init !> @brief calculates local change in damage field !-------------------------------------------------------------------------------------------------- function damage_local_updateState(subdt, ip, el) - use numerics, only: & - residualStiffness, & - err_damage_tolAbs, & - err_damage_tolRel - use material, only: & - material_homogenizationAt, & - mappingHomogenization, & - damageState integer, intent(in) :: & ip, & !< integration point number @@ -177,17 +157,6 @@ end function damage_local_updateState !> @brief calculates homogenized local damage driving forces !-------------------------------------------------------------------------------------------------- subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el) - use material, only: & - homogenization_Ngrains, & - material_homogenizationAt, & - phaseAt, & - phasememberAt, & - phase_source, & - phase_Nsources, & - SOURCE_damage_isoBrittle_ID, & - SOURCE_damage_isoDuctile_ID, & - SOURCE_damage_anisoBrittle_ID, & - SOURCE_damage_anisoDuctile_ID use source_damage_isoBrittle, only: & source_damage_isobrittle_getRateAndItsTangent use source_damage_isoDuctile, only: & @@ -244,15 +213,11 @@ subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el end subroutine damage_local_getSourceAndItsTangent + !-------------------------------------------------------------------------------------------------- !> @brief return array of damage results !-------------------------------------------------------------------------------------------------- function damage_local_postResults(ip,el) - use material, only: & - material_homogenizationAt, & - damage_typeInstance, & - damageMapping, & - damage integer, intent(in) :: & ip, & !< integration point diff --git a/src/damage_none.f90 b/src/damage_none.f90 index aa2995ef5..5ffdba030 100644 --- a/src/damage_none.f90 +++ b/src/damage_none.f90 @@ -3,6 +3,8 @@ !> @brief material subroutine for constant damage field !-------------------------------------------------------------------------------------------------- module damage_none + use config + use material implicit none private @@ -15,18 +17,8 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all neccessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- -subroutine damage_none_init() - use config, only: & - config_homogenization - use material, only: & - damage_initialPhi, & - damage, & - damage_type, & - material_homogenizationAt, & - damageState, & - DAMAGE_NONE_LABEL, & - DAMAGE_NONE_ID - +subroutine damage_none_init + integer :: & homog, & NofMyHomog diff --git a/src/damage_nonlocal.f90 b/src/damage_nonlocal.f90 index 9398b328a..8e61b619b 100644 --- a/src/damage_nonlocal.f90 +++ b/src/damage_nonlocal.f90 @@ -4,39 +4,50 @@ !> @details to be done !-------------------------------------------------------------------------------------------------- module damage_nonlocal - use prec + use prec + use material + use numerics + use config + use crystallite + use lattice + use mesh + use source_damage_isoBrittle + use source_damage_isoDuctile + use source_damage_anisoBrittle + use source_damage_anisoDuctile - implicit none - private - integer, dimension(:,:), allocatable, target, public :: & - damage_nonlocal_sizePostResult !< size of each post result output + implicit none + private + + integer, dimension(:,:), allocatable, target, public :: & + damage_nonlocal_sizePostResult !< size of each post result output - character(len=64), dimension(:,:), allocatable, target, public :: & - damage_nonlocal_output !< name of each post result output - - integer, dimension(:), allocatable, target, public :: & - damage_nonlocal_Noutput !< number of outputs per instance of this damage + character(len=64), dimension(:,:), allocatable, target, public :: & + damage_nonlocal_output !< name of each post result output + + integer, dimension(:), allocatable, target, public :: & + damage_nonlocal_Noutput !< number of outputs per instance of this damage - enum, bind(c) - enumerator :: undefined_ID, & - damage_ID - end enum + enum, bind(c) + enumerator :: undefined_ID, & + damage_ID + end enum - type, private :: tParameters - integer(kind(undefined_ID)), dimension(:), allocatable :: & - outputID - end type tParameters - - type(tparameters), dimension(:), allocatable, private :: & - param + type :: tParameters + integer(kind(undefined_ID)), dimension(:), allocatable :: & + outputID + end type tParameters + + type(tparameters), dimension(:), allocatable :: & + param - public :: & - damage_nonlocal_init, & - damage_nonlocal_getSourceAndItsTangent, & - damage_nonlocal_getDiffusion33, & - damage_nonlocal_getMobility, & - damage_nonlocal_putNonLocalDamage, & - damage_nonlocal_postResults + public :: & + damage_nonlocal_init, & + damage_nonlocal_getSourceAndItsTangent, & + damage_nonlocal_getDiffusion33, & + damage_nonlocal_getMobility, & + damage_nonlocal_putNonLocalDamage, & + damage_nonlocal_postResults contains @@ -45,283 +56,228 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine damage_nonlocal_init - use material, only: & - damage_type, & - damage_typeInstance, & - homogenization_Noutput, & - DAMAGE_nonlocal_label, & - DAMAGE_nonlocal_ID, & - material_homogenizationAt, & - mappingHomogenization, & - damageState, & - damageMapping, & - damage, & - damage_initialPhi - use config, only: & - config_homogenization + integer :: maxNinstance,homog,instance,o,i + integer :: sizeState + integer :: NofMyHomog, h + integer(kind(undefined_ID)) :: & + outputID + character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] + character(len=65536), dimension(:), allocatable :: & + outputs - integer :: maxNinstance,homog,instance,o,i - integer :: sizeState - integer :: NofMyHomog, h - integer(kind(undefined_ID)) :: & - outputID - character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] - character(len=65536), dimension(:), allocatable :: & - outputs - - write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_nonlocal_label//' init -+>>>' - - maxNinstance = int(count(damage_type == DAMAGE_nonlocal_ID)) - if (maxNinstance == 0) return - - allocate(damage_nonlocal_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0) - allocate(damage_nonlocal_output (maxval(homogenization_Noutput),maxNinstance)) - damage_nonlocal_output = '' - allocate(damage_nonlocal_Noutput (maxNinstance), source=0) - - allocate(param(maxNinstance)) + write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_nonlocal_label//' init -+>>>' - do h = 1, size(damage_type) - if (damage_type(h) /= DAMAGE_NONLOCAL_ID) cycle - associate(prm => param(damage_typeInstance(h)), & - config => config_homogenization(h)) - - instance = damage_typeInstance(h) - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) + maxNinstance = count(damage_type == DAMAGE_nonlocal_ID) + if (maxNinstance == 0) return + + allocate(damage_nonlocal_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0) + allocate(damage_nonlocal_output (maxval(homogenization_Noutput),maxNinstance)) + damage_nonlocal_output = '' + allocate(damage_nonlocal_Noutput (maxNinstance), source=0) + + allocate(param(maxNinstance)) - do i=1, size(outputs) - outputID = undefined_ID - select case(outputs(i)) - - case ('damage') - damage_nonlocal_output(i,damage_typeInstance(h)) = outputs(i) - damage_nonlocal_Noutput(instance) = damage_nonlocal_Noutput(instance) + 1 - damage_nonlocal_sizePostResult(i,damage_typeInstance(h)) = 1 - prm%outputID = [prm%outputID , damage_ID] - end select - - enddo + do h = 1, size(damage_type) + if (damage_type(h) /= DAMAGE_NONLOCAL_ID) cycle + associate(prm => param(damage_typeInstance(h)), & + config => config_homogenization(h)) + + instance = damage_typeInstance(h) + outputs = config%getStrings('(output)',defaultVal=emptyStringArray) + allocate(prm%outputID(0)) + + do i=1, size(outputs) + outputID = undefined_ID + select case(outputs(i)) + + case ('damage') + damage_nonlocal_output(i,damage_typeInstance(h)) = outputs(i) + damage_nonlocal_Noutput(instance) = damage_nonlocal_Noutput(instance) + 1 + damage_nonlocal_sizePostResult(i,damage_typeInstance(h)) = 1 + prm%outputID = [prm%outputID , damage_ID] + end select + + enddo - homog = h + homog = h - NofMyHomog = count(material_homogenizationAt == homog) - instance = damage_typeInstance(homog) + NofMyHomog = count(material_homogenizationAt == homog) + instance = damage_typeInstance(homog) -! allocate state arrays - sizeState = 1 - damageState(homog)%sizeState = sizeState - damageState(homog)%sizePostResults = sum(damage_nonlocal_sizePostResult(:,instance)) - allocate(damageState(homog)%state0 (sizeState,NofMyHomog), source=damage_initialPhi(homog)) - allocate(damageState(homog)%subState0(sizeState,NofMyHomog), source=damage_initialPhi(homog)) - allocate(damageState(homog)%state (sizeState,NofMyHomog), source=damage_initialPhi(homog)) +! allocate state arrays + sizeState = 1 + damageState(homog)%sizeState = sizeState + damageState(homog)%sizePostResults = sum(damage_nonlocal_sizePostResult(:,instance)) + allocate(damageState(homog)%state0 (sizeState,NofMyHomog), source=damage_initialPhi(homog)) + allocate(damageState(homog)%subState0(sizeState,NofMyHomog), source=damage_initialPhi(homog)) + allocate(damageState(homog)%state (sizeState,NofMyHomog), source=damage_initialPhi(homog)) - nullify(damageMapping(homog)%p) - damageMapping(homog)%p => mappingHomogenization(1,:,:) - deallocate(damage(homog)%p) - damage(homog)%p => damageState(homog)%state(1,:) - - end associate - enddo + nullify(damageMapping(homog)%p) + damageMapping(homog)%p => mappingHomogenization(1,:,:) + deallocate(damage(homog)%p) + damage(homog)%p => damageState(homog)%state(1,:) + + end associate + enddo end subroutine damage_nonlocal_init + !-------------------------------------------------------------------------------------------------- !> @brief calculates homogenized damage driving forces !-------------------------------------------------------------------------------------------------- subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el) - use material, only: & - homogenization_Ngrains, & - material_homogenizationAt, & - phaseAt, & - phasememberAt, & - phase_source, & - phase_Nsources, & - SOURCE_damage_isoBrittle_ID, & - SOURCE_damage_isoDuctile_ID, & - SOURCE_damage_anisoBrittle_ID, & - SOURCE_damage_anisoDuctile_ID - use source_damage_isoBrittle, only: & - source_damage_isobrittle_getRateAndItsTangent - use source_damage_isoDuctile, only: & - source_damage_isoductile_getRateAndItsTangent - use source_damage_anisoBrittle, only: & - source_damage_anisobrittle_getRateAndItsTangent - use source_damage_anisoDuctile, only: & - source_damage_anisoductile_getRateAndItsTangent - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - phi - integer :: & - phase, & - grain, & - source, & - constituent - real(pReal) :: & - phiDot, dPhiDot_dPhi, localphiDot, dLocalphiDot_dPhi + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + phi + integer :: & + phase, & + grain, & + source, & + constituent + real(pReal) :: & + phiDot, dPhiDot_dPhi, localphiDot, dLocalphiDot_dPhi - phiDot = 0.0_pReal - dPhiDot_dPhi = 0.0_pReal - do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) - phase = phaseAt(grain,ip,el) - constituent = phasememberAt(grain,ip,el) - do source = 1, phase_Nsources(phase) - select case(phase_source(source,phase)) - case (SOURCE_damage_isoBrittle_ID) - call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) + phiDot = 0.0_pReal + dPhiDot_dPhi = 0.0_pReal + do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) + phase = phaseAt(grain,ip,el) + constituent = phasememberAt(grain,ip,el) + do source = 1, phase_Nsources(phase) + select case(phase_source(source,phase)) + case (SOURCE_damage_isoBrittle_ID) + call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - case (SOURCE_damage_isoDuctile_ID) - call source_damage_isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) + case (SOURCE_damage_isoDuctile_ID) + call source_damage_isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - case (SOURCE_damage_anisoBrittle_ID) - call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) + case (SOURCE_damage_anisoBrittle_ID) + call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - case (SOURCE_damage_anisoDuctile_ID) - call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) + case (SOURCE_damage_anisoDuctile_ID) + call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - case default - localphiDot = 0.0_pReal - dLocalphiDot_dPhi = 0.0_pReal + case default + localphiDot = 0.0_pReal + dLocalphiDot_dPhi = 0.0_pReal - end select - phiDot = phiDot + localphiDot - dPhiDot_dPhi = dPhiDot_dPhi + dLocalphiDot_dPhi - enddo - enddo - - phiDot = phiDot/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) - dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) + end select + phiDot = phiDot + localphiDot + dPhiDot_dPhi = dPhiDot_dPhi + dLocalphiDot_dPhi + enddo + enddo + + phiDot = phiDot/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) + dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) end subroutine damage_nonlocal_getSourceAndItsTangent + !-------------------------------------------------------------------------------------------------- !> @brief returns homogenized non local damage diffusion tensor in reference configuration !-------------------------------------------------------------------------------------------------- function damage_nonlocal_getDiffusion33(ip,el) - use numerics, only: & - charLength - use lattice, only: & - lattice_DamageDiffusion33 - use material, only: & - homogenization_Ngrains, & - material_phase, & - material_homogenizationAt - use crystallite, only: & - crystallite_push33ToRef - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), dimension(3,3) :: & - damage_nonlocal_getDiffusion33 - integer :: & - homog, & - grain - - homog = material_homogenizationAt(el) - damage_nonlocal_getDiffusion33 = 0.0_pReal - do grain = 1, homogenization_Ngrains(homog) - damage_nonlocal_getDiffusion33 = damage_nonlocal_getDiffusion33 + & - crystallite_push33ToRef(grain,ip,el,lattice_DamageDiffusion33(1:3,1:3,material_phase(grain,ip,el))) - enddo + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), dimension(3,3) :: & + damage_nonlocal_getDiffusion33 + integer :: & + homog, & + grain + + homog = material_homogenizationAt(el) + damage_nonlocal_getDiffusion33 = 0.0_pReal + do grain = 1, homogenization_Ngrains(homog) + damage_nonlocal_getDiffusion33 = damage_nonlocal_getDiffusion33 + & + crystallite_push33ToRef(grain,ip,el,lattice_DamageDiffusion33(1:3,1:3,material_phase(grain,ip,el))) + enddo - damage_nonlocal_getDiffusion33 = & - charLength**2*damage_nonlocal_getDiffusion33/real(homogenization_Ngrains(homog),pReal) + damage_nonlocal_getDiffusion33 = & + charLength**2*damage_nonlocal_getDiffusion33/real(homogenization_Ngrains(homog),pReal) end function damage_nonlocal_getDiffusion33 + !-------------------------------------------------------------------------------------------------- !> @brief Returns homogenized nonlocal damage mobility !-------------------------------------------------------------------------------------------------- real(pReal) function damage_nonlocal_getMobility(ip,el) - use mesh, only: & - mesh_element - use lattice, only: & - lattice_damageMobility - use material, only: & - material_phase, & - homogenization_Ngrains - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - integer :: & - ipc - - damage_nonlocal_getMobility = 0.0_pReal - - do ipc = 1, homogenization_Ngrains(mesh_element(3,el)) - damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_DamageMobility(material_phase(ipc,ip,el)) - enddo + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + integer :: & + ipc + + damage_nonlocal_getMobility = 0.0_pReal + + do ipc = 1, homogenization_Ngrains(mesh_element(3,el)) + damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_DamageMobility(material_phase(ipc,ip,el)) + enddo - damage_nonlocal_getMobility = damage_nonlocal_getMobility/& - real(homogenization_Ngrains(mesh_element(3,el)),pReal) + damage_nonlocal_getMobility = damage_nonlocal_getMobility/& + real(homogenization_Ngrains(mesh_element(3,el)),pReal) end function damage_nonlocal_getMobility + !-------------------------------------------------------------------------------------------------- !> @brief updated nonlocal damage field with solution from damage phase field PDE !-------------------------------------------------------------------------------------------------- subroutine damage_nonlocal_putNonLocalDamage(phi,ip,el) - use material, only: & - material_homogenizationAt, & - damageMapping, & - damage - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - phi - integer :: & - homog, & - offset - - homog = material_homogenizationAt(el) - offset = damageMapping(homog)%p(ip,el) - damage(homog)%p(offset) = phi + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + phi + integer :: & + homog, & + offset + + homog = material_homogenizationAt(el) + offset = damageMapping(homog)%p(ip,el) + damage(homog)%p(offset) = phi end subroutine damage_nonlocal_putNonLocalDamage - + + !-------------------------------------------------------------------------------------------------- !> @brief return array of damage results !-------------------------------------------------------------------------------------------------- function damage_nonlocal_postResults(ip,el) - use material, only: & - material_homogenizationAt, & - damage_typeInstance, & - damageMapping, & - damage - integer, intent(in) :: & - ip, & !< integration point - el !< element - real(pReal), dimension(sum(damage_nonlocal_sizePostResult(:,damage_typeInstance(material_homogenizationAt(el))))) :: & - damage_nonlocal_postResults + integer, intent(in) :: & + ip, & !< integration point + el !< element + real(pReal), dimension(sum(damage_nonlocal_sizePostResult(:,damage_typeInstance(material_homogenizationAt(el))))) :: & + damage_nonlocal_postResults - integer :: & - instance, homog, offset, o, c - - homog = material_homogenizationAt(el) - offset = damageMapping(homog)%p(ip,el) - instance = damage_typeInstance(homog) - associate(prm => param(instance)) - c = 0 + integer :: & + instance, homog, offset, o, c + + homog = material_homogenizationAt(el) + offset = damageMapping(homog)%p(ip,el) + instance = damage_typeInstance(homog) + associate(prm => param(instance)) + c = 0 - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) - - case (damage_ID) - damage_nonlocal_postResults(c+1) = damage(homog)%p(offset) - c = c + 1 - end select - enddo outputsLoop + outputsLoop: do o = 1,size(prm%outputID) + select case(prm%outputID(o)) + + case (damage_ID) + damage_nonlocal_postResults(c+1) = damage(homog)%p(offset) + c = c + 1 + end select + enddo outputsLoop - end associate + end associate end function damage_nonlocal_postResults end module damage_nonlocal diff --git a/src/debug.f90 b/src/debug.f90 index 4f9566c05..ff084b133 100644 --- a/src/debug.f90 +++ b/src/debug.f90 @@ -6,12 +6,12 @@ !> @brief Reading in and interpretating the debugging settings for the various modules !-------------------------------------------------------------------------------------------------- module debug - use prec, only: & - pInt, & - pReal + use prec + use IO implicit none private + integer(pInt), parameter, public :: & debug_LEVELSELECTIVE = 2_pInt**0_pInt, & debug_LEVELBASIC = 2_pInt**1_pInt, & @@ -78,19 +78,7 @@ contains !> @brief reads in parameters from debug.config and allocates arrays !-------------------------------------------------------------------------------------------------- subroutine debug_init - use prec, only: & - pStringLen - use IO, only: & - IO_read_ASCII, & - IO_error, & - IO_isBlank, & - IO_stringPos, & - IO_stringValue, & - IO_lc, & - IO_floatValue, & - IO_intValue - implicit none character(len=pStringLen), dimension(:), allocatable :: fileContent integer :: i, what, j @@ -253,8 +241,6 @@ end subroutine debug_init !-------------------------------------------------------------------------------------------------- subroutine debug_reset - implicit none - debug_stressMaxLocation = 0_pInt debug_stressMinLocation = 0_pInt debug_jacobianMaxLocation = 0_pInt @@ -272,8 +258,6 @@ end subroutine debug_reset !-------------------------------------------------------------------------------------------------- subroutine debug_info - implicit none - !$OMP CRITICAL (write2out) debugOutputCPFEM: if (iand(debug_level(debug_CPFEM),debug_LEVELBASIC) /= 0 & .and. any(debug_stressMinLocation /= 0_pInt) & diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 3210f02d4..9287cc4bf 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -16,6 +16,12 @@ module homogenization use crystallite use mesh use FEsolving + use thermal_isothermal + use thermal_adiabatic + use thermal_conduction + use damage_none + use damage_local + use damage_nonlocal #if defined(PETSc) || defined(DAMASK_HDF5) use results use HDF5_utilities @@ -131,12 +137,6 @@ contains !> @brief module initialization !-------------------------------------------------------------------------------------------------- subroutine homogenization_init - use thermal_isothermal - use thermal_adiabatic - use thermal_conduction - use damage_none - use damage_local - use damage_nonlocal integer, parameter :: FILEUNIT = 200 integer :: e,i,p @@ -668,10 +668,6 @@ end subroutine partitionDeformation !> "happy" with result !-------------------------------------------------------------------------------------------------- function updateState(ip,el) - use thermal_adiabatic, only: & - thermal_adiabatic_updateState - use damage_local, only: & - damage_local_updateState integer, intent(in) :: & ip, & !< integration point @@ -753,14 +749,6 @@ end subroutine averageStressAndItsTangent !> if homogenization_sizePostResults(i,e) > 0 !! !-------------------------------------------------------------------------------------------------- function postResults(ip,el) - use thermal_adiabatic, only: & - thermal_adiabatic_postResults - use thermal_conduction, only: & - thermal_conduction_postResults - use damage_local, only: & - damage_local_postResults - use damage_nonlocal, only: & - damage_nonlocal_postResults integer, intent(in) :: & ip, & !< integration point diff --git a/src/kinematics_cleavage_opening.f90 b/src/kinematics_cleavage_opening.f90 index 349551d4d..39bfbf340 100644 --- a/src/kinematics_cleavage_opening.f90 +++ b/src/kinematics_cleavage_opening.f90 @@ -5,44 +5,51 @@ !> @details to be done !-------------------------------------------------------------------------------------------------- module kinematics_cleavage_opening - use prec + use prec + use IO + use config + use debug + use math + use lattice + use material - implicit none - private - integer, dimension(:), allocatable, private :: kinematics_cleavage_opening_instance + implicit none + private - type, private :: tParameters !< container type for internal constitutive parameters - integer :: & - totalNcleavage - integer, dimension(:), allocatable :: & - Ncleavage !< active number of cleavage systems per family - real(pReal) :: & - sdot0, & - n - real(pReal), dimension(:), allocatable :: & - critDisp, & - critLoad - end type + integer, dimension(:), allocatable :: kinematics_cleavage_opening_instance + + type :: tParameters !< container type for internal constitutive parameters + integer :: & + totalNcleavage + integer, dimension(:), allocatable :: & + Ncleavage !< active number of cleavage systems per family + real(pReal) :: & + sdot0, & + n + real(pReal), dimension(:), allocatable :: & + critDisp, & + critLoad + end type ! Begin Deprecated - integer, dimension(:), allocatable, private :: & - kinematics_cleavage_opening_totalNcleavage !< total number of cleavage systems - - integer, dimension(:,:), allocatable, private :: & - kinematics_cleavage_opening_Ncleavage !< number of cleavage systems per family - - real(pReal), dimension(:), allocatable, private :: & - kinematics_cleavage_opening_sdot_0, & - kinematics_cleavage_opening_N + integer, dimension(:), allocatable :: & + kinematics_cleavage_opening_totalNcleavage !< total number of cleavage systems + + integer, dimension(:,:), allocatable :: & + kinematics_cleavage_opening_Ncleavage !< number of cleavage systems per family + + real(pReal), dimension(:), allocatable :: & + kinematics_cleavage_opening_sdot_0, & + kinematics_cleavage_opening_N - real(pReal), dimension(:,:), allocatable, private :: & - kinematics_cleavage_opening_critDisp, & - kinematics_cleavage_opening_critLoad + real(pReal), dimension(:,:), allocatable :: & + kinematics_cleavage_opening_critDisp, & + kinematics_cleavage_opening_critLoad ! End Deprecated - public :: & - kinematics_cleavage_opening_init, & - kinematics_cleavage_opening_LiAndItsTangent + public :: & + kinematics_cleavage_opening_init, & + kinematics_cleavage_opening_LiAndItsTangent contains @@ -51,170 +58,144 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine kinematics_cleavage_opening_init() - use debug, only: & - debug_level,& - debug_constitutive,& - debug_levelBasic - use config, only: & - config_phase - use IO, only: & - IO_error - use material, only: & - phase_kinematics, & - KINEMATICS_cleavage_opening_label, & - KINEMATICS_cleavage_opening_ID - use lattice, only: & - lattice_maxNcleavageFamily, & - lattice_NcleavageSystem +subroutine kinematics_cleavage_opening_init - integer, allocatable, dimension(:) :: tempInt - real(pReal), allocatable, dimension(:) :: tempFloat + integer, allocatable, dimension(:) :: tempInt + real(pReal), allocatable, dimension(:) :: tempFloat - integer :: maxNinstance,p,instance,kinematics + integer :: maxNinstance,p,instance - write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_cleavage_opening_LABEL//' init -+>>>' + write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_cleavage_opening_LABEL//' init -+>>>' - maxNinstance = int(count(phase_kinematics == KINEMATICS_cleavage_opening_ID)) - if (maxNinstance == 0) return - - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & - write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance - - allocate(kinematics_cleavage_opening_instance(size(config_phase)), source=0) - do p = 1, size(config_phase) - kinematics_cleavage_opening_instance(p) = count(phase_kinematics(:,1:p) == kinematics_cleavage_opening_ID) ! ToDo: count correct? - enddo - - allocate(kinematics_cleavage_opening_critDisp(lattice_maxNcleavageFamily,maxNinstance), source=0.0_pReal) - allocate(kinematics_cleavage_opening_critLoad(lattice_maxNcleavageFamily,maxNinstance), source=0.0_pReal) - allocate(kinematics_cleavage_opening_Ncleavage(lattice_maxNcleavageFamily,maxNinstance), source=0) - allocate(kinematics_cleavage_opening_totalNcleavage(maxNinstance), source=0) - allocate(kinematics_cleavage_opening_sdot_0(maxNinstance), source=0.0_pReal) - allocate(kinematics_cleavage_opening_N(maxNinstance), source=0.0_pReal) + maxNinstance = count(phase_kinematics == KINEMATICS_cleavage_opening_ID) + if (maxNinstance == 0) return + + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & + write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance + + allocate(kinematics_cleavage_opening_instance(size(config_phase)), source=0) + do p = 1, size(config_phase) + kinematics_cleavage_opening_instance(p) = count(phase_kinematics(:,1:p) == kinematics_cleavage_opening_ID) ! ToDo: count correct? + enddo + + allocate(kinematics_cleavage_opening_critDisp(lattice_maxNcleavageFamily,maxNinstance), source=0.0_pReal) + allocate(kinematics_cleavage_opening_critLoad(lattice_maxNcleavageFamily,maxNinstance), source=0.0_pReal) + allocate(kinematics_cleavage_opening_Ncleavage(lattice_maxNcleavageFamily,maxNinstance), source=0) + allocate(kinematics_cleavage_opening_totalNcleavage(maxNinstance), source=0) + allocate(kinematics_cleavage_opening_sdot_0(maxNinstance), source=0.0_pReal) + allocate(kinematics_cleavage_opening_N(maxNinstance), source=0.0_pReal) - do p = 1, size(config_phase) - if (all(phase_kinematics(:,p) /= KINEMATICS_cleavage_opening_ID)) cycle - instance = kinematics_cleavage_opening_instance(p) - kinematics_cleavage_opening_sdot_0(instance) = config_phase(p)%getFloat('anisobrittle_sdot0') - kinematics_cleavage_opening_N(instance) = config_phase(p)%getFloat('anisobrittle_ratesensitivity') - tempInt = config_phase(p)%getInts('ncleavage') - kinematics_cleavage_opening_Ncleavage(1:size(tempInt),instance) = tempInt + do p = 1, size(config_phase) + if (all(phase_kinematics(:,p) /= KINEMATICS_cleavage_opening_ID)) cycle + instance = kinematics_cleavage_opening_instance(p) + kinematics_cleavage_opening_sdot_0(instance) = config_phase(p)%getFloat('anisobrittle_sdot0') + kinematics_cleavage_opening_N(instance) = config_phase(p)%getFloat('anisobrittle_ratesensitivity') + tempInt = config_phase(p)%getInts('ncleavage') + kinematics_cleavage_opening_Ncleavage(1:size(tempInt),instance) = tempInt - tempFloat = config_phase(p)%getFloats('anisobrittle_criticaldisplacement',requiredSize=size(tempInt)) - kinematics_cleavage_opening_critDisp(1:size(tempInt),instance) = tempFloat + tempFloat = config_phase(p)%getFloats('anisobrittle_criticaldisplacement',requiredSize=size(tempInt)) + kinematics_cleavage_opening_critDisp(1:size(tempInt),instance) = tempFloat - tempFloat = config_phase(p)%getFloats('anisobrittle_criticalload',requiredSize=size(tempInt)) - kinematics_cleavage_opening_critLoad(1:size(tempInt),instance) = tempFloat + tempFloat = config_phase(p)%getFloats('anisobrittle_criticalload',requiredSize=size(tempInt)) + kinematics_cleavage_opening_critLoad(1:size(tempInt),instance) = tempFloat - kinematics_cleavage_opening_Ncleavage(1:lattice_maxNcleavageFamily,instance) = & - min(lattice_NcleavageSystem(1:lattice_maxNcleavageFamily,p),& ! limit active cleavage systems per family to min of available and requested - kinematics_cleavage_opening_Ncleavage(1:lattice_maxNcleavageFamily,instance)) - kinematics_cleavage_opening_totalNcleavage(instance) = sum(kinematics_cleavage_opening_Ncleavage(:,instance)) ! how many cleavage systems altogether - if (kinematics_cleavage_opening_sdot_0(instance) <= 0.0_pReal) & - call IO_error(211,el=instance,ext_msg='sdot_0 ('//KINEMATICS_cleavage_opening_LABEL//')') - if (any(kinematics_cleavage_opening_critDisp(1:size(tempInt),instance) < 0.0_pReal)) & - call IO_error(211,el=instance,ext_msg='critical_displacement ('//KINEMATICS_cleavage_opening_LABEL//')') - if (any(kinematics_cleavage_opening_critLoad(1:size(tempInt),instance) < 0.0_pReal)) & - call IO_error(211,el=instance,ext_msg='critical_load ('//KINEMATICS_cleavage_opening_LABEL//')') - if (kinematics_cleavage_opening_N(instance) <= 0.0_pReal) & - call IO_error(211,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_cleavage_opening_LABEL//')') - enddo + kinematics_cleavage_opening_Ncleavage(1:lattice_maxNcleavageFamily,instance) = & + min(lattice_NcleavageSystem(1:lattice_maxNcleavageFamily,p),& ! limit active cleavage systems per family to min of available and requested + kinematics_cleavage_opening_Ncleavage(1:lattice_maxNcleavageFamily,instance)) + kinematics_cleavage_opening_totalNcleavage(instance) = sum(kinematics_cleavage_opening_Ncleavage(:,instance)) ! how many cleavage systems altogether + if (kinematics_cleavage_opening_sdot_0(instance) <= 0.0_pReal) & + call IO_error(211,el=instance,ext_msg='sdot_0 ('//KINEMATICS_cleavage_opening_LABEL//')') + if (any(kinematics_cleavage_opening_critDisp(1:size(tempInt),instance) < 0.0_pReal)) & + call IO_error(211,el=instance,ext_msg='critical_displacement ('//KINEMATICS_cleavage_opening_LABEL//')') + if (any(kinematics_cleavage_opening_critLoad(1:size(tempInt),instance) < 0.0_pReal)) & + call IO_error(211,el=instance,ext_msg='critical_load ('//KINEMATICS_cleavage_opening_LABEL//')') + if (kinematics_cleavage_opening_N(instance) <= 0.0_pReal) & + call IO_error(211,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_cleavage_opening_LABEL//')') + enddo end subroutine kinematics_cleavage_opening_init - + !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) - use math, only: & - math_mul33xx33 - use material, only: & - material_phase, & - material_homogenizationAt, & - damage, & - damageMapping - use lattice, only: & - lattice_Scleavage, & - lattice_maxNcleavageFamily, & - lattice_NcleavageSystem - integer, intent(in) :: & - ipc, & !< grain number - ip, & !< integration point number - el !< element number - real(pReal), intent(in), dimension(3,3) :: & - S - real(pReal), intent(out), dimension(3,3) :: & - Ld !< damage velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) - integer :: & - instance, phase, & - homog, damageOffset, & - f, i, index_myFamily, k, l, m, n - real(pReal) :: & - traction_d, traction_t, traction_n, traction_crit, & - udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt + integer, intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), intent(in), dimension(3,3) :: & + S + real(pReal), intent(out), dimension(3,3) :: & + Ld !< damage velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) + integer :: & + instance, phase, & + homog, damageOffset, & + f, i, index_myFamily, k, l, m, n + real(pReal) :: & + traction_d, traction_t, traction_n, traction_crit, & + udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt - phase = material_phase(ipc,ip,el) - instance = kinematics_cleavage_opening_instance(phase) - homog = material_homogenizationAt(el) - damageOffset = damageMapping(homog)%p(ip,el) - - Ld = 0.0_pReal - dLd_dTstar = 0.0_pReal - do f = 1,lattice_maxNcleavageFamily - index_myFamily = sum(lattice_NcleavageSystem(1:f-1,phase)) ! at which index starts my family - do i = 1,kinematics_cleavage_opening_Ncleavage(f,instance) ! process each (active) cleavage system in family - traction_d = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase)) - traction_t = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase)) - traction_n = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase)) - traction_crit = kinematics_cleavage_opening_critLoad(f,instance)* & - damage(homog)%p(damageOffset)*damage(homog)%p(damageOffset) - udotd = & - sign(1.0_pReal,traction_d)* & - kinematics_cleavage_opening_sdot_0(instance)* & - (max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**kinematics_cleavage_opening_N(instance) - if (abs(udotd) > tol_math_check) then - Ld = Ld + udotd*lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase) - dudotd_dt = sign(1.0_pReal,traction_d)*udotd*kinematics_cleavage_opening_N(instance)/ & - max(0.0_pReal, abs(traction_d) - traction_crit) - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + & - dudotd_dt*lattice_Scleavage(k,l,1,index_myFamily+i,phase)* & - lattice_Scleavage(m,n,1,index_myFamily+i,phase) - endif + phase = material_phase(ipc,ip,el) + instance = kinematics_cleavage_opening_instance(phase) + homog = material_homogenizationAt(el) + damageOffset = damageMapping(homog)%p(ip,el) + + Ld = 0.0_pReal + dLd_dTstar = 0.0_pReal + do f = 1,lattice_maxNcleavageFamily + index_myFamily = sum(lattice_NcleavageSystem(1:f-1,phase)) ! at which index starts my family + do i = 1,kinematics_cleavage_opening_Ncleavage(f,instance) ! process each (active) cleavage system in family + traction_d = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase)) + traction_t = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase)) + traction_n = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase)) + traction_crit = kinematics_cleavage_opening_critLoad(f,instance)* & + damage(homog)%p(damageOffset)*damage(homog)%p(damageOffset) + udotd = & + sign(1.0_pReal,traction_d)* & + kinematics_cleavage_opening_sdot_0(instance)* & + (max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**kinematics_cleavage_opening_N(instance) + if (abs(udotd) > tol_math_check) then + Ld = Ld + udotd*lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase) + dudotd_dt = sign(1.0_pReal,traction_d)*udotd*kinematics_cleavage_opening_N(instance)/ & + max(0.0_pReal, abs(traction_d) - traction_crit) + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + & + dudotd_dt*lattice_Scleavage(k,l,1,index_myFamily+i,phase)* & + lattice_Scleavage(m,n,1,index_myFamily+i,phase) + endif - udott = & - sign(1.0_pReal,traction_t)* & - kinematics_cleavage_opening_sdot_0(instance)* & - (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**kinematics_cleavage_opening_N(instance) - if (abs(udott) > tol_math_check) then - Ld = Ld + udott*lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase) - dudott_dt = sign(1.0_pReal,traction_t)*udott*kinematics_cleavage_opening_N(instance)/ & - max(0.0_pReal, abs(traction_t) - traction_crit) - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + & - dudott_dt*lattice_Scleavage(k,l,2,index_myFamily+i,phase)* & - lattice_Scleavage(m,n,2,index_myFamily+i,phase) - endif + udott = & + sign(1.0_pReal,traction_t)* & + kinematics_cleavage_opening_sdot_0(instance)* & + (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**kinematics_cleavage_opening_N(instance) + if (abs(udott) > tol_math_check) then + Ld = Ld + udott*lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase) + dudott_dt = sign(1.0_pReal,traction_t)*udott*kinematics_cleavage_opening_N(instance)/ & + max(0.0_pReal, abs(traction_t) - traction_crit) + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + & + dudott_dt*lattice_Scleavage(k,l,2,index_myFamily+i,phase)* & + lattice_Scleavage(m,n,2,index_myFamily+i,phase) + endif - udotn = & - sign(1.0_pReal,traction_n)* & - kinematics_cleavage_opening_sdot_0(instance)* & - (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**kinematics_cleavage_opening_N(instance) - if (abs(udotn) > tol_math_check) then - Ld = Ld + udotn*lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase) - dudotn_dt = sign(1.0_pReal,traction_n)*udotn*kinematics_cleavage_opening_N(instance)/ & - max(0.0_pReal, abs(traction_n) - traction_crit) - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + & - dudotn_dt*lattice_Scleavage(k,l,3,index_myFamily+i,phase)* & - lattice_Scleavage(m,n,3,index_myFamily+i,phase) - endif - enddo - enddo + udotn = & + sign(1.0_pReal,traction_n)* & + kinematics_cleavage_opening_sdot_0(instance)* & + (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**kinematics_cleavage_opening_N(instance) + if (abs(udotn) > tol_math_check) then + Ld = Ld + udotn*lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase) + dudotn_dt = sign(1.0_pReal,traction_n)*udotn*kinematics_cleavage_opening_N(instance)/ & + max(0.0_pReal, abs(traction_n) - traction_crit) + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + & + dudotn_dt*lattice_Scleavage(k,l,3,index_myFamily+i,phase)* & + lattice_Scleavage(m,n,3,index_myFamily+i,phase) + endif + enddo + enddo end subroutine kinematics_cleavage_opening_LiAndItsTangent diff --git a/src/kinematics_slipplane_opening.f90 b/src/kinematics_slipplane_opening.f90 index 7a0b2fe99..3e37e4c0d 100644 --- a/src/kinematics_slipplane_opening.f90 +++ b/src/kinematics_slipplane_opening.f90 @@ -6,12 +6,19 @@ !-------------------------------------------------------------------------------------------------- module kinematics_slipplane_opening use prec + use config + use IO + use debug + use math + use lattice + use material implicit none private - integer, dimension(:), allocatable, private :: kinematics_slipplane_opening_instance + + integer, dimension(:), allocatable :: kinematics_slipplane_opening_instance - type, private :: tParameters !< container type for internal constitutive parameters + type :: tParameters !< container type for internal constitutive parameters integer :: & totalNslip integer, dimension(:), allocatable :: & @@ -19,7 +26,7 @@ module kinematics_slipplane_opening real(pReal) :: & sdot0, & n - real(pReal), dimension(:), allocatable :: & + real(pReal), dimension(:), allocatable :: & critLoad real(pReal), dimension(:,:), allocatable :: & slip_direction, & @@ -27,7 +34,8 @@ module kinematics_slipplane_opening slip_transverse end type tParameters - type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) + type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) + public :: & kinematics_slipplane_opening_init, & kinematics_slipplane_opening_LiAndItsTangent @@ -39,25 +47,9 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine kinematics_slipplane_opening_init() - use debug, only: & - debug_level,& - debug_constitutive,& - debug_levelBasic - use config, only: & - config_phase - use IO, only: & - IO_error - use math, only: & - math_expand - use material, only: & - phase_kinematics, & - KINEMATICS_slipplane_opening_label, & - KINEMATICS_slipplane_opening_ID - use lattice +subroutine kinematics_slipplane_opening_init - - integer :: maxNinstance,p,instance,kinematics + integer :: maxNinstance,p,instance write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_slipplane_opening_LABEL//' init -+>>>' @@ -111,14 +103,6 @@ end subroutine kinematics_slipplane_opening_init !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) - use math, only: & - math_mul33xx33, & - math_outer - use material, only: & - material_phase, & - material_homogenizationAt, & - damage, & - damageMapping integer, intent(in) :: & ipc, & !< grain number diff --git a/src/kinematics_thermal_expansion.f90 b/src/kinematics_thermal_expansion.f90 index 86932ea69..b4f23dfa7 100644 --- a/src/kinematics_thermal_expansion.f90 +++ b/src/kinematics_thermal_expansion.f90 @@ -5,11 +5,17 @@ !-------------------------------------------------------------------------------------------------- module kinematics_thermal_expansion use prec - + use IO + use config + use debug + use math + use lattice + use material + implicit none private - type, private :: tParameters + type :: tParameters real(pReal), allocatable, dimension(:,:,:) :: & expansion end type tParameters @@ -28,19 +34,9 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine kinematics_thermal_expansion_init() - use debug, only: & - debug_level,& - debug_constitutive,& - debug_levelBasic - use material, only: & - phase_kinematics, & - KINEMATICS_thermal_expansion_label, & - KINEMATICS_thermal_expansion_ID - use config, only: & - config_phase +subroutine kinematics_thermal_expansion_init - integer(pInt) :: & + integer :: & Ninstance, & p, i real(pReal), dimension(:), allocatable :: & @@ -48,14 +44,14 @@ subroutine kinematics_thermal_expansion_init() write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_thermal_expansion_LABEL//' init -+>>>' - Ninstance = int(count(phase_kinematics == KINEMATICS_thermal_expansion_ID),pInt) + Ninstance = count(phase_kinematics == KINEMATICS_thermal_expansion_ID) - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance allocate(param(Ninstance)) - do p = 1_pInt, size(phase_kinematics) + do p = 1, size(phase_kinematics) if (all(phase_kinematics(:,p) /= KINEMATICS_thermal_expansion_ID)) cycle ! ToDo: Here we need to decide how to extend the concept of instances to @@ -78,13 +74,8 @@ end subroutine kinematics_thermal_expansion_init !> @brief report initial thermal strain based on current temperature deviation from reference !-------------------------------------------------------------------------------------------------- pure function kinematics_thermal_expansion_initialStrain(homog,phase,offset) - use material, only: & - temperature - use lattice, only: & - lattice_thermalExpansion33, & - lattice_referenceTemperature - integer(pInt), intent(in) :: & + integer, intent(in) :: & phase, & homog, offset real(pReal), dimension(3,3) :: & @@ -106,17 +97,8 @@ end function kinematics_thermal_expansion_initialStrain !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, el) - use material, only: & - material_phase, & - material_homogenizationAt, & - temperature, & - temperatureRate, & - thermalMapping - use lattice, only: & - lattice_thermalExpansion33, & - lattice_referenceTemperature - integer(pInt), intent(in) :: & + integer, intent(in) :: & ipc, & !< grain number ip, & !< integration point number el !< element number @@ -124,7 +106,7 @@ subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, Li !< thermal velocity gradient real(pReal), intent(out), dimension(3,3,3,3) :: & dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero) - integer(pInt) :: & + integer :: & phase, & homog, offset real(pReal) :: & diff --git a/src/lattice.f90 b/src/lattice.f90 index 1a7508984..43fc25530 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -7,8 +7,10 @@ ! and cleavage as well as interaction among the various systems !-------------------------------------------------------------------------------------------------- module lattice - use prec, only: & - pReal + use prec + use IO + use config + use math use future implicit none @@ -28,25 +30,25 @@ module lattice !-------------------------------------------------------------------------------------------------- ! face centered cubic - integer, dimension(2), parameter, private :: & + integer, dimension(2), parameter :: & LATTICE_FCC_NSLIPSYSTEM = [12, 6] !< # of slip systems per family for fcc - integer, dimension(1), parameter, private :: & + integer, dimension(1), parameter :: & LATTICE_FCC_NTWINSYSTEM = [12] !< # of twin systems per family for fcc - integer, dimension(1), parameter, private :: & + integer, dimension(1), parameter :: & LATTICE_FCC_NTRANSSYSTEM = [12] !< # of transformation systems per family for fcc - integer, dimension(2), parameter, private :: & + integer, dimension(2), parameter :: & LATTICE_FCC_NCLEAVAGESYSTEM = [3, 4] !< # of cleavage systems per family for fcc - integer, parameter, private :: & + integer, parameter :: & LATTICE_FCC_NSLIP = sum(LATTICE_FCC_NSLIPSYSTEM), & !< total # of slip systems for fcc LATTICE_FCC_NTWIN = sum(LATTICE_FCC_NTWINSYSTEM), & !< total # of twin systems for fcc LATTICE_FCC_NTRANS = sum(LATTICE_FCC_NTRANSSYSTEM), & !< total # of transformation systems for fcc LATTICE_FCC_NCLEAVAGE = sum(LATTICE_FCC_NCLEAVAGESYSTEM) !< total # of cleavage systems for fcc - real(pReal), dimension(3+3,LATTICE_FCC_NSLIP), parameter, private :: & + real(pReal), dimension(3+3,LATTICE_FCC_NSLIP), parameter :: & LATTICE_FCC_SYSTEMSLIP = reshape(real([& ! Slip direction Plane normal ! SCHMID-BOAS notation 0, 1,-1, 1, 1, 1, & ! B2 @@ -70,11 +72,11 @@ module lattice 0, 1,-1, 0, 1, 1 & ],pReal),shape(LATTICE_FCC_SYSTEMSLIP)) !< Slip system <110>{111} directions. Sorted according to Eisenlohr & Hantcherli - character(len=*), dimension(2), parameter, private :: LATTICE_FCC_SLIPFAMILY_NAME = & + character(len=*), dimension(2), parameter :: LATTICE_FCC_SLIPFAMILY_NAME = & ['<0 1 -1>{1 1 1}', & '<0 1 -1>{0 1 1}'] - real(pReal), dimension(3+3,LATTICE_FCC_NTWIN), parameter, private :: & + real(pReal), dimension(3+3,LATTICE_FCC_NTWIN), parameter :: & LATTICE_FCC_SYSTEMTWIN = reshape(real( [& -2, 1, 1, 1, 1, 1, & 1,-2, 1, 1, 1, 1, & @@ -90,7 +92,7 @@ module lattice -1, 1, 2, -1, 1,-1 & ],pReal),shape(LATTICE_FCC_SYSTEMTWIN)) !< Twin system <112>{111} directions. Sorted according to Eisenlohr & Hantcherli - character(len=*), dimension(1), parameter, private :: LATTICE_FCC_TWINFAMILY_NAME = & + character(len=*), dimension(1), parameter :: LATTICE_FCC_TWINFAMILY_NAME = & ['<-2 1 1>{1 1 1}'] @@ -110,7 +112,7 @@ module lattice 10,11 & ],shape(LATTICE_FCC_TWINNUCLEATIONSLIPPAIR)) - real(pReal), dimension(3+3,LATTICE_FCC_NCLEAVAGE), parameter, private :: & + real(pReal), dimension(3+3,LATTICE_FCC_NCLEAVAGE), parameter :: & LATTICE_FCC_SYSTEMCLEAVAGE = reshape(real([& ! Cleavage direction Plane normal 0, 1, 0, 1, 0, 0, & @@ -124,21 +126,21 @@ module lattice !-------------------------------------------------------------------------------------------------- ! body centered cubic - integer, dimension(2), parameter, private :: & + integer, dimension(2), parameter :: & LATTICE_BCC_NSLIPSYSTEM = [12, 12] !< # of slip systems per family for bcc - integer, dimension(1), parameter, private :: & + integer, dimension(1), parameter :: & LATTICE_BCC_NTWINSYSTEM = [12] !< # of twin systems per family for bcc - integer, dimension(2), parameter, private :: & + integer, dimension(2), parameter :: & LATTICE_BCC_NCLEAVAGESYSTEM = [3, 6] !< # of cleavage systems per family for bcc - integer, parameter, private :: & + integer, parameter :: & LATTICE_BCC_NSLIP = sum(LATTICE_BCC_NSLIPSYSTEM), & !< total # of slip systems for bcc LATTICE_BCC_NTWIN = sum(LATTICE_BCC_NTWINSYSTEM), & !< total # of twin systems for bcc LATTICE_BCC_NCLEAVAGE = sum(LATTICE_BCC_NCLEAVAGESYSTEM) !< total # of cleavage systems for bcc - real(pReal), dimension(3+3,LATTICE_BCC_NSLIP), parameter, private :: & + real(pReal), dimension(3+3,LATTICE_BCC_NSLIP), parameter :: & LATTICE_BCC_SYSTEMSLIP = reshape(real([& ! Slip direction Plane normal ! Slip system <111>{110} @@ -169,11 +171,11 @@ module lattice 1, 1, 1, 1, 1,-2 & ],pReal),shape(LATTICE_BCC_SYSTEMSLIP)) - character(len=*), dimension(2), parameter, private :: LATTICE_BCC_SLIPFAMILY_NAME = & + character(len=*), dimension(2), parameter :: LATTICE_BCC_SLIPFAMILY_NAME = & ['<1 -1 1>{0 1 1}', & '<1 -1 1>{2 1 1}'] - real(pReal), dimension(3+3,LATTICE_BCC_NTWIN), parameter, private :: & + real(pReal), dimension(3+3,LATTICE_BCC_NTWIN), parameter :: & LATTICE_BCC_SYSTEMTWIN = reshape(real([& ! Twin system <111>{112} -1, 1, 1, 2, 1, 1, & @@ -190,10 +192,10 @@ module lattice 1, 1, 1, 1, 1,-2 & ],pReal),shape(LATTICE_BCC_SYSTEMTWIN)) - character(len=*), dimension(1), parameter, private :: LATTICE_BCC_TWINFAMILY_NAME = & + character(len=*), dimension(1), parameter :: LATTICE_BCC_TWINFAMILY_NAME = & ['<1 1 1>{2 1 1}'] - real(pReal), dimension(3+3,LATTICE_BCC_NCLEAVAGE), parameter, private :: & + real(pReal), dimension(3+3,LATTICE_BCC_NCLEAVAGE), parameter :: & LATTICE_BCC_SYSTEMCLEAVAGE = reshape(real([& ! Cleavage direction Plane normal 0, 1, 0, 1, 0, 0, & @@ -209,21 +211,21 @@ module lattice !-------------------------------------------------------------------------------------------------- ! hexagonal - integer, dimension(6), parameter, private :: & + integer, dimension(6), parameter :: & LATTICE_HEX_NSLIPSYSTEM = [3, 3, 3, 6, 12, 6] !< # of slip systems per family for hex - integer, dimension(4), parameter, private :: & + integer, dimension(4), parameter :: & LATTICE_HEX_NTWINSYSTEM = [6, 6, 6, 6] !< # of slip systems per family for hex - integer, dimension(1), parameter, private :: & + integer, dimension(1), parameter :: & LATTICE_HEX_NCLEAVAGESYSTEM = [3] !< # of cleavage systems per family for hex - integer, parameter, private :: & + integer, parameter :: & LATTICE_HEX_NSLIP = sum(LATTICE_HEX_NSLIPSYSTEM), & !< total # of slip systems for hex LATTICE_HEX_NTWIN = sum(LATTICE_HEX_NTWINSYSTEM), & !< total # of twin systems for hex LATTICE_HEX_NCLEAVAGE = sum(LATTICE_HEX_NCLEAVAGESYSTEM) !< total # of cleavage systems for hex - real(pReal), dimension(4+4,LATTICE_HEX_NSLIP), parameter, private :: & + real(pReal), dimension(4+4,LATTICE_HEX_NSLIP), parameter :: & LATTICE_HEX_SYSTEMSLIP = reshape(real([& ! Slip direction Plane normal ! Basal systems <11.0>{00.1} (independent of c/a-ratio, Bravais notation (4 coordinate base)) @@ -267,7 +269,7 @@ module lattice 1, 1, -2, 3, -1, -1, 2, 2 & ],pReal),shape(LATTICE_HEX_SYSTEMSLIP)) !< slip systems for hex sorted by A. Alankar & P. Eisenlohr - character(len=*), dimension(6), parameter, private :: LATTICE_HEX_SLIPFAMILY_NAME = & + character(len=*), dimension(6), parameter :: LATTICE_HEX_SLIPFAMILY_NAME = & ['<1 1 . 1>{0 0 . 1} ', & '<1 1 . 1>{1 0 . 0} ', & '<1 0 . 0>{1 1 . 0} ', & @@ -275,7 +277,7 @@ module lattice '<1 1 . 3>{-1 0 . 1} ', & '<1 1 . 3>{-1 -1 . 2}'] - real(pReal), dimension(4+4,LATTICE_HEX_NTWIN), parameter, private :: & + real(pReal), dimension(4+4,LATTICE_HEX_NTWIN), parameter :: & LATTICE_HEX_SYSTEMTWIN = reshape(real([& ! Compression or Tension =f(twinning shear=f(c/a)) for each metal ! (according to Yoo 1981) 1, -1, 0, 1, -1, 1, 0, 2, & ! <-10.1>{10.2} shear = (3-(c/a)^2)/(sqrt(3) c/a) @@ -307,13 +309,13 @@ module lattice 1, 1, -2, -3, 1, 1, -2, 2 & ],pReal),shape(LATTICE_HEX_SYSTEMTWIN)) !< twin systems for hex, order follows Prof. Tom Bieler's scheme - character(len=*), dimension(4), parameter, private :: LATTICE_HEX_TWINFAMILY_NAME = & + character(len=*), dimension(4), parameter :: LATTICE_HEX_TWINFAMILY_NAME = & ['<-1 0 . 1>{1 0 . 2} ', & '<1 1 . 6>{-1 -1 . 1}', & '<1 0 . -2>{1 0 . 1} ', & '<1 1 . -3>{1 1 . 2} '] - real(pReal), dimension(4+4,LATTICE_HEX_NCLEAVAGE), parameter, private :: & + real(pReal), dimension(4+4,LATTICE_HEX_NCLEAVAGE), parameter :: & LATTICE_HEX_SYSTEMCLEAVAGE = reshape(real([& ! Cleavage direction Plane normal 2,-1,-1, 0, 0, 0, 0, 1, & @@ -324,13 +326,13 @@ module lattice !-------------------------------------------------------------------------------------------------- ! body centered tetragonal - integer, dimension(13), parameter, private :: & + integer, dimension(13), parameter :: & LATTICE_BCT_NSLIPSYSTEM = [2, 2, 2, 4, 2, 4, 2, 2, 4, 8, 4, 8, 8 ] !< # of slip systems per family for bct (Sn) Bieler J. Electr Mater 2009 - integer, parameter, private :: & + integer, parameter :: & LATTICE_BCT_NSLIP = sum(LATTICE_BCT_NSLIPSYSTEM) !< total # of slip systems for bct - real(pReal), dimension(3+3,LATTICE_BCT_NSLIP), parameter, private :: & + real(pReal), dimension(3+3,LATTICE_BCT_NSLIP), parameter :: & LATTICE_BCT_SYSTEMSLIP = reshape(real([& ! Slip direction Plane normal ! Slip family 1 {100)<001] (Bravais notation {hkl) @brief Module initialization !-------------------------------------------------------------------------------------------------- subroutine lattice_init - use IO, only: & - IO_error - use config, only: & - config_phase integer :: Nphases character(len=65536) :: & @@ -654,15 +652,7 @@ end subroutine lattice_init !> @brief !!!!!!!DEPRECTATED!!!!!! !-------------------------------------------------------------------------------------------------- subroutine lattice_initializeStructure(myPhase,CoverA) - use prec, only: & - tol_math_check - use math, only: & - math_sym3333to66, & - math_Voigt66to3333, & - math_cross - use IO, only: & - IO_error - + integer, intent(in) :: myPhase real(pReal), intent(in) :: & CoverA @@ -690,9 +680,10 @@ subroutine lattice_initializeStructure(myPhase,CoverA) call IO_error(135,el=i,ip=myPhase,ext_msg='matrix diagonal "el"ement of phase "ip"') enddo - forall (i = 1:3) & + do i = 1,3 lattice_thermalExpansion33 (1:3,1:3,i,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& lattice_thermalExpansion33 (1:3,1:3,i,myPhase)) + enddo lattice_thermalConductivity33 (1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& lattice_thermalConductivity33 (1:3,1:3,myPhase)) @@ -763,17 +754,17 @@ pure function lattice_symmetrizeC66(struct,C66) select case(struct) case (LATTICE_iso_ID) - forall(k=1:3) + do k=1,3 forall(j=1:3) lattice_symmetrizeC66(k,j) = C66(1,2) lattice_symmetrizeC66(k,k) = C66(1,1) lattice_symmetrizeC66(k+3,k+3) = 0.5_pReal*(C66(1,1)-C66(1,2)) - end forall + enddo case (LATTICE_fcc_ID,LATTICE_bcc_ID) - forall(k=1:3) + do k=1,3 forall(j=1:3) lattice_symmetrizeC66(k,j) = C66(1,2) lattice_symmetrizeC66(k,k) = C66(1,1) lattice_symmetrizeC66(k+3,k+3) = C66(4,4) - end forall + enddo case (LATTICE_hex_ID) lattice_symmetrizeC66(1,1) = C66(1,1) lattice_symmetrizeC66(2,2) = C66(1,1) @@ -834,7 +825,9 @@ pure function lattice_symmetrize33(struct,T33) select case(struct) case (LATTICE_iso_ID,LATTICE_fcc_ID,LATTICE_bcc_ID) - forall(k=1:3) lattice_symmetrize33(k,k) = T33(1,1) + do k=1,3 + lattice_symmetrize33(k,k) = T33(1,1) + enddo case (LATTICE_hex_ID) lattice_symmetrize33(1,1) = T33(1,1) lattice_symmetrize33(2,2) = T33(1,1) @@ -854,10 +847,6 @@ end function lattice_symmetrize33 !> @brief figures whether unit quat falls into stereographic standard triangle !-------------------------------------------------------------------------------------------------- logical pure function lattice_qInSST(Q, struct) - use, intrinsic :: & - IEEE_arithmetic - use math, only: & - math_qToRodrig real(pReal), dimension(4), intent(in) :: Q ! orientation integer(kind(LATTICE_undefined_ID)), intent(in) :: struct ! lattice structure @@ -888,11 +877,6 @@ end function lattice_qInSST !> @brief calculates the disorientation for 2 unit quaternions !-------------------------------------------------------------------------------------------------- pure function lattice_qDisorientation(Q1, Q2, struct) - use prec, only: & - tol_math_check - use math, only: & - math_qMul, & - math_qConj real(pReal), dimension(4) :: lattice_qDisorientation real(pReal), dimension(4), intent(in) :: & @@ -998,8 +982,6 @@ end function lattice_qDisorientation !> @brief Characteristic shear for twinning !-------------------------------------------------------------------------------------------------- function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(characteristicShear) - use IO, only: & - IO_error integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family character(len=*), intent(in) :: structure !< lattice structure @@ -1077,14 +1059,6 @@ end function lattice_characteristicShear_Twin !> @brief Rotated elasticity matrices for twinning in 66-vector notation !-------------------------------------------------------------------------------------------------- function lattice_C66_twin(Ntwin,C66,structure,CoverA) - use IO, only: & - IO_error - use math, only: & - PI, & - math_axisAngleToR, & - math_sym3333to66, & - math_66toSym3333, & - math_rotate_forward3333 integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family character(len=*), intent(in) :: structure !< lattice structure @@ -1125,17 +1099,6 @@ end function lattice_C66_twin !-------------------------------------------------------------------------------------------------- function lattice_C66_trans(Ntrans,C_parent66,structure_target, & CoverA_trans,a_bcc,a_fcc) - use prec, only: & - tol_math_check - use IO, only: & - IO_error - use math, only: & - INRAD, & - MATH_I3, & - math_axisAngleToR, & - math_sym3333to66, & - math_66toSym3333, & - math_rotate_forward3333 integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family character(len=*), intent(in) :: structure_target !< lattice structure @@ -1196,13 +1159,6 @@ function lattice_C66_trans(Ntrans,C_parent66,structure_target, & ! Gröger et al. 2008, Acta Materialia 56 (2008) 5412–5425, table 1 !-------------------------------------------------------------------------------------------------- function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSchmidMatrix) - use IO, only: & - IO_error - use math, only: & - INRAD, & - math_outer, & - math_cross, & - math_axisAngleToR integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family real(pReal), dimension(:), intent(in) :: nonSchmidCoefficients !< non-Schmid coefficients for projections @@ -1246,9 +1202,7 @@ end function lattice_nonSchmidMatrix !> details only active slip systems are considered !-------------------------------------------------------------------------------------------------- function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) result(interactionMatrix) - use IO, only: & - IO_error - + integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family real(pReal), dimension(:), intent(in) :: interactionValues !< values for slip-slip interaction character(len=*), intent(in) :: structure !< lattice structure @@ -1468,8 +1422,6 @@ end function lattice_interaction_SlipBySlip !> details only active twin systems are considered !-------------------------------------------------------------------------------------------------- function lattice_interaction_TwinByTwin(Ntwin,interactionValues,structure) result(interactionMatrix) - use IO, only: & - IO_error integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family real(pReal), dimension(:), intent(in) :: interactionValues !< values for twin-twin interaction @@ -1571,8 +1523,6 @@ end function lattice_interaction_TwinByTwin !> details only active trans systems are considered !-------------------------------------------------------------------------------------------------- function lattice_interaction_TransByTrans(Ntrans,interactionValues,structure) result(interactionMatrix) - use IO, only: & - IO_error integer, dimension(:), intent(in) :: Ntrans !< number of active trans systems per family real(pReal), dimension(:), intent(in) :: interactionValues !< values for trans-trans interaction @@ -1618,8 +1568,6 @@ end function lattice_interaction_TransByTrans !> details only active slip and twin systems are considered !-------------------------------------------------------------------------------------------------- function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,structure) result(interactionMatrix) - use IO, only: & - IO_error integer, dimension(:), intent(in) :: Nslip, & !< number of active slip systems per family Ntwin !< number of active twin systems per family @@ -1760,8 +1708,6 @@ end function lattice_interaction_SlipByTwin !> details only active slip and trans systems are considered !-------------------------------------------------------------------------------------------------- function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,structure) result(interactionMatrix) - use IO, only: & - IO_error integer, dimension(:), intent(in) :: Nslip, & !< number of active slip systems per family Ntrans !< number of active trans systems per family @@ -1818,8 +1764,6 @@ function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,structur !> details only active twin and slip systems are considered !-------------------------------------------------------------------------------------------------- function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,structure) result(interactionMatrix) - use IO, only: & - IO_error integer, dimension(:), intent(in) :: Ntwin, & !< number of active twin systems per family Nslip !< number of active slip systems per family @@ -1898,13 +1842,6 @@ end function lattice_interaction_TwinBySlip !> details only active slip systems are considered !-------------------------------------------------------------------------------------------------- function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) result(SchmidMatrix) - use prec, only: & - tol_math_check - use IO, only: & - IO_error - use math, only: & - math_trace33, & - math_outer integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family character(len=*), intent(in) :: structure !< lattice structure @@ -1957,13 +1894,6 @@ end function lattice_SchmidMatrix_slip !> details only active twin systems are considered !-------------------------------------------------------------------------------------------------- function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) result(SchmidMatrix) - use prec, only: & - tol_math_check - use IO, only: & - IO_error - use math, only: & - math_trace33, & - math_outer integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family character(len=*), intent(in) :: structure !< lattice structure @@ -2013,8 +1943,6 @@ function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) result(SchmidMatrix) !> details only active twin systems are considered !-------------------------------------------------------------------------------------------------- function lattice_SchmidMatrix_trans(Ntrans,structure_target,cOverA,a_bcc,a_fcc) result(SchmidMatrix) - use IO, only: & - IO_error integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family real(pReal), intent(in) :: cOverA !< c/a ratio @@ -2041,11 +1969,7 @@ end function lattice_SchmidMatrix_trans !> details only active cleavage systems are considered !-------------------------------------------------------------------------------------------------- function lattice_SchmidMatrix_cleavage(Ncleavage,structure,cOverA) result(SchmidMatrix) - use math, only: & - math_outer - use IO, only: & - IO_error - + integer, dimension(:), intent(in) :: Ncleavage !< number of active cleavage systems per family character(len=*), intent(in) :: structure !< lattice structure real(pReal), intent(in) :: cOverA !< c/a ratio @@ -2154,8 +2078,6 @@ end function lattice_slip_transverse !> @details: This projection is used to calculate forest hardening for edge dislocations !-------------------------------------------------------------------------------------------------- function slipProjection_transverse(Nslip,structure,cOverA) result(projection) - use math, only: & - math_inner integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family character(len=*), intent(in) :: structure !< lattice structure @@ -2179,8 +2101,6 @@ end function slipProjection_transverse !> @details: This projection is used to calculate forest hardening for screw dislocations !-------------------------------------------------------------------------------------------------- function slipProjection_direction(Nslip,structure,cOverA) result(projection) - use math, only: & - math_inner integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family character(len=*), intent(in) :: structure !< lattice structure @@ -2204,9 +2124,7 @@ end function slipProjection_direction !> @details Order: Direction, plane (normal), and common perpendicular !-------------------------------------------------------------------------------------------------- function coordinateSystem_slip(Nslip,structure,cOverA) result(coordinateSystem) - use IO, only: & - IO_error - + integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family character(len=*), intent(in) :: structure !< lattice structure real(pReal), intent(in) :: cOverA !< c/a ratio @@ -2249,8 +2167,6 @@ end function coordinateSystem_slip !> @brief Populates reduced interaction matrix !-------------------------------------------------------------------------------------------------- function buildInteraction(reacting_used,acting_used,reacting_max,acting_max,values,matrix) - use IO, only: & - IO_error integer, dimension(:), intent(in) :: & reacting_used, & !< # of reacting systems per family as specified in material.config @@ -2295,10 +2211,6 @@ end function buildInteraction !> @details Order: Direction, plane (normal), and common perpendicular !-------------------------------------------------------------------------------------------------- function buildCoordinateSystem(active,complete,system,structure,cOverA) - use IO, only: & - IO_error - use math, only: & - math_cross integer, dimension(:), intent(in) :: & active, & @@ -2370,16 +2282,6 @@ end function buildCoordinateSystem ! set a_bcc = 0.0 for fcc -> hex transformation !-------------------------------------------------------------------------------------------------- subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc) - use prec, only: & - dEq0 - use math, only: & - math_cross, & - math_outer, & - math_axisAngleToR, & - INRAD, & - MATH_I3 - use IO, only: & - IO_error integer, dimension(:), intent(in) :: & Ntrans diff --git a/src/list.f90 b/src/list.f90 index be80b151d..79eafc964 100644 --- a/src/list.f90 +++ b/src/list.f90 @@ -3,8 +3,8 @@ !> @brief linked list !-------------------------------------------------------------------------------------------------- module list - use prec, only: & - pReal + use prec + use IO implicit none private @@ -65,10 +65,6 @@ contains !! to lower case. The data is not stored in the new element but in the current. !-------------------------------------------------------------------------------------------------- subroutine add(this,string) - use IO, only: & - IO_isBlank, & - IO_lc, & - IO_stringPos class(tPartitionedStringList), target, intent(in) :: this character(len=*), intent(in) :: string @@ -157,8 +153,6 @@ end subroutine finalizeArray !> @brief reports wether a given key (string value at first position) exists in the list !-------------------------------------------------------------------------------------------------- logical function keyExists(this,key) - use IO, only: & - IO_stringValue class(tPartitionedStringList), target, intent(in) :: this character(len=*), intent(in) :: key @@ -180,8 +174,6 @@ end function keyExists !> @details traverses list and counts each occurrence of specified key !-------------------------------------------------------------------------------------------------- integer function countKeys(this,key) - use IO, only: & - IO_stringValue class(tPartitionedStringList), target, intent(in) :: this character(len=*), intent(in) :: key @@ -205,10 +197,6 @@ end function countKeys !! error unless default is given !-------------------------------------------------------------------------------------------------- real(pReal) function getFloat(this,key,defaultVal) - use IO, only : & - IO_error, & - IO_stringValue, & - IO_FloatValue class(tPartitionedStringList), target, intent(in) :: this character(len=*), intent(in) :: key @@ -241,10 +229,6 @@ end function getFloat !! error unless default is given !-------------------------------------------------------------------------------------------------- integer function getInt(this,key,defaultVal) - use IO, only: & - IO_error, & - IO_stringValue, & - IO_IntValue class(tPartitionedStringList), target, intent(in) :: this character(len=*), intent(in) :: key @@ -278,9 +262,6 @@ end function getInt !! the individual chunks are returned !-------------------------------------------------------------------------------------------------- character(len=65536) function getString(this,key,defaultVal,raw) - use IO, only: & - IO_error, & - IO_stringValue class(tPartitionedStringList), target, intent(in) :: this character(len=*), intent(in) :: key @@ -327,10 +308,6 @@ end function getString !! values from the last occurrence. If key is not found exits with error unless default is given. !-------------------------------------------------------------------------------------------------- function getFloats(this,key,defaultVal,requiredSize) - use IO, only: & - IO_error, & - IO_stringValue, & - IO_FloatValue real(pReal), dimension(:), allocatable :: getFloats class(tPartitionedStringList), target, intent(in) :: this @@ -376,10 +353,6 @@ end function getFloats !! values from the last occurrence. If key is not found exits with error unless default is given. !-------------------------------------------------------------------------------------------------- function getInts(this,key,defaultVal,requiredSize) - use IO, only: & - IO_error, & - IO_stringValue, & - IO_IntValue integer, dimension(:), allocatable :: getInts class(tPartitionedStringList), target, intent(in) :: this @@ -426,9 +399,6 @@ end function getInts !! If raw is true, the the complete string is returned, otherwise the individual chunks are returned !-------------------------------------------------------------------------------------------------- function getStrings(this,key,defaultVal,raw) - use IO, only: & - IO_error, & - IO_StringValue character(len=65536),dimension(:), allocatable :: getStrings class(tPartitionedStringList),target, intent(in) :: this diff --git a/src/math.f90 b/src/math.f90 index 1740ebdb7..4a32be274 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -10,12 +10,20 @@ module math use future implicit none - real(pReal), parameter, public :: PI = acos(-1.0_pReal) !< ratio of a circle's circumference to its diameter - real(pReal), parameter, public :: INDEG = 180.0_pReal/PI !< conversion from radian into degree - real(pReal), parameter, public :: INRAD = PI/180.0_pReal !< conversion from degree into radian - complex(pReal), parameter, public :: TWOPIIMG = cmplx(0.0_pReal,2.0_pReal*PI) !< Re(0.0), Im(2xPi) + public +#if __INTEL_COMPILER >= 1900 + ! do not make use associated entities available to other modules + private :: & + prec, & + future +#endif - real(pReal), dimension(3,3), parameter, public :: & + real(pReal), parameter :: PI = acos(-1.0_pReal) !< ratio of a circle's circumference to its diameter + real(pReal), parameter :: INDEG = 180.0_pReal/PI !< conversion from radian into degree + real(pReal), parameter :: INRAD = PI/180.0_pReal !< conversion from degree into radian + complex(pReal), parameter :: TWOPIIMG = cmplx(0.0_pReal,2.0_pReal*PI) !< Re(0.0), Im(2xPi) + + real(pReal), dimension(3,3), parameter :: & MATH_I3 = reshape([& 1.0_pReal,0.0_pReal,0.0_pReal, & 0.0_pReal,1.0_pReal,0.0_pReal, & @@ -75,7 +83,7 @@ module math !--------------------------------------------------------------------------------------------------- private :: & - math_check + unitTest contains @@ -116,14 +124,15 @@ subroutine math_init write(6,'(a,4(/,26x,f17.14),/)') ' start of random sequence: ', randTest call random_seed(put = randInit) - call math_check + call unitTest end subroutine math_init + !-------------------------------------------------------------------------------------------------- !> @brief check correctness of (some) math functions !-------------------------------------------------------------------------------------------------- -subroutine math_check +subroutine unitTest use IO, only: IO_error character(len=64) :: error_msg @@ -145,7 +154,7 @@ subroutine math_check call IO_error(401,ext_msg=error_msg) endif -end subroutine math_check +end subroutine unitTest !-------------------------------------------------------------------------------------------------- @@ -274,6 +283,7 @@ pure function math_identity2nd(dimen) end function math_identity2nd + !-------------------------------------------------------------------------------------------------- !> @brief symmetric fourth rank identity tensor of specified dimension ! from http://en.wikipedia.org/wiki/Tensor_derivative_(continuum_mechanics)#Derivative_of_a_second-order_tensor_with_respect_to_itself @@ -626,6 +636,7 @@ pure function math_skew33(m) end function math_skew33 + !-------------------------------------------------------------------------------------------------- !> @brief hydrostatic part of a 33 matrix !-------------------------------------------------------------------------------------------------- diff --git a/src/plastic_dislotwin.f90 b/src/plastic_dislotwin.f90 index cb4b3cbae..4d84b503e 100644 --- a/src/plastic_dislotwin.f90 +++ b/src/plastic_dislotwin.f90 @@ -8,17 +8,26 @@ !> @details to be done !-------------------------------------------------------------------------------------------------- module plastic_dislotwin - use prec, only: & - pReal + use prec + use debug + use math + use IO + use material + use config + use lattice +#if defined(PETSc) || defined(DAMASK_HDF5) + use results +#endif implicit none private + integer, dimension(:,:), allocatable, target, public :: & plastic_dislotwin_sizePostResult !< size of each post result output character(len=64), dimension(:,:), allocatable, target, public :: & plastic_dislotwin_output !< name of each post result output - real(pReal), parameter, private :: & + real(pReal), parameter :: & kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin enum, bind(c) @@ -39,7 +48,7 @@ module plastic_dislotwin f_tr_ID end enum - type, private :: tParameters + type :: tParameters real(pReal) :: & mu, & nu, & @@ -119,7 +128,7 @@ module plastic_dislotwin dipoleFormation !< flag indicating consideration of dipole formation end type !< container type for internal constitutive parameters - type, private :: tDislotwinState + type :: tDislotwinState real(pReal), dimension(:,:), pointer :: & rho_mob, & rho_dip, & @@ -128,7 +137,7 @@ module plastic_dislotwin f_tr end type tDislotwinState - type, private :: tDislotwinMicrostructure + type :: tDislotwinMicrostructure real(pReal), dimension(:,:), allocatable :: & Lambda_sl, & !* mean free path between 2 obstacles seen by a moving dislocation Lambda_tw, & !* mean free path between 2 obstacles seen by a growing twin @@ -144,11 +153,11 @@ module plastic_dislotwin !-------------------------------------------------------------------------------------------------- ! containers for parameters and state - type(tParameters), allocatable, dimension(:), private :: param - type(tDislotwinState), allocatable, dimension(:), private :: & + type(tParameters), allocatable, dimension(:) :: param + type(tDislotwinState), allocatable, dimension(:) :: & dotState, & state - type(tDislotwinMicrostructure), allocatable, dimension(:), private :: dependentState + type(tDislotwinMicrostructure), allocatable, dimension(:) :: dependentState public :: & plastic_dislotwin_init, & @@ -158,10 +167,6 @@ module plastic_dislotwin plastic_dislotwin_dotState, & plastic_dislotwin_postResults, & plastic_dislotwin_results - private :: & - kinetics_slip, & - kinetics_twin, & - kinetics_trans contains @@ -171,24 +176,6 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine plastic_dislotwin_init - use prec, only: & - pStringLen, & - dEq0, & - dNeq0, & - dNeq - use debug, only: & - debug_level,& - debug_constitutive,& - debug_levelBasic - use math, only: & - math_expand,& - PI - use IO, only: & - IO_error - use material - use config, only: & - config_phase - use lattice integer :: & Ninstance, & @@ -591,10 +578,6 @@ end subroutine plastic_dislotwin_init !> @brief returns the homogenized elasticity matrix !-------------------------------------------------------------------------------------------------- function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC) - use material, only: & - material_phase, & - phase_plasticityInstance, & - phasememberAt real(pReal), dimension(6,6) :: & homogenizedC @@ -634,14 +617,6 @@ end function plastic_dislotwin_homogenizedC !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) - use prec, only: & - tol_math_check, & - dNeq0 - use math, only: & - math_eigenValuesVectorsSym, & - math_outer, & - math_symmetric33, & - math_mul33xx33 real(pReal), dimension(3,3), intent(out) :: Lp real(pReal), dimension(3,3,3,3), intent(out) :: dLp_dMp @@ -757,13 +732,6 @@ end subroutine plastic_dislotwin_LpAndItsTangent !> @brief calculates the rate of change of microstructure !-------------------------------------------------------------------------------------------------- subroutine plastic_dislotwin_dotState(Mp,T,instance,of) - use prec, only: & - tol_math_check, & - dEq0 - use math, only: & - math_clip, & - math_mul33xx33, & - PI real(pReal), dimension(3,3), intent(in):: & Mp !< Mandel stress @@ -854,8 +822,6 @@ end subroutine plastic_dislotwin_dotState !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- subroutine plastic_dislotwin_dependentState(T,instance,of) - use math, only: & - PI integer, intent(in) :: & instance, & @@ -868,13 +834,13 @@ subroutine plastic_dislotwin_dependentState(T,instance,of) real(pReal) :: & sumf_twin,SFE,sumf_trans real(pReal), dimension(param(instance)%sum_N_sl) :: & - inv_lambda_sl_sl, & !* 1/mean free distance between 2 forest dislocations seen by a moving dislocation - inv_lambda_sl_tw, & !* 1/mean free distance between 2 twin stacks from different systems seen by a moving dislocation - inv_lambda_sl_tr !* 1/mean free distance between 2 martensite lamellar from different systems seen by a moving dislocation + inv_lambda_sl_sl, & !< 1/mean free distance between 2 forest dislocations seen by a moving dislocation + inv_lambda_sl_tw, & !< 1/mean free distance between 2 twin stacks from different systems seen by a moving dislocation + inv_lambda_sl_tr !< 1/mean free distance between 2 martensite lamellar from different systems seen by a moving dislocation real(pReal), dimension(param(instance)%sum_N_tw) :: & - inv_lambda_tw_tw !* 1/mean free distance between 2 twin stacks from different systems seen by a growing twin + inv_lambda_tw_tw !< 1/mean free distance between 2 twin stacks from different systems seen by a growing twin real(pReal), dimension(param(instance)%sum_N_tr) :: & - inv_lambda_tr_tr !* 1/mean free distance between 2 martensite stacks from different systems seen by a growing martensite (1/lambda_trans) + inv_lambda_tr_tr !< 1/mean free distance between 2 martensite stacks from different systems seen by a growing martensite real(pReal), dimension(:), allocatable :: & x0, & @@ -967,12 +933,6 @@ end subroutine plastic_dislotwin_dependentState !> @brief return array of constitutive results !-------------------------------------------------------------------------------------------------- function plastic_dislotwin_postResults(Mp,T,instance,of) result(postResults) - use prec, only: & - tol_math_check, & - dEq0 - use math, only: & - PI, & - math_mul33xx33 real(pReal), dimension(3,3),intent(in) :: & Mp !< 2nd Piola Kirchhoff stress tensor in Mandel notation @@ -1050,8 +1010,6 @@ end function plastic_dislotwin_postResults !-------------------------------------------------------------------------------------------------- subroutine plastic_dislotwin_results(instance,group) #if defined(PETSc) || defined(DAMASK_HDF5) - use results, only: & - results_writeDataset integer, intent(in) :: instance character(len=*) :: group @@ -1112,11 +1070,6 @@ end subroutine plastic_dislotwin_results !-------------------------------------------------------------------------------------------------- pure subroutine kinetics_slip(Mp,T,instance,of, & dot_gamma_sl,ddot_gamma_dtau_slip,tau_slip) - use prec, only: & - tol_math_check, & - dNeq0 - use math, only: & - math_mul33xx33 real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -1190,11 +1143,6 @@ end subroutine kinetics_slip !-------------------------------------------------------------------------------------------------- pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,of,& dot_gamma_twin,ddot_gamma_dtau_twin) - use prec, only: & - tol_math_check, & - dNeq0 - use math, only: & - math_mul33xx33 real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -1261,11 +1209,6 @@ end subroutine kinetics_twin !-------------------------------------------------------------------------------------------------- pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,& dot_gamma_tr,ddot_gamma_dtau_trans) - use prec, only: & - tol_math_check, & - dNeq0 - use math, only: & - math_mul33xx33 real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress diff --git a/src/plastic_isotropic.f90 b/src/plastic_isotropic.f90 index c572f0ded..46d0905dc 100644 --- a/src/plastic_isotropic.f90 +++ b/src/plastic_isotropic.f90 @@ -8,11 +8,19 @@ !! untextured polycrystal !-------------------------------------------------------------------------------------------------- module plastic_isotropic - use prec, only: & - pReal + use prec + use debug + use math + use IO + use material + use config +#if defined(PETSc) || defined(DAMASK_HDF5) + use results +#endif implicit none private + integer, dimension(:,:), allocatable, target, public :: & plastic_isotropic_sizePostResult !< size of each post result output character(len=64), dimension(:,:), allocatable, target, public :: & @@ -25,7 +33,7 @@ module plastic_isotropic dot_gamma_ID end enum - type, private :: tParameters + type :: tParameters real(pReal) :: & M, & !< Taylor factor xi_0, & !< initial critical stress @@ -49,7 +57,7 @@ module plastic_isotropic dilatation end type tParameters - type, private :: tIsotropicState + type :: tIsotropicState real(pReal), pointer, dimension(:) :: & xi, & gamma @@ -57,8 +65,8 @@ module plastic_isotropic !-------------------------------------------------------------------------------------------------- ! containers for parameters and state - type(tParameters), allocatable, dimension(:), private :: param - type(tIsotropicState), allocatable, dimension(:), private :: & + type(tParameters), allocatable, dimension(:) :: param + type(tIsotropicState), allocatable, dimension(:) :: & dotState, & state @@ -77,25 +85,7 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine plastic_isotropic_init - use prec, only: & - pStringLen - use debug, only: & -#ifdef DEBUG - debug_e, & - debug_i, & - debug_g, & - debug_levelExtensive, & -#endif - debug_level, & - debug_constitutive, & - debug_levelBasic - use IO, only: & - IO_error - use material - use config, only: & - config_phase - use lattice - + integer :: & Ninstance, & p, i, & @@ -235,16 +225,6 @@ end subroutine plastic_isotropic_init !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) -#ifdef DEBUG - use debug, only: & - debug_level, & - debug_constitutive,& - debug_levelExtensive, & - debug_levelSelective -#endif - use math, only: & - math_deviatoric33, & - math_mul33xx33 real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient @@ -307,10 +287,6 @@ end subroutine plastic_isotropic_LpAndItsTangent ! ToDo: Rename Tstar to Mi? !-------------------------------------------------------------------------------------------------- subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar,Tstar,instance,of) - use math, only: & - math_I3, & - math_spherical33, & - math_mul33xx33 real(pReal), dimension(3,3), intent(out) :: & Li !< inleastic velocity gradient @@ -362,11 +338,6 @@ subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar,Tstar,instance,of) !> @brief calculates the rate of change of microstructure !-------------------------------------------------------------------------------------------------- subroutine plastic_isotropic_dotState(Mp,instance,of) - use prec, only: & - dEq0 - use math, only: & - math_mul33xx33, & - math_deviatoric33 real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -416,9 +387,6 @@ end subroutine plastic_isotropic_dotState !> @brief return array of constitutive results !-------------------------------------------------------------------------------------------------- function plastic_isotropic_postResults(Mp,instance,of) result(postResults) - use math, only: & - math_mul33xx33, & - math_deviatoric33 real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -468,7 +436,6 @@ end function plastic_isotropic_postResults !-------------------------------------------------------------------------------------------------- subroutine plastic_isotropic_results(instance,group) #if defined(PETSc) || defined(DAMASKHDF5) - use results integer, intent(in) :: instance character(len=*), intent(in) :: group diff --git a/src/plastic_kinematichardening.f90 b/src/plastic_kinematichardening.f90 index 861b98da3..ab68eb176 100644 --- a/src/plastic_kinematichardening.f90 +++ b/src/plastic_kinematichardening.f90 @@ -6,11 +6,20 @@ !! and a Voce-type kinematic hardening rule !-------------------------------------------------------------------------------------------------- module plastic_kinehardening - use prec, only: & - pReal + use prec + use debug + use math + use IO + use material + use config + use lattice +#if defined(PETSc) || defined(DAMASK_HDF5) + use results +#endif implicit none private + integer, dimension(:,:), allocatable, target, public :: & plastic_kinehardening_sizePostResult !< size of each post result output character(len=64), dimension(:,:), allocatable, target, public :: & @@ -29,7 +38,7 @@ module plastic_kinehardening resolvedstress_ID end enum - type, private :: tParameters + type :: tParameters real(pReal) :: & gdot0, & !< reference shear strain rate for slip n, & !< stress exponent for slip @@ -59,7 +68,7 @@ module plastic_kinehardening outputID !< ID of each post result output end type tParameters - type, private :: tKinehardeningState + type :: tKinehardeningState real(pReal), pointer, dimension(:,:) :: & !< vectors along NipcMyInstance crss, & !< critical resolved stress crss_back, & !< critical resolved back stress @@ -71,8 +80,8 @@ module plastic_kinehardening !-------------------------------------------------------------------------------------------------- ! containers for parameters and state - type(tParameters), allocatable, dimension(:), private :: param - type(tKinehardeningState), allocatable, dimension(:), private :: & + type(tParameters), allocatable, dimension(:) :: param + type(tKinehardeningState), allocatable, dimension(:) :: & dotState, & deltaState, & state @@ -84,8 +93,6 @@ module plastic_kinehardening plastic_kinehardening_deltaState, & plastic_kinehardening_postResults, & plastic_kinehardening_results - private :: & - kinetics contains @@ -95,27 +102,6 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine plastic_kinehardening_init - use prec, only: & - dEq0, & - pStringLen - use debug, only: & -#ifdef DEBUG - debug_e, & - debug_i, & - debug_g, & - debug_levelExtensive, & -#endif - debug_level, & - debug_constitutive,& - debug_levelBasic - use math, only: & - math_expand - use IO, only: & - IO_error - use material - use config, only: & - config_phase - use lattice integer :: & Ninstance, & @@ -417,16 +403,6 @@ end subroutine plastic_kinehardening_dotState !> @brief calculates (instantaneous) incremental change of microstructure !-------------------------------------------------------------------------------------------------- subroutine plastic_kinehardening_deltaState(Mp,instance,of) - use prec, only: & - dNeq, & - dEq0 -#ifdef DEBUG - use debug, only: & - debug_level, & - debug_constitutive,& - debug_levelExtensive, & - debug_levelSelective -#endif real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -475,8 +451,6 @@ end subroutine plastic_kinehardening_deltaState !> @brief return array of constitutive results !-------------------------------------------------------------------------------------------------- function plastic_kinehardening_postResults(Mp,instance,of) result(postResults) - use math, only: & - math_mul33xx33 real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -535,8 +509,6 @@ end function plastic_kinehardening_postResults !-------------------------------------------------------------------------------------------------- subroutine plastic_kinehardening_results(instance,group) #if defined(PETSc) || defined(DAMASK_HDF5) - use results, only: & - results_writeDataset integer, intent(in) :: instance character(len=*) :: group @@ -585,10 +557,6 @@ end subroutine plastic_kinehardening_results !-------------------------------------------------------------------------------------------------- pure subroutine kinetics(Mp,instance,of, & gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) - use prec, only: & - dNeq0 - use math, only: & - math_mul33xx33 real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress diff --git a/src/plastic_none.f90 b/src/plastic_none.f90 index 4b14266f1..894cc9a40 100644 --- a/src/plastic_none.f90 +++ b/src/plastic_none.f90 @@ -5,6 +5,8 @@ !> @brief Dummy plasticity for purely elastic material !-------------------------------------------------------------------------------------------------- module plastic_none + use material + use debug implicit none private @@ -19,11 +21,6 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine plastic_none_init - use debug, only: & - debug_level, & - debug_constitutive, & - debug_levelBasic - use material integer :: & Ninstance, & diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index 196129f64..a31891573 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -5,11 +5,20 @@ !> @brief phenomenological crystal plasticity formulation using a powerlaw fitting !-------------------------------------------------------------------------------------------------- module plastic_phenopowerlaw - use prec, only: & - pReal + use prec + use debug + use math + use IO + use material + use config + use lattice +#if defined(PETSc) || defined(DAMASK_HDF5) + use results +#endif implicit none private + integer, dimension(:,:), allocatable, target, public :: & plastic_phenopowerlaw_sizePostResult !< size of each post result output character(len=64), dimension(:,:), allocatable, target, public :: & @@ -28,7 +37,7 @@ module plastic_phenopowerlaw resolvedstress_twin_ID end enum - type, private :: tParameters + type :: tParameters real(pReal) :: & gdot0_slip, & !< reference shear strain rate for slip gdot0_twin, & !< reference shear strain rate for twin @@ -73,7 +82,7 @@ module plastic_phenopowerlaw outputID !< ID of each post result output end type tParameters - type, private :: tPhenopowerlawState + type :: tPhenopowerlawState real(pReal), pointer, dimension(:,:) :: & xi_slip, & xi_twin, & @@ -83,8 +92,8 @@ module plastic_phenopowerlaw !-------------------------------------------------------------------------------------------------- ! containers for parameters and state - type(tParameters), allocatable, dimension(:), private :: param - type(tPhenopowerlawState), allocatable, dimension(:), private :: & + type(tParameters), allocatable, dimension(:) :: param + type(tPhenopowerlawState), allocatable, dimension(:) :: & dotState, & state @@ -94,9 +103,6 @@ module plastic_phenopowerlaw plastic_phenopowerlaw_dotState, & plastic_phenopowerlaw_postResults, & plastic_phenopowerlaw_results - private :: & - kinetics_slip, & - kinetics_twin contains @@ -106,20 +112,6 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine plastic_phenopowerlaw_init - use prec, only: & - pStringLen - use debug, only: & - debug_level, & - debug_constitutive,& - debug_levelBasic - use math, only: & - math_expand - use IO, only: & - IO_error - use material - use config, only: & - config_phase - use lattice integer :: & Ninstance, & @@ -484,8 +476,6 @@ end subroutine plastic_phenopowerlaw_dotState !> @brief return array of constitutive results !-------------------------------------------------------------------------------------------------- function plastic_phenopowerlaw_postResults(Mp,instance,of) result(postResults) - use math, only: & - math_mul33xx33 real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -552,8 +542,6 @@ end function plastic_phenopowerlaw_postResults !-------------------------------------------------------------------------------------------------- subroutine plastic_phenopowerlaw_results(instance,group) #if defined(PETSc) || defined(DAMASK_HDF5) - use results, only: & - results_writeDataset integer, intent(in) :: instance character(len=*), intent(in) :: group @@ -598,10 +586,6 @@ end subroutine plastic_phenopowerlaw_results !-------------------------------------------------------------------------------------------------- pure subroutine kinetics_slip(Mp,instance,of, & gdot_slip_pos,gdot_slip_neg,dgdot_dtau_slip_pos,dgdot_dtau_slip_neg) - use prec, only: & - dNeq0 - use math, only: & - math_mul33xx33 real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -674,10 +658,6 @@ end subroutine kinetics_slip !-------------------------------------------------------------------------------------------------- pure subroutine kinetics_twin(Mp,instance,of,& gdot_twin,dgdot_dtau_twin) - use prec, only: & - dNeq0 - use math, only: & - math_mul33xx33 real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress diff --git a/src/quaternions.f90 b/src/quaternions.f90 index fa9c13f38..dc894bdfa 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -3,27 +3,27 @@ ! Modified 2017-2019, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH ! All rights reserved. ! -! Redistribution and use in source and binary forms, with or without modification, are +! Redistribution and use in source and binary forms, with or without modification, are ! permitted provided that the following conditions are met: ! -! - Redistributions of source code must retain the above copyright notice, this list +! - Redistributions of source code must retain the above copyright notice, this list ! of conditions and the following disclaimer. -! - Redistributions in binary form must reproduce the above copyright notice, this -! list of conditions and the following disclaimer in the documentation and/or +! - Redistributions in binary form must reproduce the above copyright notice, this +! list of conditions and the following disclaimer in the documentation and/or ! other materials provided with the distribution. -! - Neither the names of Marc De Graef, Carnegie Mellon University nor the names -! of its contributors may be used to endorse or promote products derived from +! - Neither the names of Marc De Graef, Carnegie Mellon University nor the names +! of its contributors may be used to endorse or promote products derived from ! this software without specific prior written permission. ! -! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" -! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE -! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE -! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL -! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER -! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, -! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE ! USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! ################################################################### @@ -34,58 +34,57 @@ !> @details w is the real part, (x, y, z) are the imaginary parts. !--------------------------------------------------------------------------------------------------- module quaternions - use prec, only: & - pReal - use future + use prec + use future - implicit none - public - - real(pReal), parameter, public :: P = -1.0_pReal !< parameter for orientation conversion. - - type, public :: quaternion - real(pReal) :: w = 0.0_pReal - real(pReal) :: x = 0.0_pReal - real(pReal) :: y = 0.0_pReal - real(pReal) :: z = 0.0_pReal - + implicit none + public - contains - procedure, private :: add__ - procedure, private :: pos__ - generic, public :: operator(+) => add__,pos__ + real(pReal), parameter, public :: P = -1.0_pReal !< parameter for orientation conversion. - procedure, private :: sub__ - procedure, private :: neg__ - generic, public :: operator(-) => sub__,neg__ + type, public :: quaternion + real(pReal) :: w = 0.0_pReal + real(pReal) :: x = 0.0_pReal + real(pReal) :: y = 0.0_pReal + real(pReal) :: z = 0.0_pReal - procedure, private :: mul_quat__ - procedure, private :: mul_scal__ - generic, public :: operator(*) => mul_quat__, mul_scal__ - procedure, private :: div_quat__ - procedure, private :: div_scal__ - generic, public :: operator(/) => div_quat__, div_scal__ + contains + procedure, private :: add__ + procedure, private :: pos__ + generic, public :: operator(+) => add__,pos__ - procedure, private :: eq__ - generic, public :: operator(==) => eq__ + procedure, private :: sub__ + procedure, private :: neg__ + generic, public :: operator(-) => sub__,neg__ - procedure, private :: neq__ - generic, public :: operator(/=) => neq__ + procedure, private :: mul_quat__ + procedure, private :: mul_scal__ + generic, public :: operator(*) => mul_quat__, mul_scal__ - procedure, private :: pow_quat__ - procedure, private :: pow_scal__ - generic, public :: operator(**) => pow_quat__, pow_scal__ + procedure, private :: div_quat__ + procedure, private :: div_scal__ + generic, public :: operator(/) => div_quat__, div_scal__ - procedure, public :: abs__ - procedure, public :: dot_product__ - procedure, public :: conjg__ - procedure, public :: exp__ - procedure, public :: log__ + procedure, private :: eq__ + generic, public :: operator(==) => eq__ - procedure, public :: homomorphed => quat_homomorphed + procedure, private :: neq__ + generic, public :: operator(/=) => neq__ - end type + procedure, private :: pow_quat__ + procedure, private :: pow_scal__ + generic, public :: operator(**) => pow_quat__, pow_scal__ + + procedure, public :: abs__ + procedure, public :: dot_product__ + procedure, public :: conjg__ + procedure, public :: exp__ + procedure, public :: log__ + + procedure, public :: homomorphed => quat_homomorphed + + end type interface assignment (=) module procedure assign_quat__ @@ -124,12 +123,12 @@ contains !--------------------------------------------------------------------------------------------------- type(quaternion) pure function init__(array) - real(pReal), intent(in), dimension(4) :: array - - init__%w=array(1) - init__%x=array(2) - init__%y=array(3) - init__%z=array(4) + real(pReal), intent(in), dimension(4) :: array + + init__%w=array(1) + init__%x=array(2) + init__%y=array(3) + init__%z=array(4) end function init__ @@ -139,14 +138,14 @@ end function init__ !--------------------------------------------------------------------------------------------------- elemental subroutine assign_quat__(self,other) - type(quaternion), intent(out) :: self - type(quaternion), intent(in) :: other - - self%w = other%w - self%x = other%x - self%y = other%y - self%z = other%z - + type(quaternion), intent(out) :: self + type(quaternion), intent(in) :: other + + self%w = other%w + self%x = other%x + self%y = other%y + self%z = other%z + end subroutine assign_quat__ @@ -155,14 +154,14 @@ end subroutine assign_quat__ !--------------------------------------------------------------------------------------------------- pure subroutine assign_vec__(self,other) - type(quaternion), intent(out) :: self - real(pReal), intent(in), dimension(4) :: other - - self%w = other(1) - self%x = other(2) - self%y = other(3) - self%z = other(4) - + type(quaternion), intent(out) :: self + real(pReal), intent(in), dimension(4) :: other + + self%w = other(1) + self%x = other(2) + self%y = other(3) + self%z = other(4) + end subroutine assign_vec__ @@ -171,13 +170,13 @@ end subroutine assign_vec__ !--------------------------------------------------------------------------------------------------- type(quaternion) elemental function add__(self,other) - class(quaternion), intent(in) :: self,other - - add__%w = self%w + other%w - add__%x = self%x + other%x - add__%y = self%y + other%y - add__%z = self%z + other%z - + class(quaternion), intent(in) :: self,other + + add__%w = self%w + other%w + add__%x = self%x + other%x + add__%y = self%y + other%y + add__%z = self%z + other%z + end function add__ @@ -186,13 +185,13 @@ end function add__ !--------------------------------------------------------------------------------------------------- type(quaternion) elemental function pos__(self) - class(quaternion), intent(in) :: self - - pos__%w = self%w - pos__%x = self%x - pos__%y = self%y - pos__%z = self%z - + class(quaternion), intent(in) :: self + + pos__%w = self%w + pos__%x = self%x + pos__%y = self%y + pos__%z = self%z + end function pos__ @@ -201,13 +200,13 @@ end function pos__ !--------------------------------------------------------------------------------------------------- type(quaternion) elemental function sub__(self,other) - class(quaternion), intent(in) :: self,other - - sub__%w = self%w - other%w - sub__%x = self%x - other%x - sub__%y = self%y - other%y - sub__%z = self%z - other%z - + class(quaternion), intent(in) :: self,other + + sub__%w = self%w - other%w + sub__%x = self%x - other%x + sub__%y = self%y - other%y + sub__%z = self%z - other%z + end function sub__ @@ -216,13 +215,13 @@ end function sub__ !--------------------------------------------------------------------------------------------------- type(quaternion) elemental function neg__(self) - class(quaternion), intent(in) :: self - - neg__%w = -self%w - neg__%x = -self%x - neg__%y = -self%y - neg__%z = -self%z - + class(quaternion), intent(in) :: self + + neg__%w = -self%w + neg__%x = -self%x + neg__%y = -self%y + neg__%z = -self%z + end function neg__ @@ -231,13 +230,13 @@ end function neg__ !--------------------------------------------------------------------------------------------------- type(quaternion) elemental function mul_quat__(self,other) - class(quaternion), intent(in) :: self, other + class(quaternion), intent(in) :: self, other + + mul_quat__%w = self%w*other%w - self%x*other%x - self%y*other%y - self%z*other%z + mul_quat__%x = self%w*other%x + self%x*other%w + P * (self%y*other%z - self%z*other%y) + mul_quat__%y = self%w*other%y + self%y*other%w + P * (self%z*other%x - self%x*other%z) + mul_quat__%z = self%w*other%z + self%z*other%w + P * (self%x*other%y - self%y*other%x) - mul_quat__%w = self%w*other%w - self%x*other%x - self%y*other%y - self%z*other%z - mul_quat__%x = self%w*other%x + self%x*other%w + P * (self%y*other%z - self%z*other%y) - mul_quat__%y = self%w*other%y + self%y*other%w + P * (self%z*other%x - self%x*other%z) - mul_quat__%z = self%w*other%z + self%z*other%w + P * (self%x*other%y - self%y*other%x) - end function mul_quat__ @@ -246,14 +245,14 @@ end function mul_quat__ !--------------------------------------------------------------------------------------------------- type(quaternion) elemental function mul_scal__(self,scal) - class(quaternion), intent(in) :: self - real(pReal), intent(in) :: scal + class(quaternion), intent(in) :: self + real(pReal), intent(in) :: scal + + mul_scal__%w = self%w*scal + mul_scal__%x = self%x*scal + mul_scal__%y = self%y*scal + mul_scal__%z = self%z*scal - mul_scal__%w = self%w*scal - mul_scal__%x = self%x*scal - mul_scal__%y = self%y*scal - mul_scal__%z = self%z*scal - end function mul_scal__ @@ -262,9 +261,9 @@ end function mul_scal__ !--------------------------------------------------------------------------------------------------- type(quaternion) elemental function div_quat__(self,other) - class(quaternion), intent(in) :: self, other + class(quaternion), intent(in) :: self, other - div_quat__ = self * (conjg(other)/(abs(other)**2.0_pReal)) + div_quat__ = self * (conjg(other)/(abs(other)**2.0_pReal)) end function div_quat__ @@ -274,10 +273,10 @@ end function div_quat__ !--------------------------------------------------------------------------------------------------- type(quaternion) elemental function div_scal__(self,scal) - class(quaternion), intent(in) :: self - real(pReal), intent(in) :: scal + class(quaternion), intent(in) :: self + real(pReal), intent(in) :: scal - div_scal__ = [self%w,self%x,self%y,self%z]/scal + div_scal__ = [self%w,self%x,self%y,self%z]/scal end function div_scal__ @@ -286,14 +285,12 @@ end function div_scal__ !> equality of two quaternions !--------------------------------------------------------------------------------------------------- logical elemental function eq__(self,other) - use prec, only: & - dEq - class(quaternion), intent(in) :: self,other + class(quaternion), intent(in) :: self,other + + eq__ = all(dEq([ self%w, self%x, self%y, self%z], & + [other%w,other%x,other%y,other%z])) - eq__ = all(dEq([ self%w, self%x, self%y, self%z], & - [other%w,other%x,other%y,other%z])) - end function eq__ @@ -302,10 +299,10 @@ end function eq__ !--------------------------------------------------------------------------------------------------- logical elemental function neq__(self,other) - class(quaternion), intent(in) :: self,other + class(quaternion), intent(in) :: self,other + + neq__ = .not. self%eq__(other) - neq__ = .not. self%eq__(other) - end function neq__ @@ -314,11 +311,11 @@ end function neq__ !--------------------------------------------------------------------------------------------------- type(quaternion) elemental function pow_scal__(self,expon) - class(quaternion), intent(in) :: self - real(pReal), intent(in) :: expon - - pow_scal__ = exp(log(self)*expon) - + class(quaternion), intent(in) :: self + real(pReal), intent(in) :: expon + + pow_scal__ = exp(log(self)*expon) + end function pow_scal__ @@ -327,11 +324,11 @@ end function pow_scal__ !--------------------------------------------------------------------------------------------------- type(quaternion) elemental function pow_quat__(self,expon) - class(quaternion), intent(in) :: self - type(quaternion), intent(in) :: expon - - pow_quat__ = exp(log(self)*expon) - + class(quaternion), intent(in) :: self + type(quaternion), intent(in) :: expon + + pow_quat__ = exp(log(self)*expon) + end function pow_quat__ @@ -341,15 +338,15 @@ end function pow_quat__ !--------------------------------------------------------------------------------------------------- type(quaternion) elemental function exp__(self) - class(quaternion), intent(in) :: self - real(pReal) :: absImag + class(quaternion), intent(in) :: self + real(pReal) :: absImag - absImag = norm2([self%x, self%y, self%z]) + absImag = norm2([self%x, self%y, self%z]) - exp__ = exp(self%w) * [ cos(absImag), & - self%x/absImag * sin(absImag), & - self%y/absImag * sin(absImag), & - self%z/absImag * sin(absImag)] + exp__ = exp(self%w) * [ cos(absImag), & + self%x/absImag * sin(absImag), & + self%y/absImag * sin(absImag), & + self%z/absImag * sin(absImag)] end function exp__ @@ -360,16 +357,16 @@ end function exp__ !--------------------------------------------------------------------------------------------------- type(quaternion) elemental function log__(self) - class(quaternion), intent(in) :: self - real(pReal) :: absImag + class(quaternion), intent(in) :: self + real(pReal) :: absImag - absImag = norm2([self%x, self%y, self%z]) + absImag = norm2([self%x, self%y, self%z]) + + log__ = [log(abs(self)), & + self%x/absImag * acos(self%w/abs(self)), & + self%y/absImag * acos(self%w/abs(self)), & + self%z/absImag * acos(self%w/abs(self))] - log__ = [log(abs(self)), & - self%x/absImag * acos(self%w/abs(self)), & - self%y/absImag * acos(self%w/abs(self)), & - self%z/absImag * acos(self%w/abs(self))] - end function log__ @@ -378,10 +375,10 @@ end function log__ !--------------------------------------------------------------------------------------------------- real(pReal) elemental function abs__(a) - class(quaternion), intent(in) :: a + class(quaternion), intent(in) :: a + + abs__ = norm2([a%w,a%x,a%y,a%z]) - abs__ = norm2([a%w,a%x,a%y,a%z]) - end function abs__ @@ -390,10 +387,10 @@ end function abs__ !--------------------------------------------------------------------------------------------------- real(pReal) elemental function dot_product__(a,b) - class(quaternion), intent(in) :: a,b + class(quaternion), intent(in) :: a,b + + dot_product__ = a%w*b%w + a%x*b%x + a%y*b%y + a%z*b%z - dot_product__ = a%w*b%w + a%x*b%x + a%y*b%y + a%z*b%z - end function dot_product__ @@ -402,10 +399,10 @@ end function dot_product__ !--------------------------------------------------------------------------------------------------- type(quaternion) elemental function conjg__(a) - class(quaternion), intent(in) :: a + class(quaternion), intent(in) :: a + + conjg__ = quaternion([a%w, -a%x, -a%y, -a%z]) - conjg__ = quaternion([a%w, -a%x, -a%y, -a%z]) - end function conjg__ @@ -414,10 +411,10 @@ end function conjg__ !--------------------------------------------------------------------------------------------------- type(quaternion) elemental function quat_homomorphed(a) - class(quaternion), intent(in) :: a + class(quaternion), intent(in) :: a + + quat_homomorphed = quaternion(-[a%w,a%x,a%y,a%z]) - quat_homomorphed = quaternion(-[a%w,a%x,a%y,a%z]) - end function quat_homomorphed end module quaternions diff --git a/src/results.f90 b/src/results.f90 index 05db831f7..cee86c7da 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -5,6 +5,9 @@ !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !-------------------------------------------------------------------------------------------------- module results + use DAMASK_interface + use rotations + use numerics use HDF5_utilities #ifdef PETSc use PETSC @@ -55,8 +58,6 @@ module results contains subroutine results_init - use DAMASK_interface, only: & - getSolverJobName character(len=pStringLen) :: commandLine @@ -83,9 +84,6 @@ end subroutine results_init !> @brief opens the results file to append data !-------------------------------------------------------------------------------------------------- subroutine results_openJobFile - use DAMASK_interface, only: & - getSolverJobName - resultsFile = HDF5_openFile(trim(getSolverJobName())//'.hdf5','a',.true.) @@ -396,8 +394,6 @@ end subroutine results_writeTensorDataset_int !> @brief stores a scalar dataset in a group !-------------------------------------------------------------------------------------------------- subroutine results_writeScalarDataset_rotation(group,dataset,label,description,lattice_structure) - use rotations, only: & - rotation character(len=*), intent(in) :: label,group,description character(len=*), intent(in), optional :: lattice_structure @@ -428,9 +424,6 @@ end subroutine results_writeScalarDataset_rotation !> @brief adds the unique mapping from spatial position and constituent ID to results !-------------------------------------------------------------------------------------------------- subroutine results_mapping_constituent(phaseAt,memberAt,label) - use numerics, only: & - worldrank, & - worldsize integer, dimension(:,:), intent(in) :: phaseAt !< phase section at (constituent,element) integer, dimension(:,:,:), intent(in) :: memberAt !< phase member at (constituent,IP,element) @@ -566,9 +559,6 @@ end subroutine results_mapping_constituent !> @brief adds the unique mapping from spatial position and constituent ID to results !-------------------------------------------------------------------------------------------------- subroutine results_mapping_materialpoint(homogenizationAt,memberAt,label) - use numerics, only: & - worldrank, & - worldsize integer, dimension(:), intent(in) :: homogenizationAt !< homogenization section at (element) integer, dimension(:,:), intent(in) :: memberAt !< homogenization member at (IP,element) diff --git a/src/rotations.f90 b/src/rotations.f90 index 69529ed24..3a64f27b9 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -46,12 +46,15 @@ !--------------------------------------------------------------------------------------------------- module rotations - use prec, only: & - pReal + use prec + use IO + use math + use Lambert use quaternions implicit none private + type, public :: rotation type(quaternion), private :: q contains @@ -148,8 +151,6 @@ end subroutine !> @details: rotation is based on unit quaternion or rotation matrix (fallback) !--------------------------------------------------------------------------------------------------- function rotVector(self,v,active) - use prec, only: & - dEq real(pReal), dimension(3) :: rotVector class(rotation), intent(in) :: self @@ -260,10 +261,6 @@ end function qu2om !> @brief convert unit quaternion to Euler angles !--------------------------------------------------------------------------------------------------- pure function qu2eu(qu) result(eu) - use prec, only: & - dEq0 - use math, only: & - PI type(quaternion), intent(in) :: qu real(pReal), dimension(3) :: eu @@ -294,12 +291,6 @@ end function qu2eu !> @brief convert unit quaternion to axis angle pair !--------------------------------------------------------------------------------------------------- pure function qu2ax(qu) result(ax) - use prec, only: & - dEq0, & - dNeq0 - use math, only: & - PI, & - math_clip type(quaternion), intent(in) :: qu real(pReal), dimension(4) :: ax @@ -324,13 +315,6 @@ end function qu2ax !> @brief convert unit quaternion to Rodrigues vector !--------------------------------------------------------------------------------------------------- pure function qu2ro(qu) result(ro) - use, intrinsic :: IEEE_ARITHMETIC, only: & - IEEE_value, & - IEEE_positive_inf - use prec, only: & - dEq0 - use math, only: & - math_clip type(quaternion), intent(in) :: qu real(pReal), dimension(4) :: ro @@ -358,10 +342,6 @@ end function qu2ro !> @brief convert unit quaternion to homochoric !--------------------------------------------------------------------------------------------------- pure function qu2ho(qu) result(ho) - use prec, only: & - dEq0 - use math, only: & - math_clip type(quaternion), intent(in) :: qu real(pReal), dimension(3) :: ho @@ -415,8 +395,6 @@ end function om2qu !> @brief orientation matrix to Euler angles !--------------------------------------------------------------------------------------------------- pure function om2eu(om) result(eu) - use math, only: & - PI real(pReal), intent(in), dimension(3,3) :: om real(pReal), dimension(3) :: eu @@ -441,15 +419,6 @@ end function om2eu !> @brief convert orientation matrix to axis angle pair !--------------------------------------------------------------------------------------------------- function om2ax(om) result(ax) - use prec, only: & - dEq0, & - cEq, & - dNeq0 - use IO, only: & - IO_error - use math, only: & - math_clip, & - math_trace33 real(pReal), intent(in) :: om(3,3) real(pReal) :: ax(4) @@ -560,8 +529,6 @@ end function eu2qu !> @brief Euler angles to orientation matrix !--------------------------------------------------------------------------------------------------- pure function eu2om(eu) result(om) - use prec, only: & - dEq0 real(pReal), intent(in), dimension(3) :: eu real(pReal), dimension(3,3) :: om @@ -591,11 +558,6 @@ end function eu2om !> @brief convert euler to axis angle !--------------------------------------------------------------------------------------------------- pure function eu2ax(eu) result(ax) - use prec, only: & - dEq0, & - dEq - use math, only: & - PI real(pReal), intent(in), dimension(3) :: eu real(pReal), dimension(4) :: ax @@ -625,13 +587,6 @@ end function eu2ax !> @brief Euler angles to Rodrigues vector !--------------------------------------------------------------------------------------------------- pure function eu2ro(eu) result(ro) - use prec, only: & - dEq0 - use, intrinsic :: IEEE_ARITHMETIC, only: & - IEEE_value, & - IEEE_positive_inf - use math, only: & - PI real(pReal), intent(in), dimension(3) :: eu real(pReal), dimension(4) :: ro @@ -681,8 +636,6 @@ end function eu2cu !> @brief convert axis angle pair to quaternion !--------------------------------------------------------------------------------------------------- pure function ax2qu(ax) result(qu) - use prec, only: & - dEq0 real(pReal), intent(in), dimension(4) :: ax type(quaternion) :: qu @@ -755,13 +708,6 @@ end function ax2eu !> @brief convert axis angle pair to Rodrigues vector !--------------------------------------------------------------------------------------------------- pure function ax2ro(ax) result(ro) - use, intrinsic :: IEEE_ARITHMETIC, only: & - IEEE_value, & - IEEE_positive_inf - use prec, only: & - dEq0 - use math, only: & - PI real(pReal), intent(in), dimension(4) :: ax real(pReal), dimension(4) :: ro @@ -858,12 +804,6 @@ end function ro2eu !> @brief convert Rodrigues vector to axis angle pair !--------------------------------------------------------------------------------------------------- pure function ro2ax(ro) result(ax) - use, intrinsic :: IEEE_ARITHMETIC, only: & - IEEE_is_finite - use prec, only: & - dEq0 - use math, only: & - PI real(pReal), intent(in), dimension(4) :: ro real(pReal), dimension(4) :: ax @@ -890,12 +830,6 @@ end function ro2ax !> @brief convert Rodrigues vector to homochoric !--------------------------------------------------------------------------------------------------- pure function ro2ho(ro) result(ho) - use, intrinsic :: IEEE_ARITHMETIC, only: & - IEEE_is_finite - use prec, only: & - dEq0 - use math, only: & - PI real(pReal), intent(in), dimension(4) :: ro real(pReal), dimension(3) :: ho @@ -973,8 +907,6 @@ end function ho2eu !> @brief convert homochoric to axis angle pair !--------------------------------------------------------------------------------------------------- pure function ho2ax(ho) result(ax) - use prec, only: & - dEq0 real(pReal), intent(in), dimension(3) :: ho real(pReal), dimension(4) :: ax @@ -1029,13 +961,11 @@ end function ho2ro !> @brief convert homochoric to cubochoric !--------------------------------------------------------------------------------------------------- function ho2cu(ho) result(cu) - use Lambert, only: & - LambertBallToCube real(pReal), intent(in), dimension(3) :: ho real(pReal), dimension(3) :: cu - cu = LambertBallToCube(ho) + cu = Lambert_BallToCube(ho) end function ho2cu @@ -1115,13 +1045,11 @@ end function cu2ro !> @brief convert cubochoric to homochoric !--------------------------------------------------------------------------------------------------- function cu2ho(cu) result(ho) - use Lambert, only: & - LambertCubeToBall real(pReal), intent(in), dimension(3) :: cu real(pReal), dimension(3) :: ho - ho = LambertCubeToBall(cu) + ho = Lambert_CubeToBall(cu) end function cu2ho diff --git a/src/source_damage_anisoBrittle.f90 b/src/source_damage_anisoBrittle.f90 index 494bbc6f0..ccad7c6b0 100644 --- a/src/source_damage_anisoBrittle.f90 +++ b/src/source_damage_anisoBrittle.f90 @@ -5,55 +5,62 @@ !> @details to be done !-------------------------------------------------------------------------------------------------- module source_damage_anisoBrittle - use prec + use prec + use debug + use IO + use math + use material + use config + use lattice - implicit none - private - integer, dimension(:), allocatable, public, protected :: & - source_damage_anisoBrittle_offset, & !< which source is my current source mechanism? - source_damage_anisoBrittle_instance !< instance of source mechanism + implicit none + private - integer, dimension(:,:), allocatable, target, public :: & - source_damage_anisoBrittle_sizePostResult !< size of each post result output + integer, dimension(:), allocatable, public, protected :: & + source_damage_anisoBrittle_offset, & !< which source is my current source mechanism? + source_damage_anisoBrittle_instance !< instance of source mechanism - character(len=64), dimension(:,:), allocatable, target, public :: & - source_damage_anisoBrittle_output !< name of each post result output - - integer, dimension(:,:), allocatable, private :: & - source_damage_anisoBrittle_Ncleavage !< number of cleavage systems per family + integer, dimension(:,:), allocatable, target, public :: & + source_damage_anisoBrittle_sizePostResult !< size of each post result output - enum, bind(c) - enumerator :: undefined_ID, & - damage_drivingforce_ID - end enum + character(len=64), dimension(:,:), allocatable, target, public :: & + source_damage_anisoBrittle_output !< name of each post result output + + integer, dimension(:,:), allocatable :: & + source_damage_anisoBrittle_Ncleavage !< number of cleavage systems per family + + enum, bind(c) + enumerator :: undefined_ID, & + damage_drivingforce_ID + end enum - type, private :: tParameters !< container type for internal constitutive parameters - real(pReal) :: & - aTol, & - sdot_0, & - N - real(pReal), dimension(:), allocatable :: & - critDisp, & - critLoad - real(pReal), dimension(:,:,:,:), allocatable :: & - cleavage_systems - integer :: & - totalNcleavage - integer, dimension(:), allocatable :: & - Ncleavage - integer(kind(undefined_ID)), allocatable, dimension(:) :: & - outputID !< ID of each post result output - end type tParameters + type :: tParameters !< container type for internal constitutive parameters + real(pReal) :: & + aTol, & + sdot_0, & + N + real(pReal), dimension(:), allocatable :: & + critDisp, & + critLoad + real(pReal), dimension(:,:,:,:), allocatable :: & + cleavage_systems + integer :: & + totalNcleavage + integer, dimension(:), allocatable :: & + Ncleavage + integer(kind(undefined_ID)), allocatable, dimension(:) :: & + outputID !< ID of each post result output + end type tParameters - type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) + type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) - public :: & - source_damage_anisoBrittle_init, & - source_damage_anisoBrittle_dotState, & - source_damage_anisobrittle_getRateAndItsTangent, & - source_damage_anisoBrittle_postResults + public :: & + source_damage_anisoBrittle_init, & + source_damage_anisoBrittle_dotState, & + source_damage_anisobrittle_getRateAndItsTangent, & + source_damage_anisoBrittle_postResults contains @@ -63,266 +70,230 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine source_damage_anisoBrittle_init - use debug, only: & - debug_level,& - debug_constitutive,& - debug_levelBasic - use IO, only: & - IO_error - use math, only: & - math_expand - use material, only: & - material_allocateSourceState, & - phase_source, & - phase_Nsources, & - phase_Noutput, & - SOURCE_damage_anisoBrittle_label, & - SOURCE_damage_anisoBrittle_ID, & - material_phase, & - sourceState - use config, only: & - config_phase, & - material_Nphase - use lattice, only: & - lattice_SchmidMatrix_cleavage, & - lattice_maxNcleavageFamily - integer :: Ninstance,phase,instance,source,sourceOffset - integer :: NofMyPhase,p ,i - integer, dimension(0), parameter :: emptyIntArray = [integer::] - character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] - integer(kind(undefined_ID)) :: & - outputID + integer :: Ninstance,phase,instance,source,sourceOffset + integer :: NofMyPhase,p ,i + integer, dimension(0), parameter :: emptyIntArray = [integer::] + character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] + integer(kind(undefined_ID)) :: & + outputID - character(len=pStringLen) :: & - extmsg = '' - character(len=65536), dimension(:), allocatable :: & - outputs + character(len=pStringLen) :: & + extmsg = '' + character(len=65536), dimension(:), allocatable :: & + outputs - write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//' init -+>>>' + write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//' init -+>>>' - Ninstance = int(count(phase_source == SOURCE_damage_anisoBrittle_ID)) - if (Ninstance == 0) return - - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & - write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - - allocate(source_damage_anisoBrittle_offset(material_Nphase), source=0) - allocate(source_damage_anisoBrittle_instance(material_Nphase), source=0) - do phase = 1, material_Nphase - source_damage_anisoBrittle_instance(phase) = count(phase_source(:,1:phase) == source_damage_anisoBrittle_ID) - do source = 1, phase_Nsources(phase) - if (phase_source(source,phase) == source_damage_anisoBrittle_ID) & - source_damage_anisoBrittle_offset(phase) = source - enddo - enddo - - allocate(source_damage_anisoBrittle_sizePostResult(maxval(phase_Noutput),Ninstance), source=0) - allocate(source_damage_anisoBrittle_output(maxval(phase_Noutput),Ninstance)) - source_damage_anisoBrittle_output = '' + Ninstance = count(phase_source == SOURCE_damage_anisoBrittle_ID) + if (Ninstance == 0) return + + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & + write(6,'(a16,1x,i5,/)') '# instances:',Ninstance + + allocate(source_damage_anisoBrittle_offset(material_Nphase), source=0) + allocate(source_damage_anisoBrittle_instance(material_Nphase), source=0) + do phase = 1, material_Nphase + source_damage_anisoBrittle_instance(phase) = count(phase_source(:,1:phase) == source_damage_anisoBrittle_ID) + do source = 1, phase_Nsources(phase) + if (phase_source(source,phase) == source_damage_anisoBrittle_ID) & + source_damage_anisoBrittle_offset(phase) = source + enddo + enddo + + allocate(source_damage_anisoBrittle_sizePostResult(maxval(phase_Noutput),Ninstance), source=0) + allocate(source_damage_anisoBrittle_output(maxval(phase_Noutput),Ninstance)) + source_damage_anisoBrittle_output = '' - allocate(source_damage_anisoBrittle_Ncleavage(lattice_maxNcleavageFamily,Ninstance), source=0) + allocate(source_damage_anisoBrittle_Ncleavage(lattice_maxNcleavageFamily,Ninstance), source=0) - allocate(param(Ninstance)) - - do p=1, size(config_phase) - if (all(phase_source(:,p) /= SOURCE_DAMAGE_ANISOBRITTLE_ID)) cycle - associate(prm => param(source_damage_anisoBrittle_instance(p)), & - config => config_phase(p)) - - prm%aTol = config%getFloat('anisobrittle_atol',defaultVal = 1.0e-3_pReal) + allocate(param(Ninstance)) + + do p=1, size(config_phase) + if (all(phase_source(:,p) /= SOURCE_DAMAGE_ANISOBRITTLE_ID)) cycle + associate(prm => param(source_damage_anisoBrittle_instance(p)), & + config => config_phase(p)) + + prm%aTol = config%getFloat('anisobrittle_atol',defaultVal = 1.0e-3_pReal) - prm%N = config%getFloat('anisobrittle_ratesensitivity') - prm%sdot_0 = config%getFloat('anisobrittle_sdot0') - - ! sanity checks - if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_atol' - - if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_ratesensitivity' - if (prm%sdot_0 <= 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_sdot0' - - prm%Ncleavage = config%getInts('ncleavage',defaultVal=emptyIntArray) + prm%N = config%getFloat('anisobrittle_ratesensitivity') + prm%sdot_0 = config%getFloat('anisobrittle_sdot0') + + ! sanity checks + if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_atol' + + if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_ratesensitivity' + if (prm%sdot_0 <= 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_sdot0' + + prm%Ncleavage = config%getInts('ncleavage',defaultVal=emptyIntArray) - prm%critDisp = config%getFloats('anisobrittle_criticaldisplacement',requiredSize=size(prm%Ncleavage)) - prm%critLoad = config%getFloats('anisobrittle_criticalload', requiredSize=size(prm%Ncleavage)) - - prm%cleavage_systems = lattice_SchmidMatrix_cleavage (prm%Ncleavage,config%getString('lattice_structure'),& - config%getFloat('c/a',defaultVal=0.0_pReal)) + prm%critDisp = config%getFloats('anisobrittle_criticaldisplacement',requiredSize=size(prm%Ncleavage)) + prm%critLoad = config%getFloats('anisobrittle_criticalload', requiredSize=size(prm%Ncleavage)) + + prm%cleavage_systems = lattice_SchmidMatrix_cleavage (prm%Ncleavage,config%getString('lattice_structure'),& + config%getFloat('c/a',defaultVal=0.0_pReal)) - ! expand: family => system - prm%critDisp = math_expand(prm%critDisp, prm%Ncleavage) - prm%critLoad = math_expand(prm%critLoad, prm%Ncleavage) - - if (any(prm%critLoad < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_criticalload' - if (any(prm%critDisp < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_criticaldisplacement' + ! expand: family => system + prm%critDisp = math_expand(prm%critDisp, prm%Ncleavage) + prm%critLoad = math_expand(prm%critLoad, prm%Ncleavage) + + if (any(prm%critLoad < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_criticalload' + if (any(prm%critDisp < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_criticaldisplacement' !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range - if (extmsg /= '') & - call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//')') + if (extmsg /= '') & + call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//')') !-------------------------------------------------------------------------------------------------- ! output pararameters - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - do i=1, size(outputs) - outputID = undefined_ID - select case(outputs(i)) - - case ('anisobrittle_drivingforce') - source_damage_anisoBrittle_sizePostResult(i,source_damage_anisoBrittle_instance(p)) = 1 - source_damage_anisoBrittle_output(i,source_damage_anisoBrittle_instance(p)) = outputs(i) - prm%outputID = [prm%outputID, damage_drivingforce_ID] + outputs = config%getStrings('(output)',defaultVal=emptyStringArray) + allocate(prm%outputID(0)) + do i=1, size(outputs) + outputID = undefined_ID + select case(outputs(i)) + + case ('anisobrittle_drivingforce') + source_damage_anisoBrittle_sizePostResult(i,source_damage_anisoBrittle_instance(p)) = 1 + source_damage_anisoBrittle_output(i,source_damage_anisoBrittle_instance(p)) = outputs(i) + prm%outputID = [prm%outputID, damage_drivingforce_ID] - end select + end select - enddo + enddo - end associate - - phase = p - NofMyPhase=count(material_phase==phase) - instance = source_damage_anisoBrittle_instance(phase) - sourceOffset = source_damage_anisoBrittle_offset(phase) + end associate + + phase = p + NofMyPhase=count(material_phase==phase) + instance = source_damage_anisoBrittle_instance(phase) + sourceOffset = source_damage_anisoBrittle_offset(phase) - call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,0) - sourceState(phase)%p(sourceOffset)%sizePostResults = sum(source_damage_anisoBrittle_sizePostResult(:,instance)) - sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol + call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,0) + sourceState(phase)%p(sourceOffset)%sizePostResults = sum(source_damage_anisoBrittle_sizePostResult(:,instance)) + sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol - source_damage_anisoBrittle_Ncleavage(1:size(param(instance)%Ncleavage),instance) = param(instance)%Ncleavage - enddo + source_damage_anisoBrittle_Ncleavage(1:size(param(instance)%Ncleavage),instance) = param(instance)%Ncleavage + enddo - end subroutine source_damage_anisoBrittle_init + !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) - use math, only: & - math_mul33xx33 - use material, only: & - phaseAt, phasememberAt, & - sourceState, & - material_homogenizationAt, & - damage, & - damageMapping - use lattice, only: & - lattice_Scleavage, & - lattice_maxNcleavageFamily, & - lattice_NcleavageSystem - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), intent(in), dimension(3,3) :: & - S - integer :: & - phase, & - constituent, & - instance, & - sourceOffset, & - damageOffset, & - homog, & - f, i, index_myFamily, index - real(pReal) :: & - traction_d, traction_t, traction_n, traction_crit + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + S + integer :: & + phase, & + constituent, & + instance, & + sourceOffset, & + damageOffset, & + homog, & + f, i, index_myFamily, index + real(pReal) :: & + traction_d, traction_t, traction_n, traction_crit - phase = phaseAt(ipc,ip,el) - constituent = phasememberAt(ipc,ip,el) - instance = source_damage_anisoBrittle_instance(phase) - sourceOffset = source_damage_anisoBrittle_offset(phase) - homog = material_homogenizationAt(el) - damageOffset = damageMapping(homog)%p(ip,el) - - sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal - - index = 1 - do f = 1,lattice_maxNcleavageFamily - index_myFamily = sum(lattice_NcleavageSystem(1:f-1,phase)) ! at which index starts my family - do i = 1,source_damage_anisoBrittle_Ncleavage(f,instance) ! process each (active) cleavage system in family + phase = phaseAt(ipc,ip,el) + constituent = phasememberAt(ipc,ip,el) + instance = source_damage_anisoBrittle_instance(phase) + sourceOffset = source_damage_anisoBrittle_offset(phase) + homog = material_homogenizationAt(el) + damageOffset = damageMapping(homog)%p(ip,el) + + sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal + + index = 1 + do f = 1,lattice_maxNcleavageFamily + index_myFamily = sum(lattice_NcleavageSystem(1:f-1,phase)) ! at which index starts my family + do i = 1,source_damage_anisoBrittle_Ncleavage(f,instance) ! process each (active) cleavage system in family - traction_d = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase)) - traction_t = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase)) - traction_n = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase)) - - traction_crit = param(instance)%critLoad(index)* & - damage(homog)%p(damageOffset)*damage(homog)%p(damageOffset) + traction_d = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase)) + traction_t = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase)) + traction_n = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase)) + + traction_crit = param(instance)%critLoad(index)* & + damage(homog)%p(damageOffset)*damage(homog)%p(damageOffset) - sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = & - sourceState(phase)%p(sourceOffset)%dotState(1,constituent) + & - param(instance)%sdot_0* & - ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**param(instance)%N + & - (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**param(instance)%N + & - (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**param(instance)%N)/ & - param(instance)%critDisp(index) + sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = & + sourceState(phase)%p(sourceOffset)%dotState(1,constituent) + & + param(instance)%sdot_0* & + ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**param(instance)%N + & + (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**param(instance)%N + & + (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**param(instance)%N)/ & + param(instance)%critDisp(index) - index = index + 1 - enddo - enddo + index = index + 1 + enddo + enddo end subroutine source_damage_anisoBrittle_dotState + !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - use material, only: & - sourceState - integer, intent(in) :: & - phase, & - constituent - real(pReal), intent(in) :: & - phi - real(pReal), intent(out) :: & - localphiDot, & - dLocalphiDot_dPhi - integer :: & - sourceOffset + integer, intent(in) :: & + phase, & + constituent + real(pReal), intent(in) :: & + phi + real(pReal), intent(out) :: & + localphiDot, & + dLocalphiDot_dPhi + integer :: & + sourceOffset - sourceOffset = source_damage_anisoBrittle_offset(phase) - - localphiDot = 1.0_pReal & - - sourceState(phase)%p(sourceOffset)%state(1,constituent)*phi - - dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent) + sourceOffset = source_damage_anisoBrittle_offset(phase) + + localphiDot = 1.0_pReal & + - sourceState(phase)%p(sourceOffset)%state(1,constituent)*phi + + dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent) end subroutine source_damage_anisobrittle_getRateAndItsTangent - + + !-------------------------------------------------------------------------------------------------- !> @brief return array of local damage results !-------------------------------------------------------------------------------------------------- function source_damage_anisoBrittle_postResults(phase, constituent) - use material, only: & - sourceState - integer, intent(in) :: & - phase, & - constituent - real(pReal), dimension(sum(source_damage_anisoBrittle_sizePostResult(:, & - source_damage_anisoBrittle_instance(phase)))) :: & - source_damage_anisoBrittle_postResults + integer, intent(in) :: & + phase, & + constituent - integer :: & - instance, sourceOffset, o, c - - instance = source_damage_anisoBrittle_instance(phase) - sourceOffset = source_damage_anisoBrittle_offset(phase) + real(pReal), dimension(sum(source_damage_anisoBrittle_sizePostResult(:, & + source_damage_anisoBrittle_instance(phase)))) :: & + source_damage_anisoBrittle_postResults - c = 0 + integer :: & + instance, sourceOffset, o, c + + instance = source_damage_anisoBrittle_instance(phase) + sourceOffset = source_damage_anisoBrittle_offset(phase) - do o = 1,size(param(instance)%outputID) - select case(param(instance)%outputID(o)) - case (damage_drivingforce_ID) - source_damage_anisoBrittle_postResults(c+1) = & - sourceState(phase)%p(sourceOffset)%state(1,constituent) - c = c + 1 + c = 0 - end select - enddo + do o = 1,size(param(instance)%outputID) + select case(param(instance)%outputID(o)) + case (damage_drivingforce_ID) + source_damage_anisoBrittle_postResults(c+1) = & + sourceState(phase)%p(sourceOffset)%state(1,constituent) + c = c + 1 + + end select + enddo end function source_damage_anisoBrittle_postResults end module source_damage_anisoBrittle diff --git a/src/source_damage_isoBrittle.f90 b/src/source_damage_isoBrittle.f90 index 90aa5089f..3e0e94f82 100644 --- a/src/source_damage_isoBrittle.f90 +++ b/src/source_damage_isoBrittle.f90 @@ -84,7 +84,7 @@ subroutine source_damage_isoBrittle_init write(6,'(/,a)') ' <<<+- source_'//SOURCE_DAMAGE_ISOBRITTLE_LABEL//' init -+>>>' - Ninstance = int(count(phase_source == SOURCE_damage_isoBrittle_ID)) + Ninstance = count(phase_source == SOURCE_damage_isoBrittle_ID) if (Ninstance == 0) return if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & diff --git a/src/source_thermal_dissipation.f90 b/src/source_thermal_dissipation.f90 index 94452eb47..e8464edd0 100644 --- a/src/source_thermal_dissipation.f90 +++ b/src/source_thermal_dissipation.f90 @@ -5,27 +5,30 @@ !> @details to be done !-------------------------------------------------------------------------------------------------- module source_thermal_dissipation - use prec, only: & - pReal + use prec + use debug + use material + use config implicit none private + integer, dimension(:), allocatable, public, protected :: & - source_thermal_dissipation_offset, & !< which source is my current thermal dissipation mechanism? - source_thermal_dissipation_instance !< instance of thermal dissipation source mechanism + source_thermal_dissipation_offset, & !< which source is my current thermal dissipation mechanism? + source_thermal_dissipation_instance !< instance of thermal dissipation source mechanism integer, dimension(:,:), allocatable, target, public :: & - source_thermal_dissipation_sizePostResult !< size of each post result output + source_thermal_dissipation_sizePostResult !< size of each post result output character(len=64), dimension(:,:), allocatable, target, public :: & - source_thermal_dissipation_output !< name of each post result output + source_thermal_dissipation_output !< name of each post result output - type, private :: tParameters !< container type for internal constitutive parameters + type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & kappa end type tParameters - type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) + type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) public :: & @@ -40,21 +43,6 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine source_thermal_dissipation_init - use debug, only: & - debug_level,& - debug_constitutive,& - debug_levelBasic - use material, only: & - material_allocateSourceState, & - phase_source, & - phase_Nsources, & - phase_Noutput, & - SOURCE_thermal_dissipation_label, & - SOURCE_thermal_dissipation_ID, & - material_phase - use config, only: & - config_phase, & - material_Nphase integer :: Ninstance,instance,source,sourceOffset integer :: NofMyPhase,p diff --git a/src/source_thermal_externalheat.f90 b/src/source_thermal_externalheat.f90 index 699902ad3..99d9a6f1f 100644 --- a/src/source_thermal_externalheat.f90 +++ b/src/source_thermal_externalheat.f90 @@ -5,11 +5,14 @@ !> @brief material subroutine for variable heat source !-------------------------------------------------------------------------------------------------- module source_thermal_externalheat - use prec, only: & - pReal + use prec + use debug + use material + use config implicit none private + integer, dimension(:), allocatable, public, protected :: & source_thermal_externalheat_offset, & !< which source is my current thermal dissipation mechanism? source_thermal_externalheat_instance !< instance of thermal dissipation source mechanism @@ -23,7 +26,7 @@ module source_thermal_externalheat integer, dimension(:), allocatable, target, public :: & source_thermal_externalheat_Noutput !< number of outputs per instance of this source - type, private :: tParameters !< container type for internal constitutive parameters + type :: tParameters !< container type for internal constitutive parameters real(pReal), dimension(:), allocatable :: & time, & heat_rate @@ -31,7 +34,7 @@ module source_thermal_externalheat nIntervals end type tParameters - type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) + type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) public :: & @@ -47,22 +50,6 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine source_thermal_externalheat_init - use debug, only: & - debug_level,& - debug_constitutive,& - debug_levelBasic - use material, only: & - material_allocateSourceState, & - material_phase, & - phase_source, & - phase_Nsources, & - phase_Noutput, & - SOURCE_thermal_externalheat_label, & - SOURCE_thermal_externalheat_ID - use config, only: & - config_phase, & - material_Nphase - integer :: maxNinstance,instance,source,sourceOffset,NofMyPhase,p @@ -116,8 +103,6 @@ end subroutine source_thermal_externalheat_init !> @details state only contains current time to linearly interpolate given heat powers !-------------------------------------------------------------------------------------------------- subroutine source_thermal_externalheat_dotState(phase, of) - use material, only: & - sourceState integer, intent(in) :: & phase, & @@ -135,8 +120,6 @@ end subroutine source_thermal_externalheat_dotState !> @brief returns local heat generation rate !-------------------------------------------------------------------------------------------------- subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phase, of) - use material, only: & - sourceState integer, intent(in) :: & phase, & diff --git a/src/thermal_adiabatic.f90 b/src/thermal_adiabatic.f90 index bfc34d1c4..3c9fd0c6e 100644 --- a/src/thermal_adiabatic.f90 +++ b/src/thermal_adiabatic.f90 @@ -3,9 +3,16 @@ !> @brief material subroutine for adiabatic temperature evolution !-------------------------------------------------------------------------------------------------- module thermal_adiabatic - use prec, only: & - pReal - + use prec + use config + use numerics + use material + use source_thermal_dissipation + use source_thermal_externalheat + use crystallite + use lattice + use mesh + implicit none private @@ -21,7 +28,7 @@ module thermal_adiabatic enumerator :: undefined_ID, & temperature_ID end enum - integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & + integer(kind(undefined_ID)), dimension(:,:), allocatable :: & thermal_adiabatic_outputID !< ID of each post result output @@ -41,21 +48,6 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine thermal_adiabatic_init - use material, only: & - thermal_type, & - thermal_typeInstance, & - homogenization_Noutput, & - THERMAL_ADIABATIC_label, & - THERMAL_adiabatic_ID, & - material_homogenizationAt, & - mappingHomogenization, & - thermalState, & - thermalMapping, & - thermal_initialT, & - temperature, & - temperatureRate - use config, only: & - config_homogenization integer :: maxNinstance,section,instance,i,sizeState,NofMyHomog character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] @@ -112,16 +104,6 @@ end subroutine thermal_adiabatic_init !> @brief calculates adiabatic change in temperature based on local heat generation model !-------------------------------------------------------------------------------------------------- function thermal_adiabatic_updateState(subdt, ip, el) - use numerics, only: & - err_thermal_tolAbs, & - err_thermal_tolRel - use material, only: & - material_homogenizationAt, & - mappingHomogenization, & - thermalState, & - temperature, & - temperatureRate, & - thermalMapping integer, intent(in) :: & ip, & !< integration point number @@ -156,28 +138,11 @@ function thermal_adiabatic_updateState(subdt, ip, el) end function thermal_adiabatic_updateState + !-------------------------------------------------------------------------------------------------- !> @brief returns heat generation rate !-------------------------------------------------------------------------------------------------- subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) - use material, only: & - homogenization_Ngrains, & - material_homogenizationAt, & - mappingHomogenization, & - phaseAt, & - phasememberAt, & - thermal_typeInstance, & - phase_Nsources, & - phase_source, & - SOURCE_thermal_dissipation_ID, & - SOURCE_thermal_externalheat_ID - use source_thermal_dissipation, only: & - source_thermal_dissipation_getRateAndItsTangent - use source_thermal_externalheat, only: & - source_thermal_externalheat_getRateAndItsTangent - use crystallite, only: & - crystallite_S, & - crystallite_Lp integer, intent(in) :: & ip, & !< integration point number @@ -230,18 +195,12 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) dTdot_dT = dTdot_dT/real(homogenization_Ngrains(homog),pReal) end subroutine thermal_adiabatic_getSourceAndItsTangent - + + !-------------------------------------------------------------------------------------------------- !> @brief returns homogenized specific heat capacity !-------------------------------------------------------------------------------------------------- function thermal_adiabatic_getSpecificHeat(ip,el) - use lattice, only: & - lattice_specificHeat - use material, only: & - homogenization_Ngrains, & - material_phase - use mesh, only: & - mesh_element integer, intent(in) :: & ip, & !< integration point number @@ -270,13 +229,6 @@ end function thermal_adiabatic_getSpecificHeat !> @brief returns homogenized mass density !-------------------------------------------------------------------------------------------------- function thermal_adiabatic_getMassDensity(ip,el) - use lattice, only: & - lattice_massDensity - use material, only: & - homogenization_Ngrains, & - material_phase - use mesh, only: & - mesh_element integer, intent(in) :: & ip, & !< integration point number @@ -304,8 +256,6 @@ end function thermal_adiabatic_getMassDensity !> @brief return array of thermal results !-------------------------------------------------------------------------------------------------- function thermal_adiabatic_postResults(homog,instance,of) result(postResults) - use material, only: & - temperature integer, intent(in) :: & homog, &