From 4dd99d4c39e3473e9f6146b0f4812ca14b4d0b0e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 28 Feb 2021 19:13:20 +0100 Subject: [PATCH 01/13] solver is selected in load case, not numerics.yaml --- PRIVATE | 2 +- src/grid/DAMASK_grid.f90 | 34 +++++++++----------- src/grid/grid_mech_spectral_polarisation.f90 | 2 +- 3 files changed, 18 insertions(+), 20 deletions(-) diff --git a/PRIVATE b/PRIVATE index 0289c1bbf..738526ac7 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 0289c1bbfec1a1aef77a8cbaeed134035549e738 +Subproject commit 738526ac7ae63678c885faba7621782499f2bce9 diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index d250e2f53..e0eaef892 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -78,8 +78,7 @@ program DAMASK_grid maxCutBack, & !< max number of cut backs stagItMax !< max number of field level staggered iterations character(len=pStringLen) :: & - incInfo, & - loadcase_string + incInfo type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases type(tSolutionState), allocatable, dimension(:) :: solres @@ -98,10 +97,10 @@ program DAMASK_grid quit class (tNode), pointer :: & num_grid, & - debug_grid, & ! pointer to grid debug options config_load, & load_steps, & load_step, & + solver, & step_bc, & step_mech, & step_discretization, & @@ -117,12 +116,6 @@ program DAMASK_grid print*, 'Shanthraj et al., Handbook of Mechanics of Materials, 2019' print*, 'https://doi.org/10.1007/978-981-10-6855-3_80' -!-------------------------------------------------------------------------------------------------- -! initialize field solver information - nActiveFields = 1 - if (any(thermal_type == THERMAL_conduction_ID )) nActiveFields = nActiveFields + 1 - if (any(damage_type == DAMAGE_nonlocal_ID )) nActiveFields = nActiveFields + 1 - allocate(solres(nActiveFields)) !------------------------------------------------------------------------------------------------- ! reading field paramters from numerics file and do sanity checks @@ -132,20 +125,23 @@ program DAMASK_grid if (stagItMax < 0) call IO_error(301,ext_msg='maxStaggeredIter') if (maxCutBack < 0) call IO_error(301,ext_msg='maxCutBack') + + config_load => YAML_parse_file(trim(interface_loadFile)) + solver => config_load%get('solver') !-------------------------------------------------------------------------------------------------- ! assign mechanics solver depending on selected type - debug_grid => config_debug%get('grid',defaultVal=emptyList) - select case (trim(num_grid%get_asString('solver', defaultVal = 'Basic'))) - case ('Basic') + nActiveFields = 1 + select case (solver%get_asString('mechanical')) + case ('spectral_basic') mechanical_init => grid_mechanical_spectral_basic_init mechanical_forward => grid_mechanical_spectral_basic_forward mechanical_solution => grid_mechanical_spectral_basic_solution mechanical_updateCoords => grid_mechanical_spectral_basic_updateCoords mechanical_restartWrite => grid_mechanical_spectral_basic_restartWrite - case ('Polarisation') + case ('spectral_polarization') mechanical_init => grid_mechanical_spectral_polarisation_init mechanical_forward => grid_mechanical_spectral_polarisation_forward mechanical_solution => grid_mechanical_spectral_polarisation_solution @@ -160,14 +156,17 @@ program DAMASK_grid mechanical_restartWrite => grid_mechanical_FEM_restartWrite case default - call IO_error(error_ID = 891, ext_msg = trim(num_grid%get_asString('solver'))) + call IO_error(error_ID = 891, ext_msg = trim(solver%get_asString('mechanical'))) end select +!-------------------------------------------------------------------------------------------------- +! initialize field solver information + if (any(thermal_type == THERMAL_conduction_ID )) nActiveFields = nActiveFields + 1 + if (any(damage_type == DAMAGE_nonlocal_ID )) nActiveFields = nActiveFields + 1 + allocate(solres(nActiveFields)) !-------------------------------------------------------------------------------------------------- -! reading information from load case file and to sanity checks - config_load => YAML_parse_file(trim(interface_loadFile)) load_steps => config_load%get('loadstep') allocate(loadCases(load_steps%length)) ! array of load cases @@ -232,7 +231,6 @@ program DAMASK_grid merge(.true.,.false.,l > 1)) reportAndCheck: if (worldrank == 0) then - write (loadcase_string, '(i0)' ) l print'(/,a,i0)', ' load case: ', l print*, ' estimate_rate:', loadCases(l)%estimate_rate if (loadCases(l)%deformation%myType == 'L') then @@ -292,7 +290,7 @@ program DAMASK_grid if (loadCases(l)%f_restart < huge(0)) & print'(a,i0)', ' f_restart: ', loadCases(l)%f_restart - if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message + if (errorID > 0) call IO_error(error_ID = errorID, el = l) endif reportAndCheck enddo diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 8caf41b31..31f69b4c5 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -116,7 +116,7 @@ subroutine grid_mechanical_spectral_polarisation_init num_grid, & debug_grid - print'(/,a)', ' <<<+- grid_mechanical_spectral_polarisation init -+>>>'; flush(IO_STDOUT) + print'(/,a)', ' <<<+- grid_mechanical_spectral_polarization init -+>>>'; flush(IO_STDOUT) print*, 'Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006' From b2fea6b149a74511c3c5018137abb2a96945cea7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 28 Feb 2021 19:24:44 +0100 Subject: [PATCH 02/13] solver not specific to load case number --- src/grid/DAMASK_grid.f90 | 34 ++++++++++++++++++---------------- 1 file changed, 18 insertions(+), 16 deletions(-) diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index e0eaef892..36ebf51a3 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -37,8 +37,9 @@ program DAMASK_grid f_out, & !< frequency of result writes f_restart !< frequency of restart writes logical :: estimate_rate !< follow trajectory of former loadcase - integer(kind(FIELD_UNDEFINED_ID)), allocatable :: ID(:) end type tLoadCase + + integer(kind(FIELD_UNDEFINED_ID)), allocatable :: ID(:) !-------------------------------------------------------------------------------------------------- ! variables related to information from load case and geom file @@ -165,6 +166,19 @@ program DAMASK_grid if (any(thermal_type == THERMAL_conduction_ID )) nActiveFields = nActiveFields + 1 if (any(damage_type == DAMAGE_nonlocal_ID )) nActiveFields = nActiveFields + 1 allocate(solres(nActiveFields)) + + allocate(ID(nActiveFields)) + field = 1 + ID(field) = FIELD_MECH_ID ! mechanical active by default + thermalActive: if (any(thermal_type == THERMAL_conduction_ID)) then + field = field + 1 + ID(field) = FIELD_THERMAL_ID + endif thermalActive + damageActive: if (any(damage_type == DAMAGE_nonlocal_ID)) then + field = field + 1 + ID(field) = FIELD_DAMAGE_ID + endif damageActive + !-------------------------------------------------------------------------------------------------- @@ -173,18 +187,6 @@ program DAMASK_grid do l = 1, load_steps%length - allocate(loadCases(l)%ID(nActiveFields)) - field = 1 - loadCases(l)%ID(field) = FIELD_MECH_ID ! mechanical active by default - thermalActive: if (any(thermal_type == THERMAL_conduction_ID)) then - field = field + 1 - loadCases(l)%ID(field) = FIELD_THERMAL_ID - endif thermalActive - damageActive: if (any(damage_type == DAMAGE_nonlocal_ID)) then - field = field + 1 - loadCases(l)%ID(field) = FIELD_DAMAGE_ID - endif damageActive - load_step => load_steps%get(l) step_bc => load_step%get('boundary_conditions') step_mech => step_bc%get('mechanical') @@ -299,7 +301,7 @@ program DAMASK_grid ! doing initialization depending on active solvers call spectral_Utilities_init do field = 1, nActiveFields - select case (loadCases(1)%ID(field)) + select case (ID(field)) case(FIELD_MECH_ID) call mechanical_init @@ -375,7 +377,7 @@ program DAMASK_grid !-------------------------------------------------------------------------------------------------- ! forward fields do field = 1, nActiveFields - select case(loadCases(l)%ID(field)) + select case(ID(field)) case(FIELD_MECH_ID) call mechanical_forward (& cutBack,guess,timeinc,timeIncOld,remainingLoadCaseTime, & @@ -395,7 +397,7 @@ program DAMASK_grid stagIterate = .true. do while (stagIterate) do field = 1, nActiveFields - select case(loadCases(l)%ID(field)) + select case(ID(field)) case(FIELD_MECH_ID) solres(field) = mechanical_solution(incInfo) case(FIELD_THERMAL_ID) From ef543a5b4970d7fd49b26a3bed72939bed393936 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 28 Feb 2021 20:06:31 +0100 Subject: [PATCH 03/13] thermal solver is defined in load case not in numerics.yaml --- src/CPFEM.f90 | 10 +++++----- src/grid/DAMASK_grid.f90 | 4 ++-- src/homogenization.f90 | 14 ++++++-------- 3 files changed, 13 insertions(+), 15 deletions(-) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index f0e510433..01b4c034d 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -178,11 +178,11 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS if (iand(mode, CPFEM_AGERESULTS) /= 0_pInt) call CPFEM_forward - chosenThermal1: select case (thermal_type(material_homogenizationAt(elCP))) - !case (THERMAL_conduction_ID) chosenThermal1 - ! temperature(material_homogenizationAt(elCP))%p(material_homogenizationMemberAt(ip,elCP)) = & - ! temperature_inp - end select chosenThermal1 + !chosenThermal1: select case (thermal_type(material_homogenizationAt(elCP))) + ! case (THERMAL_conduction_ID) chosenThermal1 + ! temperature(material_homogenizationAt(elCP))%p(material_homogenizationMemberAt(ip,elCP)) = & + ! temperature_inp + !end select chosenThermal1 homogenization_F0(1:3,1:3,ma) = ffn homogenization_F(1:3,1:3,ma) = ffn1 diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index 36ebf51a3..7155a5c80 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -163,14 +163,14 @@ program DAMASK_grid !-------------------------------------------------------------------------------------------------- ! initialize field solver information - if (any(thermal_type == THERMAL_conduction_ID )) nActiveFields = nActiveFields + 1 + if (solver%get_asString('thermal',defaultVal = 'n/a') == 'spectral') nActiveFields = nActiveFields + 1 if (any(damage_type == DAMAGE_nonlocal_ID )) nActiveFields = nActiveFields + 1 allocate(solres(nActiveFields)) allocate(ID(nActiveFields)) field = 1 ID(field) = FIELD_MECH_ID ! mechanical active by default - thermalActive: if (any(thermal_type == THERMAL_conduction_ID)) then + thermalActive: if (solver%get_asString('thermal',defaultVal = 'n/a') == 'spectral') then field = field + 1 ID(field) = FIELD_THERMAL_ID endif thermalActive diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 52e90a857..7c0c1531b 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -39,11 +39,11 @@ module homogenization real(pReal), dimension(:), allocatable, public, protected :: & thermal_initialT - integer(kind(THERMAL_isothermal_ID)), dimension(:), allocatable, public, protected :: & + integer(kind(THERMAL_isothermal_ID)), dimension(:), allocatable :: & thermal_type !< thermal transport model integer(kind(DAMAGE_none_ID)), dimension(:), allocatable, public, protected :: & damage_type !< nonlocal damage model - integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:), allocatable, public, protected :: & + integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:), allocatable :: & homogenization_type !< type of each homogenization type, private :: tNumerics_damage @@ -553,10 +553,10 @@ subroutine material_parseHomogenization material_homogenization => config_material%get('homogenization') - allocate(homogenization_type(size(material_name_homogenization)), source=HOMOGENIZATION_undefined_ID) - allocate(thermal_type(size(material_name_homogenization)), source=THERMAL_isothermal_ID) - allocate(damage_type (size(material_name_homogenization)), source=DAMAGE_none_ID) - allocate(thermal_initialT(size(material_name_homogenization)), source=300.0_pReal) + allocate(homogenization_type(size(material_name_homogenization)), source=HOMOGENIZATION_undefined_ID) + allocate(thermal_type(size(material_name_homogenization)), source=THERMAL_isothermal_ID) + allocate(damage_type (size(material_name_homogenization)), source=DAMAGE_none_ID) + allocate(thermal_initialT(size(material_name_homogenization)), source=300.0_pReal) do h=1, size(material_name_homogenization) homog => material_homogenization%get(h) @@ -578,8 +578,6 @@ subroutine material_parseHomogenization thermal_initialT(h) = homogThermal%get_asFloat('T_0',defaultVal=300.0_pReal) select case (homogThermal%get_asString('type')) - case('isothermal') - thermal_type(h) = THERMAL_isothermal_ID case('conduction') thermal_type(h) = THERMAL_conduction_ID case default From ae57ba9707c774ad133a22ef35ec702073dd4486 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 28 Feb 2021 20:10:12 +0100 Subject: [PATCH 04/13] got random errors 'O not found' might be related to use of pointers in openMP. Not nice, but openMP during initialization is not really required --- src/phase_mechanical.f90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index be14c57e6..c934fed72 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -279,7 +279,6 @@ module subroutine mechanical_init(materials,phases) endif - !$OMP PARALLEL DO PRIVATE(ph,me,material,constituents,constituent) do el = 1, size(material_phaseMemberAt,3); do ip = 1, size(material_phaseMemberAt,2) do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) material => materials%get(discretization_materialAt(el)) @@ -305,7 +304,6 @@ module subroutine mechanical_init(materials,phases) enddo enddo; enddo - !$OMP END PARALLEL DO ! initialize plasticity From 0cde43198f5fdc89bc817dc5243c8cd9d7da669b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 28 Feb 2021 20:49:51 +0100 Subject: [PATCH 05/13] modernizing: - 'pass' for dummy thermal homogenization - setting temperature in load case --- PRIVATE | 2 +- src/grid/DAMASK_grid.f90 | 14 ++++++++++++-- src/grid/grid_thermal_spectral.f90 | 4 +++- src/homogenization.f90 | 2 +- 4 files changed, 17 insertions(+), 5 deletions(-) diff --git a/PRIVATE b/PRIVATE index 738526ac7..23ef2dcf3 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 738526ac7ae63678c885faba7621782499f2bce9 +Subproject commit 23ef2dcf34729c8b05caba62f648149ca37a68c9 diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index 7155a5c80..0461b43d9 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -54,6 +54,7 @@ program DAMASK_grid integer, parameter :: & subStepFactor = 2 !< for each substep, divide the last time increment by 2.0 real(pReal) :: & + T_0 = 273.0_pReal, & time = 0.0_pReal, & !< elapsed time time0 = 0.0_pReal, & !< begin of interval timeinc = 1.0_pReal, & !< current time interval @@ -102,6 +103,9 @@ program DAMASK_grid load_steps, & load_step, & solver, & + initial_conditions, & + ic_thermal, & + thermal, & step_bc, & step_mech, & step_discretization, & @@ -181,7 +185,6 @@ program DAMASK_grid !-------------------------------------------------------------------------------------------------- - load_steps => config_load%get('loadstep') allocate(loadCases(load_steps%length)) ! array of load cases @@ -306,7 +309,14 @@ program DAMASK_grid call mechanical_init case(FIELD_THERMAL_ID) - call grid_thermal_spectral_init + if (solver%contains('initial_conditions')) then + initial_conditions => solver%get('initial_conditions') + if (initial_conditions%contains('thermal')) then + thermal => solver%get('thermal') + T_0 = thermal%get_asFloat('thermal',defaultVal = T_0) + endif + endif + call grid_thermal_spectral_init(T_0) case(FIELD_DAMAGE_ID) call grid_damage_spectral_init diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 6403ee955..40a84a02e 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -61,7 +61,9 @@ contains !> @brief allocates all neccessary fields and fills them with data ! ToDo: Restart not implemented !-------------------------------------------------------------------------------------------------- -subroutine grid_thermal_spectral_init +subroutine grid_thermal_spectral_init(T_0) + + real(pReal), intent(in) :: T_0 PetscInt, dimension(0:worldsize-1) :: localK integer :: i, j, k, ce diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 7c0c1531b..1572a0197 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -578,7 +578,7 @@ subroutine material_parseHomogenization thermal_initialT(h) = homogThermal%get_asFloat('T_0',defaultVal=300.0_pReal) select case (homogThermal%get_asString('type')) - case('conduction') + case('pass') thermal_type(h) = THERMAL_conduction_ID case default call IO_error(500,ext_msg=homogThermal%get_asString('type')) From 8af0c8dbc30aa067a94de19427b454cd9dba90b9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 28 Feb 2021 21:51:07 +0100 Subject: [PATCH 06/13] using initial temperature from load case file --- PRIVATE | 2 +- src/grid/DAMASK_grid.f90 | 10 +++++----- src/grid/grid_thermal_spectral.f90 | 3 ++- src/homogenization.f90 | 6 ------ src/homogenization_thermal.f90 | 2 +- 5 files changed, 9 insertions(+), 14 deletions(-) diff --git a/PRIVATE b/PRIVATE index 23ef2dcf3..89eb15f1d 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 23ef2dcf34729c8b05caba62f648149ca37a68c9 +Subproject commit 89eb15f1d771279be8ff444202c32b6797e77bd7 diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index 0461b43d9..3930e0002 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -54,7 +54,7 @@ program DAMASK_grid integer, parameter :: & subStepFactor = 2 !< for each substep, divide the last time increment by 2.0 real(pReal) :: & - T_0 = 273.0_pReal, & + T_0 = 300.0_pReal, & time = 0.0_pReal, & !< elapsed time time0 = 0.0_pReal, & !< begin of interval timeinc = 1.0_pReal, & !< current time interval @@ -309,11 +309,11 @@ program DAMASK_grid call mechanical_init case(FIELD_THERMAL_ID) - if (solver%contains('initial_conditions')) then - initial_conditions => solver%get('initial_conditions') + if (config_load%contains('initial_conditions')) then + initial_conditions => config_load%get('initial_conditions') if (initial_conditions%contains('thermal')) then - thermal => solver%get('thermal') - T_0 = thermal%get_asFloat('thermal',defaultVal = T_0) + thermal => initial_conditions%get('thermal') + T_0 = thermal%get_asFloat('T',defaultVal = T_0) endif endif call grid_thermal_spectral_init(T_0) diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 40a84a02e..5ac043a67 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -133,9 +133,10 @@ subroutine grid_thermal_spectral_init(T_0) ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 - T_current(i,j,k) = homogenization_thermal_T(ce) + T_current(i,j,k) = T_0 T_lastInc(i,j,k) = T_current(i,j,k) T_stagInc(i,j,k) = T_current(i,j,k) + call homogenization_thermal_setField(T_0,0.0_pReal,ce) enddo; enddo; enddo call DMDAVecGetArrayF90(thermal_grid,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with x_scal(xstart:xend,ystart:yend,zstart:zend) = T_current diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 1572a0197..d5fdbd5cb 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -35,10 +35,6 @@ module homogenization homogState, & damageState_h - - real(pReal), dimension(:), allocatable, public, protected :: & - thermal_initialT - integer(kind(THERMAL_isothermal_ID)), dimension(:), allocatable :: & thermal_type !< thermal transport model integer(kind(DAMAGE_none_ID)), dimension(:), allocatable, public, protected :: & @@ -556,7 +552,6 @@ subroutine material_parseHomogenization allocate(homogenization_type(size(material_name_homogenization)), source=HOMOGENIZATION_undefined_ID) allocate(thermal_type(size(material_name_homogenization)), source=THERMAL_isothermal_ID) allocate(damage_type (size(material_name_homogenization)), source=DAMAGE_none_ID) - allocate(thermal_initialT(size(material_name_homogenization)), source=300.0_pReal) do h=1, size(material_name_homogenization) homog => material_homogenization%get(h) @@ -575,7 +570,6 @@ subroutine material_parseHomogenization if(homog%contains('thermal')) then homogThermal => homog%get('thermal') - thermal_initialT(h) = homogThermal%get_asFloat('T_0',defaultVal=300.0_pReal) select case (homogThermal%get_asString('type')) case('pass') diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90 index 928bd547b..efb6e7b58 100644 --- a/src/homogenization_thermal.f90 +++ b/src/homogenization_thermal.f90 @@ -44,7 +44,7 @@ module subroutine thermal_init() allocate(current(configHomogenizations%length)) do ho = 1, configHomogenizations%length - allocate(current(ho)%T(count(material_homogenizationAt2==ho)), source=thermal_initialT(ho)) + allocate(current(ho)%T(count(material_homogenizationAt2==ho)), source=300.0_pReal) allocate(current(ho)%dot_T(count(material_homogenizationAt2==ho)), source=0.0_pReal) configHomogenization => configHomogenizations%get(ho) associate(prm => param(ho)) From 2f68c43755ca3af45dc01701789a732b868adcbd Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 1 Mar 2021 06:16:16 +0100 Subject: [PATCH 07/13] new style configuration for damage --- PRIVATE | 2 +- src/grid/DAMASK_grid.f90 | 4 ++-- src/homogenization.f90 | 6 ++---- src/material.f90 | 4 ++-- src/phase_mechanical_plastic_nonlocal.f90 | 15 ++++++--------- 5 files changed, 13 insertions(+), 18 deletions(-) diff --git a/PRIVATE b/PRIVATE index 89eb15f1d..08e226ec7 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 89eb15f1d771279be8ff444202c32b6797e77bd7 +Subproject commit 08e226ec735589123dd67f3e0a5bc447c3dda812 diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index 3930e0002..fda9d7b9a 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -168,7 +168,7 @@ program DAMASK_grid !-------------------------------------------------------------------------------------------------- ! initialize field solver information if (solver%get_asString('thermal',defaultVal = 'n/a') == 'spectral') nActiveFields = nActiveFields + 1 - if (any(damage_type == DAMAGE_nonlocal_ID )) nActiveFields = nActiveFields + 1 + if (solver%get_asString('damage', defaultVal = 'n/a') == 'spectral') nActiveFields = nActiveFields + 1 allocate(solres(nActiveFields)) allocate(ID(nActiveFields)) @@ -178,7 +178,7 @@ program DAMASK_grid field = field + 1 ID(field) = FIELD_THERMAL_ID endif thermalActive - damageActive: if (any(damage_type == DAMAGE_nonlocal_ID)) then + damageActive: if (solver%get_asString('damage',defaultVal = 'n/a') == 'spectral') then field = field + 1 ID(field) = FIELD_DAMAGE_ID endif damageActive diff --git a/src/homogenization.f90 b/src/homogenization.f90 index d5fdbd5cb..4c7200eed 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -37,7 +37,7 @@ module homogenization integer(kind(THERMAL_isothermal_ID)), dimension(:), allocatable :: & thermal_type !< thermal transport model - integer(kind(DAMAGE_none_ID)), dimension(:), allocatable, public, protected :: & + integer(kind(DAMAGE_none_ID)), dimension(:), allocatable :: & damage_type !< nonlocal damage model integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:), allocatable :: & homogenization_type !< type of each homogenization @@ -582,9 +582,7 @@ subroutine material_parseHomogenization if(homog%contains('damage')) then homogDamage => homog%get('damage') select case (homogDamage%get_asString('type')) - case('none') - damage_type(h) = DAMAGE_none_ID - case('nonlocal') + case('pass') damage_type(h) = DAMAGE_nonlocal_ID case default call IO_error(500,ext_msg=homogDamage%get_asString('type')) diff --git a/src/material.f90 b/src/material.f90 index c7d24b6f0..aecfaa1dd 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -17,7 +17,7 @@ module material private integer, dimension(:), allocatable, public, protected :: & - homogenization_Nconstituents !< number of grains in each homogenization + homogenization_Nconstituents !< number of grains in each homogenization character(len=:), public, protected, allocatable, dimension(:) :: & material_name_phase, & !< name of each phase @@ -30,7 +30,7 @@ module material material_homogenizationAt, & !< homogenization ID of each element material_homogenizationAt2, & !< per cell material_homogenizationMemberAt2 !< cell - integer, dimension(:,:), allocatable, public, protected :: & ! (ip,elem) + integer, dimension(:,:), allocatable :: & ! (ip,elem) material_homogenizationMemberAt !< position of the element within its homogenization instance integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem) material_phaseAt, & !< phase ID of each element diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 601618425..464e57339 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -1822,16 +1822,13 @@ subroutine storeGeometry(ph) integer, intent(in) :: ph integer :: ip, el, ce, co + real(pReal), dimension(:), allocatable :: V + - ce = 0 - do el = 1, size(material_homogenizationMemberAt,2) - do ip = 1, size(material_homogenizationMemberAt,1) - ce = ce + 1 - do co = 1, homogenization_maxNconstituents - if(material_phaseAt2(co,ce) == ph) then - geom(ph)%V_0(material_phaseMemberAt2(co,ce)) = IPvolume(ip,el) - endif - enddo + V = reshape(IPvolume,[product(shape(IPvolume))]) + do ce = 1, size(material_homogenizationMemberAt2,1) + do co = 1, homogenization_maxNconstituents + if(material_phaseAt2(co,ce) == ph) geom(ph)%V_0(material_phaseMemberAt2(co,ce)) = V(ce) enddo enddo From 25c556558065ac0d6508dcba52405afcceeecb70 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 1 Mar 2021 09:24:39 +0100 Subject: [PATCH 08/13] correct runtime test --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 08e226ec7..13dfa0ee9 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 08e226ec735589123dd67f3e0a5bc447c3dda812 +Subproject commit 13dfa0ee9d702782f0b7999f3f7fb2384f58d768 From 438167804c90dd4ba364d59b79c7e58bb3734511 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Thu, 4 Mar 2021 15:11:39 -0500 Subject: [PATCH 09/13] if( polishing --- src/YAML_types.f90 | 116 ++++++++++++++++++++++----------------------- 1 file changed, 58 insertions(+), 58 deletions(-) diff --git a/src/YAML_types.f90 b/src/YAML_types.f90 index b71261d9c..036d40755 100644 --- a/src/YAML_types.f90 +++ b/src/YAML_types.f90 @@ -72,7 +72,7 @@ module YAML_types getKey => tNode_getKey_byIndex procedure :: & contains => tNode_contains - + generic :: & get => tNode_get_byIndex, & tNode_get_byKey @@ -157,7 +157,7 @@ module YAML_types emptyDict type(tList), target, public :: & emptyList - + abstract interface recursive function asFormattedString(self,indent) @@ -179,7 +179,7 @@ module YAML_types public :: & YAML_types_init, & - output_asStrings, & !ToDo: Hack for GNU. Remove later + output_asStrings, & !ToDo: Hack for GNU. Remove later assignment(=) contains @@ -207,11 +207,11 @@ subroutine selfTest select type(s1) class is(tScalar) s1 = '1' - if(s1%asInt() /= 1) error stop 'tScalar_asInt' - if(dNeq(s1%asFloat(),1.0_pReal)) error stop 'tScalar_asFloat' + if (s1%asInt() /= 1) error stop 'tScalar_asInt' + if (dNeq(s1%asFloat(),1.0_pReal)) error stop 'tScalar_asFloat' s1 = 'true' - if(.not. s1%asBool()) error stop 'tScalar_asBool' - if(s1%asString() /= 'true') error stop 'tScalar_asString' + if (.not. s1%asBool()) error stop 'tScalar_asBool' + if (s1%asString() /= 'true') error stop 'tScalar_asString' end select block @@ -232,18 +232,18 @@ subroutine selfTest call l1%append(s1) call l1%append(s2) n => l1 - if(any(l1%asInts() /= [2,3])) error stop 'tList_asInts' - if(any(dNeq(l1%asFloats(),[2.0_pReal,3.0_pReal]))) error stop 'tList_asFloats' - if(n%get_asInt(1) /= 2) error stop 'byIndex_asInt' - if(dNeq(n%get_asFloat(2),3.0_pReal)) error stop 'byIndex_asFloat' + if (any(l1%asInts() /= [2,3])) error stop 'tList_asInts' + if (any(dNeq(l1%asFloats(),[2.0_pReal,3.0_pReal]))) error stop 'tList_asFloats' + if (n%get_asInt(1) /= 2) error stop 'byIndex_asInt' + if (dNeq(n%get_asFloat(2),3.0_pReal)) error stop 'byIndex_asFloat' endselect allocate(tList::l2) select type(l2) class is(tList) call l2%append(l1) - if(any(l2%get_asInts(1) /= [2,3])) error stop 'byIndex_asInts' - if(any(dNeq(l2%get_asFloats(1),[2.0_pReal,3.0_pReal]))) error stop 'byIndex_asFloats' + if (any(l2%get_asInts(1) /= [2,3])) error stop 'byIndex_asInts' + if (any(dNeq(l2%get_asFloats(1),[2.0_pReal,3.0_pReal]))) error stop 'byIndex_asFloats' n => l2 end select deallocate(n) @@ -265,10 +265,10 @@ subroutine selfTest call l1%append(s2) n => l1 - if(any(l1%asBools() .neqv. [.true., .false.])) error stop 'tList_asBools' - if(any(l1%asStrings() /= ['true ','False'])) error stop 'tList_asStrings' - if(n%get_asBool(2)) error stop 'byIndex_asBool' - if(n%get_asString(1) /= 'true') error stop 'byIndex_asString' + if (any(l1%asBools() .neqv. [.true., .false.])) error stop 'tList_asBools' + if (any(l1%asStrings() /= ['true ','False'])) error stop 'tList_asStrings' + if (n%get_asBool(2)) error stop 'byIndex_asBool' + if (n%get_asString(1) /= 'true') error stop 'byIndex_asString' end block end subroutine selfTest @@ -418,7 +418,7 @@ function tNode_get_byIndex(self,i) result(node) integer :: j self_ => self%asList() - if(i < 1 .or. i > self_%length) call IO_error(150,ext_msg='tNode_get_byIndex') + if (i < 1 .or. i > self_%length) call IO_error(150,ext_msg='tNode_get_byIndex') j = 1 item => self_%first @@ -599,7 +599,7 @@ function tNode_getKey_byIndex(self,i) result(key) dict => self%asDict() item => dict%first do j = 1, dict%length - if(j == i) then + if (j == i) then key = item%key exit else @@ -613,7 +613,7 @@ end function tNode_getKey_byIndex !------------------------------------------------------------------------------------------------- !> @brief Checks if a given key/item is present in the dict/list !------------------------------------------------------------------------------------------------- -function tNode_contains(self,k) result(exists) +function tNode_contains(self,k) result(exists) class(tNode), intent(in), target :: self character(len=*), intent(in) :: k @@ -624,18 +624,18 @@ function tNode_contains(self,k) result(exists) type(tDict), pointer :: dict exists = .false. - if(self%isDict()) then + if (self%isDict()) then dict => self%asDict() do j=1, dict%length - if(dict%getKey(j) == k) then + if (dict%getKey(j) == k) then exists = .true. return endif enddo - elseif(self%isList()) then + elseif (self%isList()) then list => self%asList() - do j =1, list%length - if(list%get_asString(j) == k) then + do j=1, list%length + if (list%get_asString(j) == k) then exists = .true. return endif @@ -663,8 +663,8 @@ function tNode_get_byKey(self,k,defaultVal) result(node) logical :: found found = present(defaultVal) - if(found) node => defaultVal - + if (found) node => defaultVal + self_ => self%asDict() j = 1 @@ -677,11 +677,11 @@ function tNode_get_byKey(self,k,defaultVal) result(node) item => item%next j = j + 1 enddo - + if (.not. found) then call IO_error(143,ext_msg=k) else - if(associated(item)) node => item%node + if (associated(item)) node => item%node endif end function tNode_get_byKey @@ -700,11 +700,11 @@ function tNode_get_byKey_asFloat(self,k,defaultVal) result(nodeAsFloat) class(tNode), pointer :: node type(tScalar), pointer :: scalar - if(self%contains(k)) then + if (self%contains(k)) then node => self%get(k) scalar => node%asScalar() nodeAsFloat = scalar%asFloat() - elseif(present(defaultVal)) then + elseif (present(defaultVal)) then nodeAsFloat = defaultVal else call IO_error(143,ext_msg=k) @@ -726,11 +726,11 @@ function tNode_get_byKey_asInt(self,k,defaultVal) result(nodeAsInt) class(tNode), pointer :: node type(tScalar), pointer :: scalar - if(self%contains(k)) then + if (self%contains(k)) then node => self%get(k) scalar => node%asScalar() nodeAsInt = scalar%asInt() - elseif(present(defaultVal)) then + elseif (present(defaultVal)) then nodeAsInt = defaultVal else call IO_error(143,ext_msg=k) @@ -752,11 +752,11 @@ function tNode_get_byKey_asBool(self,k,defaultVal) result(nodeAsBool) class(tNode), pointer :: node type(tScalar), pointer :: scalar - if(self%contains(k)) then + if (self%contains(k)) then node => self%get(k) scalar => node%asScalar() nodeAsBool = scalar%asBool() - elseif(present(defaultVal)) then + elseif (present(defaultVal)) then nodeAsBool = defaultVal else call IO_error(143,ext_msg=k) @@ -778,11 +778,11 @@ function tNode_get_byKey_asString(self,k,defaultVal) result(nodeAsString) class(tNode), pointer :: node type(tScalar), pointer :: scalar - if(self%contains(k)) then + if (self%contains(k)) then node => self%get(k) scalar => node%asScalar() nodeAsString = scalar%asString() - elseif(present(defaultVal)) then + elseif (present(defaultVal)) then nodeAsString = defaultVal else call IO_error(143,ext_msg=k) @@ -806,18 +806,18 @@ function tNode_get_byKey_asFloats(self,k,defaultVal,requiredSize) result(nodeAsF class(tNode), pointer :: node type(tList), pointer :: list - if(self%contains(k)) then + if (self%contains(k)) then node => self%get(k) list => node%asList() nodeAsFloats = list%asFloats() - elseif(present(defaultVal)) then + elseif (present(defaultVal)) then nodeAsFloats = defaultVal else call IO_error(143,ext_msg=k) endif - if(present(requiredSize)) then - if(requiredSize /= size(nodeAsFloats)) call IO_error(146,ext_msg=k) + if (present(requiredSize)) then + if (requiredSize /= size(nodeAsFloats)) call IO_error(146,ext_msg=k) endif end function tNode_get_byKey_asFloats @@ -837,18 +837,18 @@ function tNode_get_byKey_asInts(self,k,defaultVal,requiredSize) result(nodeAsInt class(tNode), pointer :: node type(tList), pointer :: list - if(self%contains(k)) then + if (self%contains(k)) then node => self%get(k) list => node%asList() nodeAsInts = list%asInts() - elseif(present(defaultVal)) then + elseif (present(defaultVal)) then nodeAsInts = defaultVal else call IO_error(143,ext_msg=k) endif - if(present(requiredSize)) then - if(requiredSize /= size(nodeAsInts)) call IO_error(146,ext_msg=k) + if (present(requiredSize)) then + if (requiredSize /= size(nodeAsInts)) call IO_error(146,ext_msg=k) endif end function tNode_get_byKey_asInts @@ -867,11 +867,11 @@ function tNode_get_byKey_asBools(self,k,defaultVal) result(nodeAsBools) class(tNode), pointer :: node type(tList), pointer :: list - if(self%contains(k)) then + if (self%contains(k)) then node => self%get(k) list => node%asList() nodeAsBools = list%asBools() - elseif(present(defaultVal)) then + elseif (present(defaultVal)) then nodeAsBools = defaultVal else call IO_error(143,ext_msg=k) @@ -893,11 +893,11 @@ function tNode_get_byKey_asStrings(self,k,defaultVal) result(nodeAsStrings) class(tNode), pointer :: node type(tList), pointer :: list - if(self%contains(k)) then + if (self%contains(k)) then node => self%get(k) list => node%asList() nodeAsStrings = list%asStrings() - elseif(present(defaultVal)) then + elseif (present(defaultVal)) then nodeAsStrings = defaultVal else call IO_error(143,ext_msg=k) @@ -925,7 +925,7 @@ function output_asStrings(self) result(output) !ToDo: SR: Rem end function output_asStrings - + !-------------------------------------------------------------------------------------------------- !> @brief Returns the index of a key in a dictionary @@ -944,7 +944,7 @@ function tNode_get_byKey_asIndex(self,key) result(keyIndex) item => dict%first keyIndex = -1 do i = 1, dict%length - if(key == item%key) then + if (key == item%key) then keyIndex = i exit else @@ -952,9 +952,9 @@ function tNode_get_byKey_asIndex(self,key) result(keyIndex) endif enddo - if(keyIndex == -1) call IO_error(140,ext_msg=key) + if (keyIndex == -1) call IO_error(140,ext_msg=key) + - end function tNode_get_byKey_asIndex @@ -985,7 +985,7 @@ recursive function tList_asFormattedString(self,indent) result(str) integer :: i, indent_ str = '' - if(present(indent)) then + if (present(indent)) then indent_ = indent else indent_ = 0 @@ -993,7 +993,7 @@ recursive function tList_asFormattedString(self,indent) result(str) item => self%first do i = 1, self%length - if(i /= 1) str = str//repeat(' ',indent_) + if (i /= 1) str = str//repeat(' ',indent_) str = str//'- '//item%node%asFormattedString(indent_+2) item => item%next end do @@ -1014,7 +1014,7 @@ recursive function tDict_asFormattedString(self,indent) result(str) integer :: i, indent_ str = '' - if(present(indent)) then + if (present(indent)) then indent_ = indent else indent_ = 0 @@ -1022,7 +1022,7 @@ recursive function tDict_asFormattedString(self,indent) result(str) item => self%first do i = 1, self%length - if(i /= 1) str = str//repeat(' ',indent_) + if (i /= 1) str = str//repeat(' ',indent_) select type(node_1 =>item%node) class is(tScalar) str = str//trim(item%key)//': '//item%node%asFormattedString(indent_+len_trim(item%key)+2) @@ -1270,7 +1270,7 @@ recursive subroutine tItem_finalize(self) type(tItem),intent(inout) :: self deallocate(self%node) - if(associated(self%next)) deallocate(self%next) + if (associated(self%next)) deallocate(self%next) end subroutine tItem_finalize From e4271537c5640f4ddb1645c6d2a2f20ede059337 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Thu, 4 Mar 2021 15:14:16 -0500 Subject: [PATCH 10/13] syntax polishing; use of YAML defaults to avoid if%contains --- src/grid/DAMASK_grid.f90 | 37 ++++++++++++++++--------------------- 1 file changed, 16 insertions(+), 21 deletions(-) diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index fda9d7b9a..0aa462e7e 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -38,7 +38,7 @@ program DAMASK_grid f_restart !< frequency of restart writes logical :: estimate_rate !< follow trajectory of former loadcase end type tLoadCase - + integer(kind(FIELD_UNDEFINED_ID)), allocatable :: ID(:) !-------------------------------------------------------------------------------------------------- @@ -116,7 +116,7 @@ program DAMASK_grid ! init DAMASK (all modules) call CPFEM_initAll - print'(/,a)', ' <<<+- DAMASK_spectral init -+>>>'; flush(IO_STDOUT) + print'(/,a)', ' <<<+- DAMASK_grid init -+>>>'; flush(IO_STDOUT) print*, 'Shanthraj et al., Handbook of Mechanics of Materials, 2019' print*, 'https://doi.org/10.1007/978-981-10-6855-3_80' @@ -130,7 +130,7 @@ program DAMASK_grid if (stagItMax < 0) call IO_error(301,ext_msg='maxStaggeredIter') if (maxCutBack < 0) call IO_error(301,ext_msg='maxCutBack') - + config_load => YAML_parse_file(trim(interface_loadFile)) solver => config_load%get('solver') @@ -169,9 +169,10 @@ program DAMASK_grid ! initialize field solver information if (solver%get_asString('thermal',defaultVal = 'n/a') == 'spectral') nActiveFields = nActiveFields + 1 if (solver%get_asString('damage', defaultVal = 'n/a') == 'spectral') nActiveFields = nActiveFields + 1 + allocate(solres(nActiveFields)) - - allocate(ID(nActiveFields)) + allocate( ID(nActiveFields)) + field = 1 ID(field) = FIELD_MECH_ID ! mechanical active by default thermalActive: if (solver%get_asString('thermal',defaultVal = 'n/a') == 'spectral') then @@ -224,16 +225,15 @@ program DAMASK_grid if (.not. allocated(loadCases(l)%deformation%myType)) call IO_error(error_ID=837,ext_msg = 'L/dot_F/F missing') step_discretization => load_step%get('discretization') - if(.not. step_discretization%contains('t')) call IO_error(error_ID=837,ext_msg = 't missing') - if(.not. step_discretization%contains('N')) call IO_error(error_ID=837,ext_msg = 'N missing') + if (.not. step_discretization%contains('t')) call IO_error(error_ID=837,ext_msg = 't missing') + if (.not. step_discretization%contains('N')) call IO_error(error_ID=837,ext_msg = 'N missing') loadCases(l)%t = step_discretization%get_asFloat('t') loadCases(l)%N = step_discretization%get_asInt ('N') loadCases(l)%r = step_discretization%get_asFloat('r', defaultVal= 1.0_pReal) loadCases(l)%f_restart = load_step%get_asInt('f_restart', defaultVal=huge(0)) loadCases(l)%f_out = load_step%get_asInt('f_out', defaultVal=1) - loadCases(l)%estimate_rate = (load_step%get_asBool('estimate_rate',defaultVal=.true.) .and. & - merge(.true.,.false.,l > 1)) + loadCases(l)%estimate_rate = (load_step%get_asBool('estimate_rate',defaultVal=.true.) .and. l>1) reportAndCheck: if (worldrank == 0) then print'(/,a,i0)', ' load case: ', l @@ -289,11 +289,11 @@ program DAMASK_grid else print'(a,f0.3)', ' r: ', loadCases(l)%r endif - print'(a,f0.3)', ' t: ', loadCases(l)%t - print'(a,i0)', ' N: ', loadCases(l)%N - print'(a,i0)', ' f_out: ', loadCases(l)%f_out + print'(a,f0.3)', ' t: ', loadCases(l)%t + print'(a,i0)', ' N: ', loadCases(l)%N + print'(a,i0)', ' f_out: ', loadCases(l)%f_out if (loadCases(l)%f_restart < huge(0)) & - print'(a,i0)', ' f_restart: ', loadCases(l)%f_restart + print'(a,i0)', ' f_restart: ', loadCases(l)%f_restart if (errorID > 0) call IO_error(error_ID = errorID, el = l) @@ -309,14 +309,9 @@ program DAMASK_grid call mechanical_init case(FIELD_THERMAL_ID) - if (config_load%contains('initial_conditions')) then - initial_conditions => config_load%get('initial_conditions') - if (initial_conditions%contains('thermal')) then - thermal => initial_conditions%get('thermal') - T_0 = thermal%get_asFloat('T',defaultVal = T_0) - endif - endif - call grid_thermal_spectral_init(T_0) + initial_conditions => config_load%get('initial_conditions',defaultVal=emptyDict) + thermal => initial_conditions%get('thermal',defaultVal=emptyDict) + call grid_thermal_spectral_init(thermal%get_asFloat('T',defaultVal = T_0)) case(FIELD_DAMAGE_ID) call grid_damage_spectral_init From e1f0d2e0a3a99ad4488cee5a673ef3d13d5017b8 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Thu, 4 Mar 2021 15:15:40 -0500 Subject: [PATCH 11/13] polishing of indentation and whitespaces; thermal_homogenize only once after all constituents --- src/homogenization.f90 | 41 ++++++++++++++--------------------------- 1 file changed, 14 insertions(+), 27 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 4c7200eed..90654ee99 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -91,7 +91,7 @@ module homogenization real(pReal), intent(in), dimension(3,3) :: & subF integer, intent(in) :: & - ce + ce end subroutine mechanical_partition module subroutine thermal_partition(ce) @@ -129,10 +129,8 @@ module homogenization module function thermal_conduction_getConductivity(ce) result(K) - integer, intent(in) :: ce real(pReal), dimension(3,3) :: K - end function thermal_conduction_getConductivity module function thermal_conduction_getSpecificHeat(ce) result(c_P) @@ -167,13 +165,12 @@ module homogenization real(pReal), intent(out) :: Tdot end subroutine thermal_conduction_getSource - module function damage_nonlocal_getMobility(ce) result(M) - integer, intent(in) :: ce - real(pReal) :: M + module function damage_nonlocal_getMobility(ce) result(M) + integer, intent(in) :: ce + real(pReal) :: M end function damage_nonlocal_getMobility module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ce) - integer, intent(in) :: ce real(pReal), intent(in) :: & phi @@ -181,21 +178,17 @@ module homogenization phiDot, dPhiDot_dPhi end subroutine damage_nonlocal_getSourceAndItsTangent - module subroutine damage_nonlocal_putNonLocalDamage(phi,ce) - integer, intent(in) :: ce real(pReal), intent(in) :: & phi - end subroutine damage_nonlocal_putNonLocalDamage module subroutine damage_nonlocal_results(ho,group) - integer, intent(in) :: ho character(len=*), intent(in) :: group - end subroutine damage_nonlocal_results + end interface public :: & @@ -238,21 +231,18 @@ subroutine homogenization_init() allocate(homogState (size(material_name_homogenization))) allocate(damageState_h (size(material_name_homogenization))) - call material_parseHomogenization - + call material_parseHomogenization() num_homog => config_numerics%get('homogenization',defaultVal=emptyDict) num_homogGeneric => num_homog%get('generic',defaultVal=emptyDict) - num%nMPstate = num_homogGeneric%get_asInt ('nMPstate', defaultVal=10) - if (num%nMPstate < 1) call IO_error(301,ext_msg='nMPstate') - + num%nMPstate = num_homogGeneric%get_asInt('nMPstate',defaultVal=10) + if (num%nMPstate < 1) call IO_error(301,ext_msg='nMPstate') call mechanical_init(num_homog) call thermal_init() call damage_init() - - call damage_nonlocal_init + call damage_nonlocal_init() end subroutine homogenization_init @@ -322,7 +312,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE enddo !$OMP END DO - if (.not. terminallyIll ) then + if (.not. terminallyIll) then !$OMP DO PRIVATE(ho,ph,ce) do el = FEsolving_execElem(1),FEsolving_execElem(2) if (terminallyIll) continue @@ -336,9 +326,9 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE if (.not. terminallyIll) & ! so first signals terminally ill... print*, ' Integration point ', ip,' at element ', el, ' terminally ill' terminallyIll = .true. ! ...and kills all others - endif - call thermal_homogenize(ip,el) + endif enddo + call thermal_homogenize(ip,el) enddo enddo !$OMP END DO @@ -567,10 +557,8 @@ subroutine material_parseHomogenization call IO_error(500,ext_msg=homogMech%get_asString('type')) end select - - if(homog%contains('thermal')) then + if (homog%contains('thermal')) then homogThermal => homog%get('thermal') - select case (homogThermal%get_asString('type')) case('pass') thermal_type(h) = THERMAL_conduction_ID @@ -579,7 +567,7 @@ subroutine material_parseHomogenization end select endif - if(homog%contains('damage')) then + if (homog%contains('damage')) then homogDamage => homog%get('damage') select case (homogDamage%get_asString('type')) case('pass') @@ -590,7 +578,6 @@ subroutine material_parseHomogenization endif enddo - end subroutine material_parseHomogenization From a7e2ed40dd183444e7500a6899ddde33ee188d98 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Thu, 4 Mar 2021 15:16:36 -0500 Subject: [PATCH 12/13] rename Nconstituents --> Nmembers --- src/phase_damage.f90 | 10 ++-- src/phase_damage_anisobrittle.f90 | 6 +- src/phase_damage_anisoductile.f90 | 6 +- src/phase_damage_isobrittle.f90 | 6 +- src/phase_damage_isoductile.f90 | 6 +- src/phase_mechanical.f90 | 38 ++++++------- ...phase_mechanical_plastic_dislotungsten.f90 | 14 ++--- src/phase_mechanical_plastic_dislotwin.f90 | 30 +++++----- src/phase_mechanical_plastic_isotropic.f90 | 6 +- ...phase_mechanical_plastic_kinehardening.f90 | 8 +-- src/phase_mechanical_plastic_nonlocal.f90 | 54 +++++++++--------- ...phase_mechanical_plastic_phenopowerlaw.f90 | 10 ++-- src/phase_thermal.f90 | 56 ++++++++----------- src/phase_thermal_dissipation.f90 | 6 +- src/phase_thermal_externalheat.f90 | 6 +- 15 files changed, 127 insertions(+), 135 deletions(-) diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index c6fd7bfab..6ad25f1e9 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -133,7 +133,7 @@ module subroutine damage_init integer :: & ph, & !< counter in phase loop - Nconstituents + Nmembers class(tNode), pointer :: & phases, & phase, & @@ -151,10 +151,10 @@ module subroutine damage_init do ph = 1,phases%length - Nconstituents = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseAt2 == ph) - allocate(current(ph)%phi(Nconstituents),source=1.0_pReal) - allocate(current(ph)%d_phi_d_dot_phi(Nconstituents),source=0.0_pReal) + allocate(current(ph)%phi(Nmembers),source=1.0_pReal) + allocate(current(ph)%d_phi_d_dot_phi(Nmembers),source=0.0_pReal) phase => phases%get(ph) sources => phase%get('damage',defaultVal=emptyList) @@ -199,7 +199,7 @@ module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, phiDot = 0.0_pReal dPhiDot_dPhi = 0.0_pReal - do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) + do co = 1, homogenization_Nmembers(material_homogenizationAt2(ce)) ph = material_phaseAt2(co,ce) me = material_phasememberAt2(co,ce) diff --git a/src/phase_damage_anisobrittle.f90 b/src/phase_damage_anisobrittle.f90 index fd82b74c4..096da6fb8 100644 --- a/src/phase_damage_anisobrittle.f90 +++ b/src/phase_damage_anisobrittle.f90 @@ -40,7 +40,7 @@ module function anisobrittle_init() result(mySources) phase, & sources, & src - integer :: Nconstituents,p + integer :: Nmembers,p integer, dimension(:), allocatable :: N_cl character(len=pStringLen) :: extmsg = '' @@ -92,8 +92,8 @@ module function anisobrittle_init() result(mySources) if (any(prm%g_crit < 0.0_pReal)) extmsg = trim(extmsg)//' g_crit' if (any(prm%s_crit < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit' - Nconstituents = count(material_phaseAt==p) * discretization_nIPs - call phase_allocateState(damageState(p),Nconstituents,1,1,0) + Nmembers = count(material_phaseAt==p) * discretization_nIPs + call phase_allocateState(damageState(p),Nmembers,1,1,0) damageState(p)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal) if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol' diff --git a/src/phase_damage_anisoductile.f90 b/src/phase_damage_anisoductile.f90 index f50c05f07..a687df594 100644 --- a/src/phase_damage_anisoductile.f90 +++ b/src/phase_damage_anisoductile.f90 @@ -35,7 +35,7 @@ module function anisoductile_init() result(mySources) pl, & sources, & src - integer :: Ninstances,Nconstituents,p + integer :: Ninstances,Nmembers,p integer, dimension(:), allocatable :: N_sl character(len=pStringLen) :: extmsg = '' @@ -78,8 +78,8 @@ module function anisoductile_init() result(mySources) if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit' - Nconstituents=count(material_phaseAt2==p) - call phase_allocateState(damageState(p),Nconstituents,1,1,0) + Nmembers=count(material_phaseAt2==p) + call phase_allocateState(damageState(p),Nmembers,1,1,0) damageState(p)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal) if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol' diff --git a/src/phase_damage_isobrittle.f90 b/src/phase_damage_isobrittle.f90 index 1cc0d7900..0cf85897a 100644 --- a/src/phase_damage_isobrittle.f90 +++ b/src/phase_damage_isobrittle.f90 @@ -31,7 +31,7 @@ module function isobrittle_init() result(mySources) phase, & sources, & src - integer :: Nconstituents,p + integer :: Nmembers,p character(len=pStringLen) :: extmsg = '' @@ -64,8 +64,8 @@ module function isobrittle_init() result(mySources) ! sanity checks if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit' - Nconstituents = count(material_phaseAt2==p) - call phase_allocateState(damageState(p),Nconstituents,1,1,1) + Nmembers = count(material_phaseAt2==p) + call phase_allocateState(damageState(p),Nmembers,1,1,1) damageState(p)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal) if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' diff --git a/src/phase_damage_isoductile.f90 b/src/phase_damage_isoductile.f90 index 247cd715a..9d00bb1a7 100644 --- a/src/phase_damage_isoductile.f90 +++ b/src/phase_damage_isoductile.f90 @@ -33,7 +33,7 @@ module function isoductile_init() result(mySources) phase, & sources, & src - integer :: Ninstances,Nconstituents,p + integer :: Ninstances,Nmembers,p character(len=pStringLen) :: extmsg = '' @@ -68,8 +68,8 @@ module function isoductile_init() result(mySources) if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit' - Nconstituents=count(material_phaseAt2==p) - call phase_allocateState(damageState(p),Nconstituents,1,1,0) + Nmembers=count(material_phaseAt2==p) + call phase_allocateState(damageState(p),Nmembers,1,1,0) damageState(p)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal) if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index c934fed72..d065fc74a 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -191,7 +191,7 @@ module subroutine mechanical_init(materials,phases) ph, & me, & stiffDegradationCtr, & - Nconstituents + Nmembers class(tNode), pointer :: & num_crystallite, & material, & @@ -229,22 +229,22 @@ module subroutine mechanical_init(materials,phases) allocate(material_orientation0(homogenization_maxNconstituents,phases%length,maxVal(material_phaseMemberAt))) do ph = 1, phases%length - Nconstituents = count(material_phaseAt == ph) * discretization_nIPs + Nmembers = count(material_phaseAt == ph) * discretization_nIPs - allocate(phase_mechanical_Fi(ph)%data(3,3,Nconstituents)) - allocate(phase_mechanical_Fe(ph)%data(3,3,Nconstituents)) - allocate(phase_mechanical_Fi0(ph)%data(3,3,Nconstituents)) - allocate(phase_mechanical_Fp(ph)%data(3,3,Nconstituents)) - allocate(phase_mechanical_Fp0(ph)%data(3,3,Nconstituents)) - allocate(phase_mechanical_Li(ph)%data(3,3,Nconstituents)) - allocate(phase_mechanical_Li0(ph)%data(3,3,Nconstituents)) - allocate(phase_mechanical_Lp0(ph)%data(3,3,Nconstituents)) - allocate(phase_mechanical_Lp(ph)%data(3,3,Nconstituents)) - allocate(phase_mechanical_S(ph)%data(3,3,Nconstituents),source=0.0_pReal) - allocate(phase_mechanical_P(ph)%data(3,3,Nconstituents),source=0.0_pReal) - allocate(phase_mechanical_S0(ph)%data(3,3,Nconstituents),source=0.0_pReal) - allocate(phase_mechanical_F(ph)%data(3,3,Nconstituents)) - allocate(phase_mechanical_F0(ph)%data(3,3,Nconstituents)) + allocate(phase_mechanical_Fi(ph)%data(3,3,Nmembers)) + allocate(phase_mechanical_Fe(ph)%data(3,3,Nmembers)) + allocate(phase_mechanical_Fi0(ph)%data(3,3,Nmembers)) + allocate(phase_mechanical_Fp(ph)%data(3,3,Nmembers)) + allocate(phase_mechanical_Fp0(ph)%data(3,3,Nmembers)) + allocate(phase_mechanical_Li(ph)%data(3,3,Nmembers)) + allocate(phase_mechanical_Li0(ph)%data(3,3,Nmembers)) + allocate(phase_mechanical_Lp0(ph)%data(3,3,Nmembers)) + allocate(phase_mechanical_Lp(ph)%data(3,3,Nmembers)) + allocate(phase_mechanical_S(ph)%data(3,3,Nmembers),source=0.0_pReal) + allocate(phase_mechanical_P(ph)%data(3,3,Nmembers),source=0.0_pReal) + allocate(phase_mechanical_S0(ph)%data(3,3,Nmembers),source=0.0_pReal) + allocate(phase_mechanical_F(ph)%data(3,3,Nmembers)) + allocate(phase_mechanical_F0(ph)%data(3,3,Nmembers)) phase => phases%get(ph) mech => phase%get('mechanics') @@ -278,9 +278,9 @@ module subroutine mechanical_init(materials,phases) enddo endif - + do el = 1, size(material_phaseMemberAt,3); do ip = 1, size(material_phaseMemberAt,2) - do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) + do co = 1, homogenization_Nmembers(material_homogenizationAt(el)) material => materials%get(discretization_materialAt(el)) constituents => material%get('constituents') constituent => constituents%get(co) @@ -1238,7 +1238,7 @@ module subroutine mechanical_restore(ce,includeL) co, ph, me - do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce)) + do co = 1,homogenization_Nmembers(material_homogenizationAt2(ce)) ph = material_phaseAt2(co,ce) me = material_phaseMemberAt2(co,ce) if (includeL) then diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index 9c8baf0a5..1fda51a58 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -79,7 +79,7 @@ module function plastic_dislotungsten_init() result(myPlasticity) logical, dimension(:), allocatable :: myPlasticity integer :: & ph, i, & - Nconstituents, & + Nmembers, & sizeState, sizeDotState, & startIndex, endIndex integer, dimension(:), allocatable :: & @@ -220,18 +220,18 @@ module function plastic_dislotungsten_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nconstituents = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseAt2 == ph) sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl sizeState = sizeDotState - call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,0) + call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0) !-------------------------------------------------------------------------------------------------- ! state aliases and initialization startIndex = 1 endIndex = prm%sum_N_sl stt%rho_mob => plasticState(ph)%state(startIndex:endIndex,:) - stt%rho_mob = spread(rho_mob_0,2,Nconstituents) + stt%rho_mob = spread(rho_mob_0,2,Nmembers) dot%rho_mob => plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho' @@ -239,7 +239,7 @@ module function plastic_dislotungsten_init() result(myPlasticity) startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl stt%rho_dip => plasticState(ph)%state(startIndex:endIndex,:) - stt%rho_dip = spread(rho_dip_0,2,Nconstituents) + stt%rho_dip = spread(rho_dip_0,2,Nmembers) dot%rho_dip => plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) @@ -251,8 +251,8 @@ module function plastic_dislotungsten_init() result(myPlasticity) ! global alias plasticState(ph)%slipRate => plasticState(ph)%dotState(startIndex:endIndex,:) - allocate(dst%Lambda_sl(prm%sum_N_sl,Nconstituents), source=0.0_pReal) - allocate(dst%threshold_stress(prm%sum_N_sl,Nconstituents), source=0.0_pReal) + allocate(dst%Lambda_sl(prm%sum_N_sl,Nmembers), source=0.0_pReal) + allocate(dst%threshold_stress(prm%sum_N_sl,Nmembers), source=0.0_pReal) plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index f479ba007..5500ae731 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -127,7 +127,7 @@ module function plastic_dislotwin_init() result(myPlasticity) logical, dimension(:), allocatable :: myPlasticity integer :: & ph, i, & - Nconstituents, & + Nmembers, & sizeState, sizeDotState, & startIndex, endIndex integer, dimension(:), allocatable :: & @@ -406,21 +406,21 @@ module function plastic_dislotwin_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nconstituents = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseAt2 == ph) sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl & + size(['f_tw']) * prm%sum_N_tw & + size(['f_tr']) * prm%sum_N_tr sizeState = sizeDotState - call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,0) + call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0) !-------------------------------------------------------------------------------------------------- ! locally defined state aliases and initialization of state0 and atol startIndex = 1 endIndex = prm%sum_N_sl stt%rho_mob=>plasticState(ph)%state(startIndex:endIndex,:) - stt%rho_mob= spread(rho_mob_0,2,Nconstituents) + stt%rho_mob= spread(rho_mob_0,2,Nmembers) dot%rho_mob=>plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho' @@ -428,7 +428,7 @@ module function plastic_dislotwin_init() result(myPlasticity) startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl stt%rho_dip=>plasticState(ph)%state(startIndex:endIndex,:) - stt%rho_dip= spread(rho_dip_0,2,Nconstituents) + stt%rho_dip= spread(rho_dip_0,2,Nmembers) dot%rho_dip=>plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) @@ -454,18 +454,18 @@ module function plastic_dislotwin_init() result(myPlasticity) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_f_tr',defaultVal=1.0e-6_pReal) if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_f_tr' - allocate(dst%Lambda_sl (prm%sum_N_sl,Nconstituents),source=0.0_pReal) - allocate(dst%tau_pass (prm%sum_N_sl,Nconstituents),source=0.0_pReal) + allocate(dst%Lambda_sl (prm%sum_N_sl,Nmembers),source=0.0_pReal) + allocate(dst%tau_pass (prm%sum_N_sl,Nmembers),source=0.0_pReal) - allocate(dst%Lambda_tw (prm%sum_N_tw,Nconstituents),source=0.0_pReal) - allocate(dst%tau_hat_tw (prm%sum_N_tw,Nconstituents),source=0.0_pReal) - allocate(dst%tau_r_tw (prm%sum_N_tw,Nconstituents),source=0.0_pReal) - allocate(dst%V_tw (prm%sum_N_tw,Nconstituents),source=0.0_pReal) + allocate(dst%Lambda_tw (prm%sum_N_tw,Nmembers),source=0.0_pReal) + allocate(dst%tau_hat_tw (prm%sum_N_tw,Nmembers),source=0.0_pReal) + allocate(dst%tau_r_tw (prm%sum_N_tw,Nmembers),source=0.0_pReal) + allocate(dst%V_tw (prm%sum_N_tw,Nmembers),source=0.0_pReal) - allocate(dst%Lambda_tr (prm%sum_N_tr,Nconstituents),source=0.0_pReal) - allocate(dst%tau_hat_tr (prm%sum_N_tr,Nconstituents),source=0.0_pReal) - allocate(dst%tau_r_tr (prm%sum_N_tr,Nconstituents),source=0.0_pReal) - allocate(dst%V_tr (prm%sum_N_tr,Nconstituents),source=0.0_pReal) + allocate(dst%Lambda_tr (prm%sum_N_tr,Nmembers),source=0.0_pReal) + allocate(dst%tau_hat_tr (prm%sum_N_tr,Nmembers),source=0.0_pReal) + allocate(dst%tau_r_tr (prm%sum_N_tr,Nmembers),source=0.0_pReal) + allocate(dst%V_tr (prm%sum_N_tr,Nmembers),source=0.0_pReal) plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally diff --git a/src/phase_mechanical_plastic_isotropic.f90 b/src/phase_mechanical_plastic_isotropic.f90 index 245f4293a..d02436fba 100644 --- a/src/phase_mechanical_plastic_isotropic.f90 +++ b/src/phase_mechanical_plastic_isotropic.f90 @@ -52,7 +52,7 @@ module function plastic_isotropic_init() result(myPlasticity) logical, dimension(:), allocatable :: myPlasticity integer :: & ph, & - Nconstituents, & + Nmembers, & sizeState, sizeDotState real(pReal) :: & xi_0 !< initial critical stress @@ -119,11 +119,11 @@ module function plastic_isotropic_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nconstituents = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseAt2 == ph) sizeDotState = size(['xi ','gamma']) sizeState = sizeDotState - call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,0) + call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0) !-------------------------------------------------------------------------------------------------- ! state aliases and initialization diff --git a/src/phase_mechanical_plastic_kinehardening.f90 b/src/phase_mechanical_plastic_kinehardening.f90 index 8acf42710..75e8d9e59 100644 --- a/src/phase_mechanical_plastic_kinehardening.f90 +++ b/src/phase_mechanical_plastic_kinehardening.f90 @@ -62,7 +62,7 @@ module function plastic_kinehardening_init() result(myPlasticity) logical, dimension(:), allocatable :: myPlasticity integer :: & ph, o, & - Nconstituents, & + Nmembers, & sizeState, sizeDeltaState, sizeDotState, & startIndex, endIndex integer, dimension(:), allocatable :: & @@ -165,19 +165,19 @@ module function plastic_kinehardening_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nconstituents = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseAt2 == ph) sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%sum_N_sl !ToDo: adjust names like in material.yaml sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl !ToDo: adjust names like in material.yaml sizeState = sizeDotState + sizeDeltaState - call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,sizeDeltaState) + call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState) !-------------------------------------------------------------------------------------------------- ! state aliases and initialization startIndex = 1 endIndex = prm%sum_N_sl stt%crss => plasticState(ph)%state (startIndex:endIndex,:) - stt%crss = spread(xi_0, 2, Nconstituents) + stt%crss = spread(xi_0, 2, Nmembers) dot%crss => plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 464e57339..d0007075b 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -177,7 +177,7 @@ module function plastic_nonlocal_init() result(myPlasticity) integer :: & Ninstances, & ph, & - Nconstituents, & + Nmembers, & sizeState, sizeDotState, sizeDependentState, sizeDeltaState, & s1, s2, & s, t, l @@ -398,7 +398,7 @@ module function plastic_nonlocal_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nconstituents = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseAt2 == ph) sizeDotState = size([ 'rhoSglEdgePosMobile ','rhoSglEdgeNegMobile ', & 'rhoSglScrewPosMobile ','rhoSglScrewNegMobile ', & 'rhoSglEdgePosImmobile ','rhoSglEdgeNegImmobile ', & @@ -412,9 +412,9 @@ module function plastic_nonlocal_init() result(myPlasticity) 'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure sizeDeltaState = sizeDotState - call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,sizeDeltaState) + call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState) - allocate(geom(ph)%V_0(Nconstituents)) + allocate(geom(ph)%V_0(Nmembers)) call storeGeometry(ph) plasticState(ph)%nonlocal = pl%get_asBool('nonlocal') @@ -486,26 +486,26 @@ module function plastic_nonlocal_init() result(myPlasticity) dot%rho_dip_scr => plasticState(ph)%dotState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:) del%rho_dip_scr => plasticState(ph)%deltaState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:) - stt%gamma => plasticState(ph)%state (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents) - dot%gamma => plasticState(ph)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents) - del%gamma => plasticState(ph)%deltaState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents) + stt%gamma => plasticState(ph)%state (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nmembers) + dot%gamma => plasticState(ph)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nmembers) + del%gamma => plasticState(ph)%deltaState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nmembers) plasticState(ph)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl ) = pl%get_asFloat('atol_gamma', defaultVal = 1.0e-2_pReal) if(any(plasticState(ph)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl) < 0.0_pReal)) & extmsg = trim(extmsg)//' atol_gamma' - plasticState(ph)%slipRate => plasticState(ph)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents) + plasticState(ph)%slipRate => plasticState(ph)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nmembers) - stt%rho_forest => plasticState(ph)%state (11*prm%sum_N_sl + 1:12*prm%sum_N_sl,1:Nconstituents) - stt%v => plasticState(ph)%state (12*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nconstituents) - stt%v_edg_pos => plasticState(ph)%state (12*prm%sum_N_sl + 1:13*prm%sum_N_sl,1:Nconstituents) - stt%v_edg_neg => plasticState(ph)%state (13*prm%sum_N_sl + 1:14*prm%sum_N_sl,1:Nconstituents) - stt%v_scr_pos => plasticState(ph)%state (14*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:Nconstituents) - stt%v_scr_neg => plasticState(ph)%state (15*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nconstituents) + stt%rho_forest => plasticState(ph)%state (11*prm%sum_N_sl + 1:12*prm%sum_N_sl,1:Nmembers) + stt%v => plasticState(ph)%state (12*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nmembers) + stt%v_edg_pos => plasticState(ph)%state (12*prm%sum_N_sl + 1:13*prm%sum_N_sl,1:Nmembers) + stt%v_edg_neg => plasticState(ph)%state (13*prm%sum_N_sl + 1:14*prm%sum_N_sl,1:Nmembers) + stt%v_scr_pos => plasticState(ph)%state (14*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:Nmembers) + stt%v_scr_neg => plasticState(ph)%state (15*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nmembers) - allocate(dst%tau_pass(prm%sum_N_sl,Nconstituents),source=0.0_pReal) - allocate(dst%tau_back(prm%sum_N_sl,Nconstituents),source=0.0_pReal) + allocate(dst%tau_pass(prm%sum_N_sl,Nmembers),source=0.0_pReal) + allocate(dst%tau_back(prm%sum_N_sl,Nmembers),source=0.0_pReal) end associate - if (Nconstituents > 0) call stateInit(ini,ph,Nconstituents) + if (Nmembers > 0) call stateInit(ini,ph,Nmembers) plasticState(ph)%state0 = plasticState(ph)%state !-------------------------------------------------------------------------------------------------- @@ -527,7 +527,7 @@ module function plastic_nonlocal_init() result(myPlasticity) if(.not. myPlasticity(ph)) cycle phase => phases%get(ph) - Nconstituents = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseAt2 == ph) l = 0 do t = 1,4 do s = 1,param(ph)%sum_N_sl @@ -1406,7 +1406,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) me, & neighbor_e, & ! element index of my neighbor neighbor_i, & ! integration point index of my neighbor - neighbor_me, & + neighbor_me, & neighbor_phase, & ns, & ! number of active slip systems s1, & ! slip system index (me) @@ -1579,13 +1579,13 @@ end subroutine plastic_nonlocal_results !-------------------------------------------------------------------------------------------------- !> @brief populates the initial dislocation density !-------------------------------------------------------------------------------------------------- -subroutine stateInit(ini,phase,Nconstituents) +subroutine stateInit(ini,phase,Nmembers) type(tInitialParameters) :: & ini integer,intent(in) :: & phase, & - Nconstituents + Nmembers integer :: & i, & e, & @@ -1602,7 +1602,7 @@ subroutine stateInit(ini,phase,Nconstituents) totalVolume, & densityBinning, & minimumIpVolume - real(pReal), dimension(Nconstituents) :: & + real(pReal), dimension(Nmembers) :: & volume @@ -1622,13 +1622,13 @@ subroutine stateInit(ini,phase,Nconstituents) meanDensity = 0.0_pReal do while(meanDensity < ini%random_rho_u) call random_number(rnd) - phasemember = nint(rnd(1)*real(Nconstituents,pReal) + 0.5_pReal) + phasemember = nint(rnd(1)*real(Nmembers,pReal) + 0.5_pReal) s = nint(rnd(2)*real(sum(ini%N_sl),pReal)*4.0_pReal + 0.5_pReal) meanDensity = meanDensity + densityBinning * volume(phasemember) / totalVolume stt%rhoSglMobile(s,phasemember) = densityBinning enddo else ! homogeneous distribution with noise - do e = 1, Nconstituents + do e = 1, Nmembers do f = 1,size(ini%N_sl,1) from = 1 + sum(ini%N_sl(1:f-1)) upto = sum(ini%N_sl(1:f)) @@ -1809,7 +1809,7 @@ pure function getRho0(ph,me) getRho0(:,mob) = max(getRho0(:,mob),0.0_pReal) getRho0(:,dip) = max(getRho0(:,dip),0.0_pReal) - where(abs(getRho0) < max(prm%rho_min/geom(ph)%V_0(me)**(2.0_pReal/3.0_pReal),prm%rho_significant)) & + where (abs(getRho0) < max(prm%rho_min/geom(ph)%V_0(me)**(2.0_pReal/3.0_pReal),prm%rho_significant)) & getRho0 = 0.0_pReal end associate @@ -1823,12 +1823,12 @@ subroutine storeGeometry(ph) integer :: ip, el, ce, co real(pReal), dimension(:), allocatable :: V - + V = reshape(IPvolume,[product(shape(IPvolume))]) do ce = 1, size(material_homogenizationMemberAt2,1) do co = 1, homogenization_maxNconstituents - if(material_phaseAt2(co,ce) == ph) geom(ph)%V_0(material_phaseMemberAt2(co,ce)) = V(ce) + if (material_phaseAt2(co,ce) == ph) geom(ph)%V_0(material_phaseMemberAt2(co,ce)) = V(ce) enddo enddo diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90 index d769431b4..ae5926c0f 100644 --- a/src/phase_mechanical_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanical_plastic_phenopowerlaw.f90 @@ -71,7 +71,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) logical, dimension(:), allocatable :: myPlasticity integer :: & ph, i, & - Nconstituents, & + Nmembers, & sizeState, sizeDotState, & startIndex, endIndex integer, dimension(:), allocatable :: & @@ -223,20 +223,20 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nconstituents = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseAt2 == ph) sizeDotState = size(['xi_sl ','gamma_sl']) * prm%sum_N_sl & + size(['xi_tw ','gamma_tw']) * prm%sum_N_tw sizeState = sizeDotState - call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,0) + call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0) !-------------------------------------------------------------------------------------------------- ! state aliases and initialization startIndex = 1 endIndex = prm%sum_N_sl stt%xi_slip => plasticState(ph)%state (startIndex:endIndex,:) - stt%xi_slip = spread(xi_0_sl, 2, Nconstituents) + stt%xi_slip = spread(xi_0_sl, 2, Nmembers) dot%xi_slip => plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' @@ -244,7 +244,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tw stt%xi_twin => plasticState(ph)%state (startIndex:endIndex,:) - stt%xi_twin = spread(xi_0_tw, 2, Nconstituents) + stt%xi_twin = spread(xi_0_tw, 2, Nmembers) dot%xi_twin => plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90 index 012554aef..21ec93c9c 100644 --- a/src/phase_thermal.f90 +++ b/src/phase_thermal.f90 @@ -15,13 +15,13 @@ submodule(phase) thermal THERMAL_EXTERNALHEAT_ID end enum - type :: tDataContainer + type :: tDataContainer ! ?? not very telling name. Better: "fieldQuantities" ?? real(pReal), dimension(:), allocatable :: T, dot_T end type tDataContainer integer(kind(THERMAL_UNDEFINED_ID)), dimension(:,:), allocatable :: & thermal_source - type(tDataContainer), dimension(:), allocatable :: current + type(tDataContainer), dimension(:), allocatable :: current ! ?? not very telling name. Better: "field" ?? integer :: thermal_source_maxSizeDotState @@ -78,41 +78,36 @@ module subroutine thermal_init(phases) integer :: & ph, so, & - Nconstituents + Nmembers print'(/,a)', ' <<<+- phase:thermal init -+>>>' allocate(current(phases%length)) - allocate(thermalState (phases%length)) + allocate(thermalState(phases%length)) allocate(thermal_Nsources(phases%length),source = 0) do ph = 1, phases%length - - Nconstituents = count(material_phaseAt2 == ph) - - allocate(current(ph)%T(Nconstituents),source=300.0_pReal) - allocate(current(ph)%dot_T(Nconstituents),source=0.0_pReal) + Nmembers = count(material_phaseAt2 == ph) + allocate(current(ph)%T(Nmembers),source=300.0_pReal) + allocate(current(ph)%dot_T(Nmembers),source=0.0_pReal) phase => phases%get(ph) - if(phase%contains('thermal')) then - thermal => phase%get('thermal') - sources => thermal%get('source',defaultVal=emptyList) - - thermal_Nsources(ph) = sources%length - endif + thermal => phase%get('thermal',defaultVal=emptyDict) + sources => thermal%get('source',defaultVal=emptyList) + thermal_Nsources(ph) = sources%length allocate(thermalstate(ph)%p(thermal_Nsources(ph))) enddo allocate(thermal_source(maxval(thermal_Nsources),phases%length), source = THERMAL_UNDEFINED_ID) - if(maxval(thermal_Nsources) /= 0) then + if (maxval(thermal_Nsources) /= 0) then where(dissipation_init (maxval(thermal_Nsources))) thermal_source = THERMAL_DISSIPATION_ID where(externalheat_init(maxval(thermal_Nsources))) thermal_source = THERMAL_EXTERNALHEAT_ID endif thermal_source_maxSizeDotState = 0 - PhaseLoop2:do ph = 1,phases%length + do ph = 1,phases%length do so = 1,thermal_Nsources(ph) thermalState(ph)%p(so)%state = thermalState(ph)%p(so)%state0 @@ -120,7 +115,7 @@ module subroutine thermal_init(phases) thermal_source_maxSizeDotState = max(thermal_source_maxSizeDotState, & maxval(thermalState(ph)%p%sizeDotState)) - enddo PhaseLoop2 + enddo end subroutine thermal_init @@ -145,18 +140,17 @@ module subroutine phase_thermal_getRate(TDot, ph,me) do so = 1, thermal_Nsources(ph) select case(thermal_source(so,ph)) case (THERMAL_DISSIPATION_ID) - call dissipation_getRate(my_Tdot, ph,me) + call dissipation_getRate(my_Tdot, ph,me) case (THERMAL_EXTERNALHEAT_ID) - call externalheat_getRate(my_Tdot, ph,me) + call externalheat_getRate(my_Tdot, ph,me) case default - my_Tdot = 0.0_pReal + my_Tdot = 0.0_pReal end select Tdot = Tdot + my_Tdot enddo - end subroutine phase_thermal_getRate @@ -185,7 +179,7 @@ function phase_thermal_collectDotState(ph,me) result(broken) end function phase_thermal_collectDotState -module function thermal_stress(Delta_t,ph,me) result(converged_) +module function thermal_stress(Delta_t,ph,me) result(converged_) ! ?? why is this called "stress" when it seems closer to "updateState" ?? real(pReal), intent(in) :: Delta_t integer, intent(in) :: ph, me @@ -212,7 +206,7 @@ function integrateThermalState(Delta_t, ph,me) result(broken) sizeDotState broken = phase_thermal_collectDotState(ph,me) - if(broken) return + if (broken) return do so = 1, thermal_Nsources(ph) sizeDotState = thermalState(ph)%p(so)%sizeDotState @@ -301,14 +295,12 @@ function thermal_active(source_label,src_length) result(active_source) allocate(active_source(src_length,phases%length), source = .false. ) do p = 1, phases%length phase => phases%get(p) - if (phase%contains('thermal')) then - thermal => phase%get('thermal',defaultVal=emptyList) - sources => thermal%get('source',defaultVal=emptyList) - do s = 1, sources%length - src => sources%get(s) - if(src%get_asString('type') == source_label) active_source(s,p) = .true. - enddo - endif + thermal => phase%get('thermal',defaultVal=emptyDict) + sources => thermal%get('source',defaultVal=emptyList) + do s = 1, sources%length + src => sources%get(s) + active_source(s,p) = src%get_asString('type') == source_label + enddo enddo diff --git a/src/phase_thermal_dissipation.f90 b/src/phase_thermal_dissipation.f90 index 9c75763dc..d3e7094a1 100644 --- a/src/phase_thermal_dissipation.f90 +++ b/src/phase_thermal_dissipation.f90 @@ -31,7 +31,7 @@ module function dissipation_init(source_length) result(mySources) phase, & sources, thermal, & src - integer :: so,Nconstituents,ph + integer :: so,Nmembers,ph mySources = thermal_active('dissipation',source_length) @@ -54,8 +54,8 @@ module function dissipation_init(source_length) result(mySources) src => sources%get(so) prm%kappa = src%get_asFloat('kappa') - Nconstituents = count(material_phaseAt2 == ph) - call phase_allocateState(thermalState(ph)%p(so),Nconstituents,0,0,0) + Nmembers = count(material_phaseAt2 == ph) + call phase_allocateState(thermalState(ph)%p(so),Nmembers,0,0,0) end associate endif diff --git a/src/phase_thermal_externalheat.f90 b/src/phase_thermal_externalheat.f90 index c7f894215..257b4e282 100644 --- a/src/phase_thermal_externalheat.f90 +++ b/src/phase_thermal_externalheat.f90 @@ -38,7 +38,7 @@ module function externalheat_init(source_length) result(mySources) phase, & sources, thermal, & src - integer :: so,Nconstituents,ph + integer :: so,Nmembers,ph mySources = thermal_active('externalheat',source_length) @@ -67,8 +67,8 @@ module function externalheat_init(source_length) result(mySources) prm%f_T = src%get_asFloats('f_T',requiredSize = size(prm%t_n)) - Nconstituents = count(material_phaseAt2 == ph) - call phase_allocateState(thermalState(ph)%p(so),Nconstituents,1,1,0) + Nmembers = count(material_phaseAt2 == ph) + call phase_allocateState(thermalState(ph)%p(so),Nmembers,1,1,0) end associate endif enddo From 5cd68dbddb85191664e3aa78a7b9472ce42bdcb6 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Thu, 4 Mar 2021 15:35:17 -0500 Subject: [PATCH 13/13] fixed stray homogenization_Nmembers --> _Nconstituents --- src/phase_damage.f90 | 2 +- src/phase_mechanical.f90 | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 6ad25f1e9..c85075288 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -199,7 +199,7 @@ module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, phiDot = 0.0_pReal dPhiDot_dPhi = 0.0_pReal - do co = 1, homogenization_Nmembers(material_homogenizationAt2(ce)) + do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) ph = material_phaseAt2(co,ce) me = material_phasememberAt2(co,ce) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index d065fc74a..e642b22ef 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -280,7 +280,7 @@ module subroutine mechanical_init(materials,phases) do el = 1, size(material_phaseMemberAt,3); do ip = 1, size(material_phaseMemberAt,2) - do co = 1, homogenization_Nmembers(material_homogenizationAt(el)) + do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) material => materials%get(discretization_materialAt(el)) constituents => material%get('constituents') constituent => constituents%get(co) @@ -1238,7 +1238,7 @@ module subroutine mechanical_restore(ce,includeL) co, ph, me - do co = 1,homogenization_Nmembers(material_homogenizationAt2(ce)) + do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce)) ph = material_phaseAt2(co,ce) me = material_phaseMemberAt2(co,ce) if (includeL) then