From c6863a612407b2e746ee54bb5638e27706b8b6fa Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 15 May 2019 23:32:23 +0200 Subject: [PATCH 01/23] also consider homogenization/materialpoint results --- python/damask/dadf5.py | 23 +++++++++++++++++++++-- 1 file changed, 21 insertions(+), 2 deletions(-) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 841029af8..eed5907e5 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -48,16 +48,24 @@ class DADF5(): for o in f['inc{:05}/constituent/{}'.format(self.increments[0]['inc'],c)].keys(): self.c_output_types.append(o) self.c_output_types = list(set(self.c_output_types)) # make unique + + self.m_output_types = [] + for m in self.materialpoints: + for o in f['inc{:05}/materialpoint/{}'.format(self.increments[0]['inc'],m)].keys(): + self.m_output_types.append(o) + self.m_output_types = list(set(self.m_output_types)) # make unique self.active= {'increments': self.increments, 'constituents': self.constituents, 'materialpoints': self.materialpoints, 'constituent': self.Nconstituents, - 'c_output_types': self.c_output_types} + 'c_output_types': self.c_output_types, + 'm_output_types': self.m_output_types} self.filename = filename self.mode = mode - + + def list_data(self): """Shows information on all datasets in the file""" with h5py.File(self.filename,'r') as f: @@ -73,6 +81,17 @@ class DADF5(): print(' {} ({})'.format(x,f[group_output_types+'/'+x].attrs['Description'].decode())) except: pass + for m in self.active['materialpoints']: + print('\n'+m) + group_materialpoint = group_inc+'/materialpoint/'+m + for t in self.active['m_output_types']: + print(' {}'.format(t)) + group_output_types = group_materialpoint+'/'+t + try: + for x in f[group_output_types].keys(): + print(' {} ({})'.format(x,f[group_output_types+'/'+x].attrs['Description'].decode())) + except: + pass def get_dataset_location(self,label): From 3c0c0a2cd10d38c472b92d84cb687f2268653b02 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 16 May 2019 00:27:06 +0200 Subject: [PATCH 02/23] more flexible in selecting data --- processing/post/DADF5_postResults.py | 28 ++++++++++++++++++++++++---- processing/post/DADF5_vtk_cells.py | 9 +++++++-- python/damask/dadf5.py | 26 +++++++++++++++++++++----- 3 files changed, 52 insertions(+), 11 deletions(-) diff --git a/processing/post/DADF5_postResults.py b/processing/post/DADF5_postResults.py index caf09d536..5ad49604b 100755 --- a/processing/post/DADF5_postResults.py +++ b/processing/post/DADF5_postResults.py @@ -22,10 +22,15 @@ parser.add_argument('filenames', nargs='+', help='DADF5 files') parser.add_argument('-d','--dir', dest='dir',default='postProc',metavar='string', help='name of subdirectory to hold output') +parser.add_argument('--mat', nargs='+', + help='labels for materialpoint/homogenization',dest='mat') +parser.add_argument('--con', nargs='+', + help='labels for constituent/crystallite/constitutive',dest='con') options = parser.parse_args() -options.labels = ['Fe','Fp','xi_sl'] +if options.mat is None: options.mat=[] +if options.con is None: options.con=[] # --- loop over input files ------------------------------------------------------------------------ @@ -57,7 +62,7 @@ for filename in options.filenames: header+=' 1_pos 2_pos 3_pos' results.active['increments'] = [inc] - for label in options.labels: + for label in options.con: for o in results.c_output_types: results.active['c_output_types'] = [o] for c in results.constituents: @@ -67,12 +72,27 @@ for filename in options.filenames: continue label = x[0].split('/')[-1] array = results.read_dataset(x,0) - d = np.product(np.shape(array)[1:]) + d = int(np.product(np.shape(array)[1:])) + array = np.reshape(array,[np.product(results.grid),d]) + data = np.concatenate((data,array),1) + + header+= ''.join([' {}_{}'.format(j+1,label) for j in range(d)]) + + for label in options.mat: + for o in results.m_output_types: + results.active['m_output_types'] = [o] + for m in results.materialpoints: + results.active['materialpoints'] = [m] + x = results.get_dataset_location(label) + if len(x) == 0: + continue + label = x[0].split('/')[-1] + array = results.read_dataset(x,0) + d = int(np.product(np.shape(array)[1:])) array = np.reshape(array,[np.product(results.grid),d]) data = np.concatenate((data,array),1) header+= ''.join([' {}_{}'.format(j+1,label) for j in range(d)]) - dirname = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir)) try: diff --git a/processing/post/DADF5_vtk_cells.py b/processing/post/DADF5_vtk_cells.py index dc1177488..0e79209e2 100755 --- a/processing/post/DADF5_vtk_cells.py +++ b/processing/post/DADF5_vtk_cells.py @@ -23,10 +23,15 @@ parser.add_argument('filenames', nargs='+', help='DADF5 files') parser.add_argument('-d','--dir', dest='dir',default='postProc',metavar='string', help='name of subdirectory to hold output') +parser.add_argument('--mat', nargs='+', + help='labels for materialpoint/homogenization',dest='mat') +parser.add_argument('--con', nargs='+', + help='labels for constituent/crystallite/constitutive',dest='con') options = parser.parse_args() -options.labels = ['Fe','Fp','xi_sl'] +if options.mat is None: options.mat=[] +if options.con is None: options.con=[] # --- loop over input files ------------------------------------------------------------------------ @@ -54,7 +59,7 @@ for filename in options.filenames: print('Output step {}/{}'.format(i+1,len(results.increments))) vtk_data = [] results.active['increments'] = [inc] - for label in options.labels: + for label in options.con: for o in results.c_output_types: results.active['c_output_types'] = [o] if o != 'generic': diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index eed5907e5..56082a1f4 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -82,7 +82,6 @@ class DADF5(): except: pass for m in self.active['materialpoints']: - print('\n'+m) group_materialpoint = group_inc+'/materialpoint/'+m for t in self.active['m_output_types']: print(' {}'.format(t)) @@ -108,6 +107,14 @@ class DADF5(): path.append(group_constituent+'/'+t+'/'+label) except: pass + for m in self.active['materialpoints']: + group_materialpoint = group_inc+'/materialpoint/'+m + for t in self.active['m_output_types']: + try: + f[group_materialpoint+'/'+t+'/'+label] + path.append(group_materialpoint+'/'+t+'/'+label) + except: + pass return path @@ -122,10 +129,19 @@ class DADF5(): dataset = np.full(shape,np.nan) for pa in path: label = pa.split('/')[2] - p = np.where(f['mapping/cellResults/constituent'][:,c]['Name'] == str.encode(label))[0] - u = (f['mapping/cellResults/constituent'][p,c]['Position']) - dataset[p,:] = f[pa][u,:] - + try: + p = np.where(f['mapping/cellResults/constituent'][:,c]['Name'] == str.encode(label))[0] + u = (f['mapping/cellResults/constituent'][p,c]['Position']) + dataset[p,:] = f[pa][u,:] # does not work for scalar datasets + except: + pass + try: + p = np.where(f['mapping/cellResults/materialpoint']['Name'] == str.encode(label))[0] + u = (f['mapping/cellResults/materialpoint'][p.tolist()]['Position']) + dataset[p,:] = f[pa][u,:] # does not work for scalar datasets + except: + pass + return dataset From 4599d1c34e219d9c1ab3d49971d4d37ab65049e3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 16 May 2019 00:42:22 +0200 Subject: [PATCH 03/23] does not match node does not make sense, is weirdly numbered --- processing/post/DADF5_postResults.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/processing/post/DADF5_postResults.py b/processing/post/DADF5_postResults.py index 5ad49604b..b918388d9 100755 --- a/processing/post/DADF5_postResults.py +++ b/processing/post/DADF5_postResults.py @@ -53,9 +53,6 @@ for filename in options.filenames: data = np.array([inc['inc'] for j in range(np.product(results.grid))]).reshape([np.product(results.grid),1]) header+= 'inc' - - data = np.concatenate((data,np.array([j+1 for j in range(np.product(results.grid))]).reshape([np.product(results.grid),1])),1) - header+=' node' coords = coords.reshape([np.product(results.grid),3]) data = np.concatenate((data,coords),1) From 9f7fa5393a8af1d784433262dfbbe8bebded5ab0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 16 May 2019 09:31:13 +0200 Subject: [PATCH 04/23] fix for scalar datasets --- processing/post/DADF5_postResults.py | 12 +++++++++--- python/damask/dadf5.py | 11 +++++++++-- 2 files changed, 18 insertions(+), 5 deletions(-) diff --git a/processing/post/DADF5_postResults.py b/processing/post/DADF5_postResults.py index b918388d9..1af8c2aa9 100755 --- a/processing/post/DADF5_postResults.py +++ b/processing/post/DADF5_postResults.py @@ -73,8 +73,11 @@ for filename in options.filenames: array = np.reshape(array,[np.product(results.grid),d]) data = np.concatenate((data,array),1) - header+= ''.join([' {}_{}'.format(j+1,label) for j in range(d)]) - + if d>1: + header+= ''.join([' {}_{}'.format(j+1,label) for j in range(d)]) + else: + header+=' '+label + for label in options.mat: for o in results.m_output_types: results.active['m_output_types'] = [o] @@ -89,7 +92,10 @@ for filename in options.filenames: array = np.reshape(array,[np.product(results.grid),d]) data = np.concatenate((data,array),1) - header+= ''.join([' {}_{}'.format(j+1,label) for j in range(d)]) + if d>1: + header+= ''.join([' {}_{}'.format(j+1,label) for j in range(d)]) + else: + header+=' '+label dirname = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir)) try: diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 56082a1f4..0abde662b 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -126,19 +126,26 @@ class DADF5(): """ with h5py.File(self.filename,'r') as f: shape = (self.Nmaterialpoints,) + np.shape(f[path[0]])[1:] + if len(shape) == 1: shape = shape +(1,) dataset = np.full(shape,np.nan) for pa in path: label = pa.split('/')[2] try: p = np.where(f['mapping/cellResults/constituent'][:,c]['Name'] == str.encode(label))[0] u = (f['mapping/cellResults/constituent'][p,c]['Position']) - dataset[p,:] = f[pa][u,:] # does not work for scalar datasets + a = np.array(f[pa]) + if len(a.shape) == 1: + a=a.reshape([a.shape[0],1]) + dataset except: pass try: p = np.where(f['mapping/cellResults/materialpoint']['Name'] == str.encode(label))[0] u = (f['mapping/cellResults/materialpoint'][p.tolist()]['Position']) - dataset[p,:] = f[pa][u,:] # does not work for scalar datasets + a = np.array(f[pa]) + if len(a.shape) == 1: + a=a.reshape([a.shape[0],1]) + dataset[p,:] = a[u,:] except: pass From f2268d055ff9889e11aeafbc8d0a6c34e247a75f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 16 May 2019 09:42:37 +0200 Subject: [PATCH 05/23] first test relying on DADF5 --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index d96bfb329..df0b7cc97 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit d96bfb32920c96a8a43958f76a209d34c6bd841a +Subproject commit df0b7cc97963cd712f6e33397937b187635c99f4 From 84a82d0878f442067ac6c3fba30276c0c7523f21 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 16 May 2019 09:46:50 +0200 Subject: [PATCH 06/23] rely on working HDF5 output for RGC homog --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index df0b7cc97..8cec21296 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit df0b7cc97963cd712f6e33397937b187635c99f4 +Subproject commit 8cec2129617f4a206193671e0e1070965c7c2e53 From 0958c4bb88bc441e598a460eb785344ca9edf9dd Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 16 May 2019 09:48:30 +0200 Subject: [PATCH 07/23] phase out postResults --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 8cec21296..639c6f4a5 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 8cec2129617f4a206193671e0e1070965c7c2e53 +Subproject commit 639c6f4a5eafc893c83c740c57f417eaaabc45ae From 39a75c201545327d62a37df8f3b8f89b21cf3a32 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 16 May 2019 10:23:23 +0200 Subject: [PATCH 08/23] phasing out postResults starting with RGC because it is rarely used and removing it here allows to go ahead with the submodule structure for homogenization --- src/homogenization.f90 | 44 +---------------- src/homogenization_mech_RGC.f90 | 88 +++------------------------------ 2 files changed, 10 insertions(+), 122 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 77a2e004f..1ff036c3a 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -120,7 +120,6 @@ subroutine homogenization_init character(len=32) :: outputName !< name of output, intermediate fix until HDF5 output is ready logical :: valid - if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call mech_none_init if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call mech_isostrain_init if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call homogenization_RGC_init @@ -139,35 +138,9 @@ subroutine homogenization_init call IO_write_jobFile(FILEUNIT,'outputHomogenization') do p = 1,size(config_homogenization) if (any(material_homogenizationAt == p)) then - i = homogenization_typeInstance(p) ! which instance of this homogenization type - valid = .true. ! assume valid - select case(homogenization_type(p)) ! split per homogenization type - case (HOMOGENIZATION_NONE_ID) - outputName = HOMOGENIZATION_NONE_label - thisOutput => null() - thisSize => null() - case (HOMOGENIZATION_ISOSTRAIN_ID) - outputName = HOMOGENIZATION_ISOSTRAIN_label - thisOutput => null() - thisSize => null() - case (HOMOGENIZATION_RGC_ID) - outputName = HOMOGENIZATION_RGC_label - thisOutput => homogenization_RGC_output - thisSize => homogenization_RGC_sizePostResult - case default - valid = .false. - end select write(FILEUNIT,'(/,a,/)') '['//trim(homogenization_name(p))//']' - if (valid) then - write(FILEUNIT,'(a)') '(type)'//char(9)//trim(outputName) - write(FILEUNIT,'(a,i4)') '(ngrains)'//char(9),homogenization_Ngrains(p) - if (homogenization_type(p) /= HOMOGENIZATION_NONE_ID .and. & - homogenization_type(p) /= HOMOGENIZATION_ISOSTRAIN_ID) then - do e = 1,size(thisOutput(:,i)) - write(FILEUNIT,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i) - enddo - endif - endif + write(FILEUNIT,'(a,i4)') '(ngrains)'//char(9),homogenization_Ngrains(p) + i = thermal_typeInstance(p) ! which instance of this thermal type valid = .true. ! assume valid select case(thermal_type(p)) ! split per thermal type @@ -837,8 +810,6 @@ end subroutine averageStressAndItsTangent function postResults(ip,el) use mesh, only: & mesh_element - use homogenization_mech_RGC, only: & - homogenization_RGC_postResults use thermal_adiabatic, only: & thermal_adiabatic_postResults use thermal_conduction, only: & @@ -861,17 +832,6 @@ function postResults(ip,el) postResults = 0.0_pReal - startPos = 1 - endPos = homogState(material_homogenizationAt(el))%sizePostResults - chosenHomogenization: select case (homogenization_type(mesh_element(3,el))) - - case (HOMOGENIZATION_RGC_ID) chosenHomogenization - instance = homogenization_typeInstance(material_homogenizationAt(el)) - of = mappingHomogenization(1,ip,el) - postResults(startPos:endPos) = homogenization_RGC_postResults(instance,of) - - end select chosenHomogenization - startPos = endPos + 1 endPos = endPos + thermalState(material_homogenizationAt(el))%sizePostResults chosenThermal: select case (thermal_type(mesh_element(3,el))) diff --git a/src/homogenization_mech_RGC.f90 b/src/homogenization_mech_RGC.f90 index d3feac1eb..d7b1b31bf 100644 --- a/src/homogenization_mech_RGC.f90 +++ b/src/homogenization_mech_RGC.f90 @@ -12,10 +12,6 @@ module homogenization_mech_RGC implicit none private - integer, dimension(:,:), allocatable,target, public :: & - homogenization_RGC_sizePostResult - character(len=64), dimension(:,:), allocatable,target, public :: & - homogenization_RGC_output ! name of each post result output enum, bind(c) enumerator :: & @@ -28,7 +24,7 @@ module homogenization_mech_RGC magnitudemismatch_ID end enum - type, private :: tParameters + type :: tParameters integer, dimension(:), allocatable :: & Nconstituents real(pReal) :: & @@ -43,7 +39,7 @@ module homogenization_mech_RGC outputID end type tParameters - type, private :: tRGCstate + type :: tRGCstate real(pReal), pointer, dimension(:) :: & work, & penaltyEnergy @@ -51,7 +47,7 @@ module homogenization_mech_RGC relaxationVector end type tRGCstate - type, private :: tRGCdependentState + type :: tRGCdependentState real(pReal), allocatable, dimension(:) :: & volumeDiscrepancy, & relaxationRate_avg, & @@ -62,12 +58,12 @@ module homogenization_mech_RGC orientation end type tRGCdependentState - type(tparameters), dimension(:), allocatable, private :: & + type(tparameters), dimension(:), allocatable :: & param - type(tRGCstate), dimension(:), allocatable, private :: & + type(tRGCstate), dimension(:), allocatable :: & state, & state0 - type(tRGCdependentState), dimension(:), allocatable, private :: & + type(tRGCdependentState), dimension(:), allocatable :: & dependentState public :: & @@ -75,16 +71,7 @@ module homogenization_mech_RGC homogenization_RGC_partitionDeformation, & homogenization_RGC_averageStressAndItsTangent, & homogenization_RGC_updateState, & - homogenization_RGC_postResults, & mech_RGC_results ! name suited for planned submodule situation - private :: & - relaxationVector, & - interfaceNormal, & - getInterface, & - grain1to3, & - grain3to1, & - interface4to1, & - interface1to4 contains @@ -111,7 +98,7 @@ subroutine homogenization_RGC_init() integer :: & Ninstance, & h, i, & - NofMyHomog, outputSize, & + NofMyHomog, & sizeState, nIntFaceTot character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] @@ -139,9 +126,6 @@ subroutine homogenization_RGC_init() allocate(state0(Ninstance)) allocate(dependentState(Ninstance)) - allocate(homogenization_RGC_sizePostResult(maxval(homogenization_Noutput),Ninstance),source=0) - allocate(homogenization_RGC_output(maxval(homogenization_Noutput),Ninstance)) - homogenization_RGC_output='' do h = 1, size(homogenization_type) if (homogenization_type(h) /= HOMOGENIZATION_RGC_ID) cycle @@ -176,28 +160,20 @@ subroutine homogenization_RGC_init() case('constitutivework') outputID = constitutivework_ID - outputSize = 1 case('penaltyenergy') outputID = penaltyenergy_ID - outputSize = 1 case('volumediscrepancy') outputID = volumediscrepancy_ID - outputSize = 1 case('averagerelaxrate') outputID = averagerelaxrate_ID - outputSize = 1 case('maximumrelaxrate') outputID = maximumrelaxrate_ID - outputSize = 1 case('magnitudemismatch') outputID = magnitudemismatch_ID - outputSize = 3 end select if (outputID /= undefined_ID) then - homogenization_RGC_output(i,homogenization_typeInstance(h)) = outputs(i) - homogenization_RGC_sizePostResult(i,homogenization_typeInstance(h)) = outputSize prm%outputID = [prm%outputID , outputID] endif @@ -211,7 +187,7 @@ subroutine homogenization_RGC_init() + size(['avg constitutive work ','average penalty energy']) homogState(h)%sizeState = sizeState - homogState(h)%sizePostResults = sum(homogenization_RGC_sizePostResult(:,homogenization_typeInstance(h))) + homogState(h)%sizePostResults = 0 allocate(homogState(h)%state0 (sizeState,NofMyHomog), source=0.0_pReal) allocate(homogState(h)%subState0(sizeState,NofMyHomog), source=0.0_pReal) allocate(homogState(h)%state (sizeState,NofMyHomog), source=0.0_pReal) @@ -1033,54 +1009,6 @@ subroutine homogenization_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF, end subroutine homogenization_RGC_averageStressAndItsTangent -!-------------------------------------------------------------------------------------------------- -!> @brief return array of homogenization results for post file inclusion -!-------------------------------------------------------------------------------------------------- -pure function homogenization_RGC_postResults(instance,of) result(postResults) - - integer, intent(in) :: & - instance, & - of - - integer :: & - o,c - real(pReal), dimension(sum(homogenization_RGC_sizePostResult(:,instance))) :: & - postResults - - associate(stt => state(instance), dst => dependentState(instance), prm => param(instance)) - - c = 0 - - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) - - case (constitutivework_ID) - postResults(c+1) = stt%work(of) - c = c + 1 - case (magnitudemismatch_ID) - postResults(c+1:c+3) = dst%mismatch(1:3,of) - c = c + 3 - case (penaltyenergy_ID) - postResults(c+1) = stt%penaltyEnergy(of) - c = c + 1 - case (volumediscrepancy_ID) - postResults(c+1) = dst%volumeDiscrepancy(of) - c = c + 1 - case (averagerelaxrate_ID) - postResults(c+1) = dst%relaxationrate_avg(of) - c = c + 1 - case (maximumrelaxrate_ID) - postResults(c+1) = dst%relaxationrate_max(of) - c = c + 1 - end select - - enddo outputsLoop - - end associate - -end function homogenization_RGC_postResults - - !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file ! ToDo: check wheter units are correct From 339b86f784bdb31c1359b486203e12d1e1d0717a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 16 May 2019 11:44:03 +0200 Subject: [PATCH 09/23] bugfix + more verbose reporting --- processing/post/DADF5_postResults.py | 2 +- processing/post/DADF5_vtk_cells.py | 2 ++ python/damask/dadf5.py | 23 +++++++++++++---------- 3 files changed, 16 insertions(+), 11 deletions(-) diff --git a/processing/post/DADF5_postResults.py b/processing/post/DADF5_postResults.py index 1af8c2aa9..fa47805bb 100755 --- a/processing/post/DADF5_postResults.py +++ b/processing/post/DADF5_postResults.py @@ -78,7 +78,7 @@ for filename in options.filenames: else: header+=' '+label - for label in options.mat: + for label in options.mat: for o in results.m_output_types: results.active['m_output_types'] = [o] for m in results.materialpoints: diff --git a/processing/post/DADF5_vtk_cells.py b/processing/post/DADF5_vtk_cells.py index 0e79209e2..ef27a973c 100755 --- a/processing/post/DADF5_vtk_cells.py +++ b/processing/post/DADF5_vtk_cells.py @@ -59,7 +59,9 @@ for filename in options.filenames: print('Output step {}/{}'.format(i+1,len(results.increments))) vtk_data = [] results.active['increments'] = [inc] + for label in options.con: + for o in results.c_output_types: results.active['c_output_types'] = [o] if o != 'generic': diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 0abde662b..4aa31ec56 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -99,22 +99,25 @@ class DADF5(): with h5py.File(self.filename,'r') as f: for i in self.active['increments']: group_inc = 'inc{:05}'.format(i['inc']) + for c in self.active['constituents']: group_constituent = group_inc+'/constituent/'+c for t in self.active['c_output_types']: try: f[group_constituent+'/'+t+'/'+label] path.append(group_constituent+'/'+t+'/'+label) - except: - pass - for m in self.active['materialpoints']: + except Exception as e: + print('unable to locate constituents dataset: '+ str(e)) + + for m in []: #self.active['materialpoints']: group_materialpoint = group_inc+'/materialpoint/'+m for t in self.active['m_output_types']: try: f[group_materialpoint+'/'+t+'/'+label] path.append(group_materialpoint+'/'+t+'/'+label) - except: - pass + except Exception as e: + print('unable to locate materialpoints dataset: '+ str(e)) + return path @@ -136,9 +139,9 @@ class DADF5(): a = np.array(f[pa]) if len(a.shape) == 1: a=a.reshape([a.shape[0],1]) - dataset - except: - pass + dataset[p,:] = a[u,:] + except Exception as e: + print('unable to read constituent: '+ str(e)) try: p = np.where(f['mapping/cellResults/materialpoint']['Name'] == str.encode(label))[0] u = (f['mapping/cellResults/materialpoint'][p.tolist()]['Position']) @@ -146,8 +149,8 @@ class DADF5(): if len(a.shape) == 1: a=a.reshape([a.shape[0],1]) dataset[p,:] = a[u,:] - except: - pass + except Exception as e: + print('unable to read materialpoint: '+ str(e)) return dataset From 6df563624dc858ed5204c7181006fcddeb0b956a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 16 May 2019 21:57:20 +0200 Subject: [PATCH 10/23] type needed for postResults --- src/homogenization.f90 | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 1ff036c3a..426b37cb8 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -138,9 +138,23 @@ subroutine homogenization_init call IO_write_jobFile(FILEUNIT,'outputHomogenization') do p = 1,size(config_homogenization) if (any(material_homogenizationAt == p)) then + i = homogenization_typeInstance(p) ! which instance of this homogenization type + valid = .true. ! assume valid + select case(homogenization_type(p)) ! split per homogenization type + case (HOMOGENIZATION_NONE_ID) + outputName = HOMOGENIZATION_NONE_label + case (HOMOGENIZATION_ISOSTRAIN_ID) + outputName = HOMOGENIZATION_ISOSTRAIN_label + case (HOMOGENIZATION_RGC_ID) + outputName = HOMOGENIZATION_RGC_label + case default + valid = .false. + end select write(FILEUNIT,'(/,a,/)') '['//trim(homogenization_name(p))//']' - write(FILEUNIT,'(a,i4)') '(ngrains)'//char(9),homogenization_Ngrains(p) - + if (valid) then + write(FILEUNIT,'(a)') '(type)'//char(9)//trim(outputName) + write(FILEUNIT,'(a,i4)') '(ngrains)'//char(9),homogenization_Ngrains(p) + endif i = thermal_typeInstance(p) ! which instance of this thermal type valid = .true. ! assume valid select case(thermal_type(p)) ! split per thermal type From 118c74a268ef9417903576fd7f886d2b4ab569c1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 16 May 2019 22:01:12 +0200 Subject: [PATCH 11/23] enable materialpoint output again accidently lost during a former commit --- python/damask/dadf5.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 4aa31ec56..f5fc81b16 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -109,7 +109,7 @@ class DADF5(): except Exception as e: print('unable to locate constituents dataset: '+ str(e)) - for m in []: #self.active['materialpoints']: + for m in self.active['materialpoints']: group_materialpoint = group_inc+'/materialpoint/'+m for t in self.active['m_output_types']: try: From dce4775c17abed2b2ca09dd8b1acd53a8e477ccb Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 17 May 2019 06:36:30 +0200 Subject: [PATCH 12/23] removal of RGC out led to undefined variable --- src/homogenization.f90 | 30 +++++++++--------------------- 1 file changed, 9 insertions(+), 21 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 426b37cb8..c83322f26 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -134,27 +134,14 @@ subroutine homogenization_init !-------------------------------------------------------------------------------------------------- ! write description file for homogenization output - mainProcess2: if (worldrank == 0) then + mainProcess: if (worldrank == 0) then call IO_write_jobFile(FILEUNIT,'outputHomogenization') do p = 1,size(config_homogenization) if (any(material_homogenizationAt == p)) then - i = homogenization_typeInstance(p) ! which instance of this homogenization type - valid = .true. ! assume valid - select case(homogenization_type(p)) ! split per homogenization type - case (HOMOGENIZATION_NONE_ID) - outputName = HOMOGENIZATION_NONE_label - case (HOMOGENIZATION_ISOSTRAIN_ID) - outputName = HOMOGENIZATION_ISOSTRAIN_label - case (HOMOGENIZATION_RGC_ID) - outputName = HOMOGENIZATION_RGC_label - case default - valid = .false. - end select write(FILEUNIT,'(/,a,/)') '['//trim(homogenization_name(p))//']' - if (valid) then - write(FILEUNIT,'(a)') '(type)'//char(9)//trim(outputName) - write(FILEUNIT,'(a,i4)') '(ngrains)'//char(9),homogenization_Ngrains(p) - endif + write(FILEUNIT,'(a)') '(type) n/a' + write(FILEUNIT,'(a,i4)') '(ngrains)'//char(9),homogenization_Ngrains(p) + i = thermal_typeInstance(p) ! which instance of this thermal type valid = .true. ! assume valid select case(thermal_type(p)) ! split per thermal type @@ -184,6 +171,7 @@ subroutine homogenization_init enddo endif endif + i = damage_typeInstance(p) ! which instance of this damage type valid = .true. ! assume valid select case(damage_type(p)) ! split per damage type @@ -216,7 +204,7 @@ subroutine homogenization_init endif enddo close(FILEUNIT) - endif mainProcess2 + endif mainProcess call config_deallocate('material.config/homogenization') @@ -842,12 +830,12 @@ function postResults(ip,el) postResults integer :: & startPos, endPos ,& - of, instance, homog + homog postResults = 0.0_pReal - startPos = endPos + 1 - endPos = endPos + thermalState(material_homogenizationAt(el))%sizePostResults + startPos = 1 + endPos = thermalState(material_homogenizationAt(el))%sizePostResults chosenThermal: select case (thermal_type(mesh_element(3,el))) case (THERMAL_adiabatic_ID) chosenThermal From ed8af98d695a670f3a8df0e66975d804d0851723 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 18 May 2019 06:54:45 +0200 Subject: [PATCH 13/23] don't clutter with use statements --- src/homogenization.f90 | 111 +++----------------------- src/homogenization_mech_RGC.f90 | 75 +++-------------- src/homogenization_mech_isostrain.f90 | 13 +-- src/homogenization_mech_none.f90 | 8 +- 4 files changed, 29 insertions(+), 178 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index c83322f26..d7b5b3fdc 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -6,7 +6,20 @@ !-------------------------------------------------------------------------------------------------- module homogenization use prec + use IO + use config + use debug + use math use material + use numerics + use constitutive + use crystallite + use mesh + use FEsolving +#if defined(PETSc) || defined(DAMASK_HDF5) + use results + use HDF5_utilities +#endif !-------------------------------------------------------------------------------------------------- ! General variables for the homogenization at a material point @@ -81,26 +94,6 @@ contains !> @brief module initialization !-------------------------------------------------------------------------------------------------- subroutine homogenization_init - use math, only: & - math_I3 - use debug, only: & - debug_level, & - debug_homogenization, & - debug_levelBasic, & - debug_e, & - debug_g - use mesh, only: & - theMesh, & - mesh_element - use constitutive, only: & - constitutive_plasticity_maxSizePostResults, & - constitutive_source_maxSizePostResults - use crystallite, only: & - crystallite_maxSizePostResults - use config, only: & - config_deallocate, & - config_homogenization, & - homogenization_name use homogenization_mech_RGC use thermal_isothermal use thermal_adiabatic @@ -108,9 +101,6 @@ subroutine homogenization_init use damage_none use damage_local use damage_nonlocal - use IO - use numerics, only: & - worldrank integer, parameter :: FILEUNIT = 200 integer :: e,i,p @@ -278,51 +268,6 @@ end subroutine homogenization_init !> @brief parallelized calculation of stress and corresponding tangent at material points !-------------------------------------------------------------------------------------------------- subroutine materialpoint_stressAndItsTangent(updateJaco,dt) - use numerics, only: & - subStepMinHomog, & - subStepSizeHomog, & - stepIncreaseHomog, & - nMPstate - use FEsolving, only: & - FEsolving_execElem, & - FEsolving_execIP, & - terminallyIll - use mesh, only: & - mesh_element - use crystallite, only: & - crystallite_F0, & - crystallite_Fp0, & - crystallite_Fp, & - crystallite_Fi0, & - crystallite_Fi, & - crystallite_Lp0, & - crystallite_Lp, & - crystallite_Li0, & - crystallite_Li, & - crystallite_S0, & - crystallite_S, & - crystallite_partionedF0, & - crystallite_partionedF, & - crystallite_partionedFp0, & - crystallite_partionedLp0, & - crystallite_partionedFi0, & - crystallite_partionedLi0, & - crystallite_partionedS0, & - crystallite_dt, & - crystallite_requested, & - crystallite_stress, & - crystallite_stressTangent, & - crystallite_orientations -#ifdef DEBUG - use debug, only: & - debug_level, & - debug_homogenization, & - debug_levelBasic, & - debug_levelExtensive, & - debug_levelSelective, & - debug_e, & - debug_i -#endif real(pReal), intent(in) :: dt !< time increment logical, intent(in) :: updateJaco !< initiating Jacobian update @@ -616,14 +561,6 @@ end subroutine materialpoint_stressAndItsTangent !> @brief parallelized calculation of result array at material points !-------------------------------------------------------------------------------------------------- subroutine materialpoint_postResults - use FEsolving, only: & - FEsolving_execElem, & - FEsolving_execIP - use mesh, only: & - mesh_element - use crystallite, only: & - crystallite_sizePostResults, & - crystallite_postResults integer :: & thePos, & @@ -673,10 +610,6 @@ end subroutine materialpoint_postResults !> @brief partition material point def grad onto constituents !-------------------------------------------------------------------------------------------------- subroutine partitionDeformation(ip,el) - use mesh, only: & - mesh_element - use crystallite, only: & - crystallite_partionedF use homogenization_mech_RGC, only: & homogenization_RGC_partitionDeformation @@ -710,13 +643,6 @@ end subroutine partitionDeformation !> "happy" with result !-------------------------------------------------------------------------------------------------- function updateState(ip,el) - use mesh, only: & - mesh_element - use crystallite, only: & - crystallite_P, & - crystallite_dPdF, & - crystallite_partionedF,& - crystallite_partionedF0 use homogenization_mech_RGC, only: & homogenization_RGC_updateState use thermal_adiabatic, only: & @@ -769,10 +695,6 @@ end function updateState !> @brief derive average stress and stiffness from constituent quantities !-------------------------------------------------------------------------------------------------- subroutine averageStressAndItsTangent(ip,el) - use mesh, only: & - mesh_element - use crystallite, only: & - crystallite_P,crystallite_dPdF use homogenization_mech_RGC, only: & homogenization_RGC_averageStressAndItsTangent @@ -810,8 +732,6 @@ end subroutine averageStressAndItsTangent !> if homogenization_sizePostResults(i,e) > 0 !! !-------------------------------------------------------------------------------------------------- function postResults(ip,el) - use mesh, only: & - mesh_element use thermal_adiabatic, only: & thermal_adiabatic_postResults use thermal_conduction, only: & @@ -868,14 +788,9 @@ end function postResults !-------------------------------------------------------------------------------------------------- subroutine homogenization_results #if defined(PETSc) || defined(DAMASK_HDF5) - use results - use homogenization_mech_RGC - use HDF5_utilities use config, only: & config_name_homogenization => homogenization_name ! anticipate logical name - use material, only: & - homogenization_typeInstance, & material_homogenization_type => homogenization_type integer :: p diff --git a/src/homogenization_mech_RGC.f90 b/src/homogenization_mech_RGC.f90 index d7b1b31bf..f27733eb5 100644 --- a/src/homogenization_mech_RGC.f90 +++ b/src/homogenization_mech_RGC.f90 @@ -8,10 +8,18 @@ !-------------------------------------------------------------------------------------------------- module homogenization_mech_RGC use prec + use IO + use config + use debug + use math use material + use numerics + use constitutive +#if defined(PETSc) || defined(DAMASK_HDF5) + use results +#endif implicit none - private enum, bind(c) enumerator :: & @@ -66,34 +74,12 @@ module homogenization_mech_RGC type(tRGCdependentState), dimension(:), allocatable :: & dependentState - public :: & - homogenization_RGC_init, & - homogenization_RGC_partitionDeformation, & - homogenization_RGC_averageStressAndItsTangent, & - homogenization_RGC_updateState, & - mech_RGC_results ! name suited for planned submodule situation - contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all necessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- -subroutine homogenization_RGC_init() - use debug, only: & -#ifdef DEBUG - debug_i, & - debug_e, & -#endif - debug_level, & - debug_homogenization, & - debug_levelBasic - use math, only: & - math_EulerToR, & - INRAD - use IO, only: & - IO_error - use config, only: & - config_homogenization +subroutine homogenization_RGC_init integer :: & Ninstance, & @@ -218,12 +204,6 @@ end subroutine homogenization_RGC_init !> @brief partitions the deformation gradient onto the constituents !-------------------------------------------------------------------------------------------------- subroutine homogenization_RGC_partitionDeformation(F,avgF,instance,of) -#ifdef DEBUG - use debug, only: & - debug_level, & - debug_homogenization, & - debug_levelExtensive -#endif real(pReal), dimension (:,:,:), intent(out) :: F !< partioned F per grain @@ -275,24 +255,6 @@ end subroutine homogenization_RGC_partitionDeformation ! "happy" with result !-------------------------------------------------------------------------------------------------- function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) -#ifdef DEBUG - use debug, only: & - debug_level, & - debug_homogenization,& - debug_levelExtensive -#endif - use math, only: & - math_invert2 - use numerics, only: & - absTol_RGC, & - relTol_RGC, & - absMax_RGC, & - relMax_RGC, & - pPert_RGC, & - maxdRelax_RGC, & - viscPower_RGC, & - viscModus_RGC, & - refRelaxRate_RGC real(pReal), dimension(:,:,:), intent(in) :: & P,& !< array of P @@ -712,10 +674,6 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) !> @brief calculate stress-like penalty due to deformation mismatch !-------------------------------------------------------------------------------------------------- subroutine stressPenalty(rPen,nMis,avgF,fDef,ip,el,instance,of) - use math, only: & - math_civita - use numerics, only: & - xSmoo_RGC real(pReal), dimension (:,:,:), intent(out) :: rPen !< stress-like penalty real(pReal), dimension (:,:), intent(out) :: nMis !< total amount of mismatch @@ -828,13 +786,6 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) !> @brief calculate stress-like penalty due to volume discrepancy !-------------------------------------------------------------------------------------------------- subroutine volumePenalty(vPen,vDiscrep,fAvg,fDef,nGrain,instance,of) - use math, only: & - math_det33, & - math_inv33 - use numerics, only: & - maxVolDiscr_RGC,& - volDiscrMod_RGC,& - volDiscrPow_RGC real(pReal), dimension (:,:,:), intent(out) :: vPen ! stress-like penalty due to volume real(pReal), intent(out) :: vDiscrep ! total volume discrepancy @@ -883,8 +834,6 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) ! deformation !-------------------------------------------------------------------------------------------------- function surfaceCorrection(avgF,instance,of) - use math, only: & - math_invert33 real(pReal), dimension(3) :: surfaceCorrection @@ -916,8 +865,6 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) !> @brief compute the equivalent shear and bulk moduli from the elasticity tensor !-------------------------------------------------------------------------------------------------- function equivalentModuli(grainID,ip,el) - use constitutive, only: & - constitutive_homogenizedC real(pReal), dimension(2) :: equivalentModuli @@ -1015,8 +962,6 @@ end subroutine homogenization_RGC_averageStressAndItsTangent !-------------------------------------------------------------------------------------------------- subroutine mech_RGC_results(instance,group) #if defined(PETSc) || defined(DAMASK_HDF5) - use results, only: & - results_writeDataset integer, intent(in) :: instance character(len=*) :: group diff --git a/src/homogenization_mech_isostrain.f90 b/src/homogenization_mech_isostrain.f90 index 7dd7bad7d..b17085b57 100644 --- a/src/homogenization_mech_isostrain.f90 +++ b/src/homogenization_mech_isostrain.f90 @@ -5,7 +5,10 @@ !> @brief Isostrain (full constraint Taylor assuption) homogenization scheme !-------------------------------------------------------------------------------------------------- submodule(homogenization) homogenization_mech_isostrain - + use config + use debug + use IO + implicit none enum, bind(c) @@ -30,14 +33,6 @@ contains !> @brief allocates all neccessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- module subroutine mech_isostrain_init - use debug, only: & - debug_HOMOGENIZATION, & - debug_level, & - debug_levelBasic - use IO, only: & - IO_error - use config, only: & - config_homogenization integer :: & Ninstance, & diff --git a/src/homogenization_mech_none.f90 b/src/homogenization_mech_none.f90 index e7a5a12e6..b8b74c267 100644 --- a/src/homogenization_mech_none.f90 +++ b/src/homogenization_mech_none.f90 @@ -5,6 +5,8 @@ !> @brief dummy homogenization homogenization scheme for 1 constituent per material point !-------------------------------------------------------------------------------------------------- submodule(homogenization) homogenization_mech_none + use config + use debug implicit none @@ -14,12 +16,6 @@ contains !> @brief allocates all neccessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- module subroutine mech_none_init - use debug, only: & - debug_HOMOGENIZATION, & - debug_level, & - debug_levelBasic - use config, only: & - config_homogenization integer :: & Ninstance, & From 2258bfb22161c60987043ef89bc2d4e4dfeac132 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 18 May 2019 07:23:46 +0200 Subject: [PATCH 14/23] RGC as submodule submodules inherit use-associated entities and implicit none/private statements --- src/homogenization.f90 | 71 ++++++++++++++++++++------- src/homogenization_mech_RGC.f90 | 70 +++++++++++--------------- src/homogenization_mech_isostrain.f90 | 5 -- src/homogenization_mech_none.f90 | 6 +-- 4 files changed, 83 insertions(+), 69 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index d7b5b3fdc..cbc6471b4 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -60,11 +60,24 @@ module homogenization module subroutine mech_isostrain_init end subroutine mech_isostrain_init + module subroutine mech_RGC_init + end subroutine mech_RGC_init + + module subroutine mech_isostrain_partitionDeformation(F,avgF) real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point end subroutine mech_isostrain_partitionDeformation + module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of) + real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient + real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point + integer, intent(in) :: & + instance, & + of + end subroutine mech_RGC_partitionDeformation + + module subroutine mech_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point @@ -73,7 +86,37 @@ module homogenization real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses integer, intent(in) :: instance end subroutine mech_isostrain_averageStressAndItsTangent + + module subroutine mech_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) + real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point + real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point + + real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses + real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses + integer, intent(in) :: instance + end subroutine mech_RGC_averageStressAndItsTangent + + + module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) +logical, dimension(2) :: mech_RGC_updateState + real(pReal), dimension(:,:,:), intent(in) :: & + P,& !< array of P + F,& !< array of F + F0 !< array of initial F + real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF !< array of current grain stiffness + real(pReal), dimension(3,3), intent(in) :: avgF !< average F + real(pReal), intent(in) :: dt !< time increment + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + end function + + module subroutine mech_RGC_results(instance,group) + + integer, intent(in) :: instance + character(len=*), intent(in) :: group + end subroutine end interface public :: & @@ -112,7 +155,7 @@ subroutine homogenization_init if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call mech_none_init if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call mech_isostrain_init - if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call homogenization_RGC_init + if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call mech_RGC_init if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init if (any(thermal_type == THERMAL_adiabatic_ID)) call thermal_adiabatic_init @@ -610,8 +653,6 @@ end subroutine materialpoint_postResults !> @brief partition material point def grad onto constituents !-------------------------------------------------------------------------------------------------- subroutine partitionDeformation(ip,el) - use homogenization_mech_RGC, only: & - homogenization_RGC_partitionDeformation integer, intent(in) :: & ip, & !< integration point @@ -628,7 +669,7 @@ subroutine partitionDeformation(ip,el) materialpoint_subF(1:3,1:3,ip,el)) case (HOMOGENIZATION_RGC_ID) chosenHomogenization - call homogenization_RGC_partitionDeformation(& + call mech_RGC_partitionDeformation(& crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & materialpoint_subF(1:3,1:3,ip,el),& ip, & @@ -643,8 +684,6 @@ end subroutine partitionDeformation !> "happy" with result !-------------------------------------------------------------------------------------------------- function updateState(ip,el) - use homogenization_mech_RGC, only: & - homogenization_RGC_updateState use thermal_adiabatic, only: & thermal_adiabatic_updateState use damage_local, only: & @@ -660,14 +699,14 @@ function updateState(ip,el) case (HOMOGENIZATION_RGC_ID) chosenHomogenization updateState = & updateState .and. & - homogenization_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & - crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & - crystallite_partionedF0(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el),& - materialpoint_subF(1:3,1:3,ip,el),& - materialpoint_subdt(ip,el), & - crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & - ip, & - el) + mech_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & + crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & + crystallite_partionedF0(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el),& + materialpoint_subF(1:3,1:3,ip,el),& + materialpoint_subdt(ip,el), & + crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & + ip, & + el) end select chosenHomogenization chosenThermal: select case (thermal_type(mesh_element(3,el))) @@ -695,8 +734,6 @@ end function updateState !> @brief derive average stress and stiffness from constituent quantities !-------------------------------------------------------------------------------------------------- subroutine averageStressAndItsTangent(ip,el) - use homogenization_mech_RGC, only: & - homogenization_RGC_averageStressAndItsTangent integer, intent(in) :: & ip, & !< integration point @@ -716,7 +753,7 @@ subroutine averageStressAndItsTangent(ip,el) homogenization_typeInstance(mesh_element(3,el))) case (HOMOGENIZATION_RGC_ID) chosenHomogenization - call homogenization_RGC_averageStressAndItsTangent(& + call mech_RGC_averageStressAndItsTangent(& materialpoint_P(1:3,1:3,ip,el), & materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),& crystallite_P(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & diff --git a/src/homogenization_mech_RGC.f90 b/src/homogenization_mech_RGC.f90 index f27733eb5..4e115a498 100644 --- a/src/homogenization_mech_RGC.f90 +++ b/src/homogenization_mech_RGC.f90 @@ -6,20 +6,7 @@ !> @brief Relaxed grain cluster (RGC) homogenization scheme !> Nconstituents is defined as p x q x r (cluster) !-------------------------------------------------------------------------------------------------- -module homogenization_mech_RGC - use prec - use IO - use config - use debug - use math - use material - use numerics - use constitutive -#if defined(PETSc) || defined(DAMASK_HDF5) - use results -#endif - - implicit none +submodule(homogenization) homogenization_mech_RGC enum, bind(c) enumerator :: & @@ -79,7 +66,7 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all necessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- -subroutine homogenization_RGC_init +module subroutine mech_RGC_init integer :: & Ninstance, & @@ -197,17 +184,17 @@ subroutine homogenization_RGC_init enddo -end subroutine homogenization_RGC_init +end subroutine mech_RGC_init !-------------------------------------------------------------------------------------------------- !> @brief partitions the deformation gradient onto the constituents !-------------------------------------------------------------------------------------------------- -subroutine homogenization_RGC_partitionDeformation(F,avgF,instance,of) +module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of) real(pReal), dimension (:,:,:), intent(out) :: F !< partioned F per grain - real(pReal), dimension (:,:), intent(in) :: avgF !< averaged F + real(pReal), dimension (3,3), intent(in) :: avgF !< averaged F integer, intent(in) :: & instance, & of @@ -247,14 +234,14 @@ subroutine homogenization_RGC_partitionDeformation(F,avgF,instance,of) end associate -end subroutine homogenization_RGC_partitionDeformation +end subroutine mech_RGC_partitionDeformation !-------------------------------------------------------------------------------------------------- !> @brief update the internal state of the homogenization scheme and tell whether "done" and ! "happy" with result !-------------------------------------------------------------------------------------------------- -function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) +module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) real(pReal), dimension(:,:,:), intent(in) :: & P,& !< array of P @@ -267,8 +254,6 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) ip, & !< integration point number el !< element number - logical, dimension(2) :: homogenization_RGC_updateState - integer, dimension(4) :: intFaceN,intFaceP,faceID integer, dimension(3) :: nGDim,iGr3N,iGr3P integer :: instance,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ipert,iGrain,nGrain, of @@ -285,7 +270,7 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) #endif zeroTimeStep: if(dEq0(dt)) then - homogenization_RGC_updateState = .true. ! pretend everything is fine and return + mech_RGC_updateState = .true. ! pretend everything is fine and return return endif zeroTimeStep @@ -404,12 +389,12 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) endif #endif - homogenization_RGC_updateState = .false. + mech_RGC_updateState = .false. !-------------------------------------------------------------------------------------------------- ! If convergence reached => done and happy if (residMax < relTol_RGC*stresMax .or. residMax < absTol_RGC) then - homogenization_RGC_updateState = .true. + mech_RGC_updateState = .true. #ifdef DEBUG if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0 & .and. prm%of_debug == of) write(6,'(1x,a55,/)')'... done and happy' @@ -451,7 +436,7 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) !-------------------------------------------------------------------------------------------------- ! if residual blows-up => done but unhappy elseif (residMax > relMax_RGC*stresMax .or. residMax > absMax_RGC) then ! try to restart when residual blows up exceeding maximum bound - homogenization_RGC_updateState = [.true.,.false.] ! with direct cut-back + mech_RGC_updateState = [.true.,.false.] ! with direct cut-back #ifdef DEBUG if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0 & @@ -648,7 +633,7 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) enddo; enddo stt%relaxationVector(:,of) = relax + drelax ! Updateing the state variable for the next iteration if (any(abs(drelax) > maxdRelax_RGC)) then ! Forcing cutback when the incremental change of relaxation vector becomes too large - homogenization_RGC_updateState = [.true.,.false.] + mech_RGC_updateState = [.true.,.false.] !$OMP CRITICAL (write2out) write(6,'(1x,a,1x,i3,1x,a,1x,i3,1x,a)')'RGC_updateState: ip',ip,'| el',el,'enforces cutback' write(6,'(1x,a,1x,e15.8)')'due to large relaxation change =',maxval(abs(drelax)) @@ -935,13 +920,13 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) end subroutine grainDeformation -end function homogenization_RGC_updateState +end function mech_RGC_updateState !-------------------------------------------------------------------------------------------------- !> @brief derive average stress and stiffness from constituent quantities !-------------------------------------------------------------------------------------------------- -subroutine homogenization_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) +module subroutine mech_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point @@ -953,7 +938,7 @@ subroutine homogenization_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF, avgP = sum(P,3) /real(product(param(instance)%Nconstituents),pReal) dAvgPdAvgF = sum(dPdF,5)/real(product(param(instance)%Nconstituents),pReal) -end subroutine homogenization_RGC_averageStressAndItsTangent +end subroutine mech_RGC_averageStressAndItsTangent !-------------------------------------------------------------------------------------------------- @@ -963,8 +948,9 @@ end subroutine homogenization_RGC_averageStressAndItsTangent subroutine mech_RGC_results(instance,group) #if defined(PETSc) || defined(DAMASK_HDF5) - integer, intent(in) :: instance - character(len=*) :: group + integer, intent(in) :: instance + character(len=*), intent(in) :: group + integer :: o associate(stt => state(instance), dst => dependentState(instance), prm => param(instance)) @@ -1135,7 +1121,7 @@ integer pure function interface4to1(iFace4D, nGDim) else interface4to1 = iFace4D(4) + nGDim(3)*(iFace4D(2)-1) & + nGDim(3)*nGDim(1)*(iFace4D(3)-1) & - + (nGDim(1)-1)*nGDim(2)*nGDim(3) ! total number of interfaces normal //e1 + + (nGDim(1)-1)*nGDim(2)*nGDim(3) ! total # of interfaces normal || e1 endif case(3) @@ -1144,8 +1130,8 @@ integer pure function interface4to1(iFace4D, nGDim) else interface4to1 = iFace4D(2) + nGDim(1)*(iFace4D(3)-1) & + nGDim(1)*nGDim(2)*(iFace4D(4)-1) & - + (nGDim(1)-1)*nGDim(2)*nGDim(3) & ! total number of interfaces normal //e1 - + nGDim(1)*(nGDim(2)-1)*nGDim(3) ! total number of interfaces normal //e2 + + (nGDim(1)-1)*nGDim(2)*nGDim(3) & ! total # of interfaces normal || e1 + + nGDim(1)*(nGDim(2)-1)*nGDim(3) ! total # of interfaces normal || e2 endif case default @@ -1169,23 +1155,23 @@ pure function interface1to4(iFace1D, nGDim) !-------------------------------------------------------------------------------------------------- ! compute the total number of interfaces, which ... - nIntFace = [(nGDim(1)-1)*nGDim(2)*nGDim(3), & ! ... normal //e1 - nGDim(1)*(nGDim(2)-1)*nGDim(3), & ! ... normal //e2 - nGDim(1)*nGDim(2)*(nGDim(3)-1)] ! ... normal //e3 + nIntFace = [(nGDim(1)-1)*nGDim(2)*nGDim(3), & ! ... normal || e1 + nGDim(1)*(nGDim(2)-1)*nGDim(3), & ! ... normal || e2 + nGDim(1)*nGDim(2)*(nGDim(3)-1)] ! ... normal || e3 !-------------------------------------------------------------------------------------------------- ! get the corresponding interface ID in 4D (normal and local position) - if (iFace1D > 0 .and. iFace1D <= nIntFace(1)) then ! interface with normal //e1 + if (iFace1D > 0 .and. iFace1D <= nIntFace(1)) then ! interface with normal || e1 interface1to4(1) = 1 interface1to4(3) = mod((iFace1D-1),nGDim(2))+1 interface1to4(4) = mod(int(real(iFace1D-1,pReal)/real(nGDim(2),pReal)),nGDim(3))+1 interface1to4(2) = int(real(iFace1D-1,pReal)/real(nGDim(2),pReal)/real(nGDim(3),pReal))+1 - elseif (iFace1D > nIntFace(1) .and. iFace1D <= (nIntFace(2) + nIntFace(1))) then ! interface with normal //e2 + elseif (iFace1D > nIntFace(1) .and. iFace1D <= (nIntFace(2) + nIntFace(1))) then ! interface with normal || e2 interface1to4(1) = 2 interface1to4(4) = mod((iFace1D-nIntFace(1)-1),nGDim(3))+1 interface1to4(2) = mod(int(real(iFace1D-nIntFace(1)-1,pReal)/real(nGDim(3),pReal)),nGDim(1))+1 interface1to4(3) = int(real(iFace1D-nIntFace(1)-1,pReal)/real(nGDim(3),pReal)/real(nGDim(1),pReal))+1 - elseif (iFace1D > nIntFace(2) + nIntFace(1) .and. iFace1D <= (nIntFace(3) + nIntFace(2) + nIntFace(1))) then ! interface with normal //e3 + elseif (iFace1D > nIntFace(2) + nIntFace(1) .and. iFace1D <= (nIntFace(3) + nIntFace(2) + nIntFace(1))) then ! interface with normal || e3 interface1to4(1) = 3 interface1to4(2) = mod((iFace1D-nIntFace(2)-nIntFace(1)-1),nGDim(1))+1 interface1to4(3) = mod(int(real(iFace1D-nIntFace(2)-nIntFace(1)-1,pReal)/real(nGDim(1),pReal)),nGDim(2))+1 @@ -1195,4 +1181,4 @@ pure function interface1to4(iFace1D, nGDim) end function interface1to4 -end module homogenization_mech_RGC +end submodule homogenization_mech_RGC diff --git a/src/homogenization_mech_isostrain.f90 b/src/homogenization_mech_isostrain.f90 index b17085b57..c1eba821c 100644 --- a/src/homogenization_mech_isostrain.f90 +++ b/src/homogenization_mech_isostrain.f90 @@ -5,11 +5,6 @@ !> @brief Isostrain (full constraint Taylor assuption) homogenization scheme !-------------------------------------------------------------------------------------------------- submodule(homogenization) homogenization_mech_isostrain - use config - use debug - use IO - - implicit none enum, bind(c) enumerator :: & diff --git a/src/homogenization_mech_none.f90 b/src/homogenization_mech_none.f90 index b8b74c267..d5b24e5d1 100644 --- a/src/homogenization_mech_none.f90 +++ b/src/homogenization_mech_none.f90 @@ -5,10 +5,6 @@ !> @brief dummy homogenization homogenization scheme for 1 constituent per material point !-------------------------------------------------------------------------------------------------- submodule(homogenization) homogenization_mech_none - use config - use debug - - implicit none contains @@ -23,7 +19,7 @@ module subroutine mech_none_init NofMyHomog write(6,'(/,a)') ' <<<+- homogenization_'//HOMOGENIZATION_NONE_label//' init -+>>>' - + Ninstance = count(homogenization_type == HOMOGENIZATION_NONE_ID) if (iand(debug_level(debug_HOMOGENIZATION),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance From 34bcd382403d1dff1c7e136d53db3e9672cb269f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 18 May 2019 07:39:55 +0200 Subject: [PATCH 15/23] cleanup --- src/commercialFEM_fileList.f90 | 2 +- src/homogenization.f90 | 53 ++++++++++++--------------------- src/homogenization_mech_RGC.f90 | 8 ++--- 3 files changed, 24 insertions(+), 39 deletions(-) diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index 0ae1323f0..5e1cd71eb 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -46,7 +46,6 @@ #include "plastic_nonlocal.f90" #include "constitutive.f90" #include "crystallite.f90" -#include "homogenization_mech_RGC.f90" #include "thermal_isothermal.f90" #include "thermal_adiabatic.f90" #include "thermal_conduction.f90" @@ -56,4 +55,5 @@ #include "homogenization.f90" #include "homogenization_mech_none.f90" #include "homogenization_mech_isostrain.f90" +#include "homogenization_mech_RGC.f90" #include "CPFEM.f90" diff --git a/src/homogenization.f90 b/src/homogenization.f90 index cbc6471b4..be07318d7 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -35,7 +35,6 @@ module homogenization materialpoint_results !< results array of material point integer, public, protected :: & materialpoint_sizeResults, & - homogenization_maxSizePostResults, & thermal_maxSizePostResults, & damage_maxSizePostResults @@ -98,25 +97,25 @@ module homogenization module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) -logical, dimension(2) :: mech_RGC_updateState - real(pReal), dimension(:,:,:), intent(in) :: & - P,& !< array of P - F,& !< array of F - F0 !< array of initial F - real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF !< array of current grain stiffness - real(pReal), dimension(3,3), intent(in) :: avgF !< average F - real(pReal), intent(in) :: dt !< time increment - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - end function - - module subroutine mech_RGC_results(instance,group) + logical, dimension(2) :: mech_RGC_updateState + real(pReal), dimension(:,:,:), intent(in) :: & + P,& !< partitioned stresses + F,& !< partitioned deformation gradients + F0 !< partitioned initial deformation gradients + real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses + real(pReal), dimension(3,3), intent(in) :: avgF !< average F + real(pReal), intent(in) :: dt !< time increment + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + end function mech_RGC_updateState - integer, intent(in) :: instance - character(len=*), intent(in) :: group - end subroutine + module subroutine mech_RGC_results(instance,group) + integer, intent(in) :: instance !< homogenization instance + character(len=*), intent(in) :: group !< group name in HDF5 file + end subroutine mech_RGC_results + end interface public :: & @@ -124,11 +123,6 @@ logical, dimension(2) :: mech_RGC_updateState materialpoint_stressAndItsTangent, & materialpoint_postResults, & homogenization_results - private :: & - partitionDeformation, & - updateState, & - averageStressAndItsTangent, & - postResults contains @@ -137,7 +131,6 @@ contains !> @brief module initialization !-------------------------------------------------------------------------------------------------- subroutine homogenization_init - use homogenization_mech_RGC use thermal_isothermal use thermal_adiabatic use thermal_conduction @@ -245,7 +238,7 @@ subroutine homogenization_init ! allocate and initialize global variables allocate(materialpoint_dPdF(3,3,3,3,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) allocate(materialpoint_F0(3,3,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) - materialpoint_F0 = spread(spread(math_I3,3,theMesh%elem%nIPs),4,theMesh%nElems) ! initialize to identity + materialpoint_F0 = spread(spread(math_I3,3,theMesh%elem%nIPs),4,theMesh%nElems) ! initialize to identity allocate(materialpoint_F(3,3,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) materialpoint_F = materialpoint_F0 ! initialize to identity allocate(materialpoint_subF0(3,3,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) @@ -260,18 +253,15 @@ subroutine homogenization_init !-------------------------------------------------------------------------------------------------- ! allocate and initialize global state and postresutls variables - homogenization_maxSizePostResults = 0 thermal_maxSizePostResults = 0 damage_maxSizePostResults = 0 do p = 1,size(config_homogenization) - homogenization_maxSizePostResults = max(homogenization_maxSizePostResults,homogState (p)%sizePostResults) thermal_maxSizePostResults = max(thermal_maxSizePostResults, thermalState (p)%sizePostResults) damage_maxSizePostResults = max(damage_maxSizePostResults ,damageState (p)%sizePostResults) enddo materialpoint_sizeResults = 1 & ! grain count - + 1 + homogenization_maxSizePostResults & ! homogSize & homogResult - + thermal_maxSizePostResults & + + 1 + thermal_maxSizePostResults & + damage_maxSizePostResults & + homogenization_maxNgrains * (1 + crystallite_maxSizePostResults & ! crystallite size & crystallite results + 1 + constitutive_plasticity_maxSizePostResults & ! constitutive size & constitutive results @@ -281,11 +271,6 @@ subroutine homogenization_init write(6,'(/,a)') ' <<<+- homogenization init -+>>>' if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0) then -#ifdef TODO - write(6,'(a32,1x,7(i8,1x))') 'homogenization_state0: ', shape(homogenization_state0) - write(6,'(a32,1x,7(i8,1x))') 'homogenization_subState0: ', shape(homogenization_subState0) - write(6,'(a32,1x,7(i8,1x))') 'homogenization_state: ', shape(homogenization_state) -#endif write(6,'(a32,1x,7(i8,1x))') 'materialpoint_dPdF: ', shape(materialpoint_dPdF) write(6,'(a32,1x,7(i8,1x))') 'materialpoint_F0: ', shape(materialpoint_F0) write(6,'(a32,1x,7(i8,1x))') 'materialpoint_F: ', shape(materialpoint_F) diff --git a/src/homogenization_mech_RGC.f90 b/src/homogenization_mech_RGC.f90 index 4e115a498..a1e92d575 100644 --- a/src/homogenization_mech_RGC.f90 +++ b/src/homogenization_mech_RGC.f90 @@ -30,7 +30,7 @@ submodule(homogenization) homogenization_mech_RGC angles integer :: & of_debug = 0 - integer(kind(undefined_ID)), dimension(:), allocatable :: & + integer(kind(undefined_ID)), dimension(:), allocatable :: & outputID end type tParameters @@ -945,7 +945,7 @@ end subroutine mech_RGC_averageStressAndItsTangent !> @brief writes results to HDF5 output file ! ToDo: check wheter units are correct !-------------------------------------------------------------------------------------------------- -subroutine mech_RGC_results(instance,group) +module subroutine mech_RGC_results(instance,group) #if defined(PETSc) || defined(DAMASK_HDF5) integer, intent(in) :: instance @@ -981,8 +981,8 @@ subroutine mech_RGC_results(instance,group) end associate #else - integer, intent(in) :: instance - character(len=*) :: group + integer, intent(in) :: instance + character(len=*), intent(in) :: group #endif end subroutine mech_RGC_results From 86205f508149e90c0e00f85af0b235c16aa0bdb7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 18 May 2019 09:47:20 +0200 Subject: [PATCH 16/23] gfortran complaints about repeated dimension attribute we need to decide whether we want to repeat the declaration of the interface or not --- src/homogenization_mech_RGC.f90 | 18 ++++-------------- 1 file changed, 4 insertions(+), 14 deletions(-) diff --git a/src/homogenization_mech_RGC.f90 b/src/homogenization_mech_RGC.f90 index a1e92d575..1cbe837d2 100644 --- a/src/homogenization_mech_RGC.f90 +++ b/src/homogenization_mech_RGC.f90 @@ -241,18 +241,7 @@ end subroutine mech_RGC_partitionDeformation !> @brief update the internal state of the homogenization scheme and tell whether "done" and ! "happy" with result !-------------------------------------------------------------------------------------------------- -module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) - - real(pReal), dimension(:,:,:), intent(in) :: & - P,& !< array of P - F,& !< array of F - F0 !< array of initial F - real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF !< array of current grain stiffness - real(pReal), dimension(3,3), intent(in) :: avgF !< average F - real(pReal), intent(in) :: dt !< time increment - integer, intent(in) :: & - ip, & !< integration point number - el !< element number +module procedure mech_RGC_updateState integer, dimension(4) :: intFaceN,intFaceP,faceID integer, dimension(3) :: nGDim,iGr3N,iGr3P @@ -579,9 +568,10 @@ module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) !-------------------------------------------------------------------------------------------------- ! ... of the numerical viscosity traction "rmatrix" allocate(rmatrix(3*nIntFaceTot,3*nIntFaceTot),source=0.0_pReal) - forall (i=1:3*nIntFaceTot) & + do i=1,3*nIntFaceTot rmatrix(i,i) = viscModus_RGC*viscPower_RGC/(refRelaxRate_RGC*dt)* & ! tangent due to numerical viscosity traction appears (abs(drelax(i))/(refRelaxRate_RGC*dt))**(viscPower_RGC - 1.0_pReal) ! only in the main diagonal term + enddo #ifdef DEBUG if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0) then @@ -920,7 +910,7 @@ module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) end subroutine grainDeformation -end function mech_RGC_updateState +end procedure mech_RGC_updateState !-------------------------------------------------------------------------------------------------- From 60ed514e0cdc013e764256773332ad7729e807ae Mon Sep 17 00:00:00 2001 From: Test User Date: Wed, 15 May 2019 23:51:29 +0200 Subject: [PATCH 17/23] [skip ci] updated version information after successful test of v2.0.3-301-g789420c9 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index acf2aa8dc..69faed072 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-297-gae084bb2 +v2.0.3-301-g789420c9 From 60c2a5fc068393033a778bdf2a823a07bdf9dcd8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 16 May 2019 18:24:54 +0200 Subject: [PATCH 18/23] loop (forall) over integration points wrong this was done for each integration point, but this was not detected for the forall loop --- src/homogenization.f90 | 119 ++++++++++++++++++++--------------------- 1 file changed, 59 insertions(+), 60 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index be07318d7..3210f02d4 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -323,43 +323,46 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) ! initialize restoration points of ... do e = FEsolving_execElem(1),FEsolving_execElem(2) myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e); do g = 1,myNgrains + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e); + do g = 1,myNgrains + + plasticState (phaseAt(g,i,e))%partionedState0(:,phasememberAt(g,i,e)) = & + plasticState (phaseAt(g,i,e))%state0( :,phasememberAt(g,i,e)) + do mySource = 1, phase_Nsources(phaseAt(g,i,e)) + sourceState(phaseAt(g,i,e))%p(mySource)%partionedState0(:,phasememberAt(g,i,e)) = & + sourceState(phaseAt(g,i,e))%p(mySource)%state0( :,phasememberAt(g,i,e)) + enddo + + crystallite_partionedFp0(1:3,1:3,g,i,e) = crystallite_Fp0(1:3,1:3,g,i,e) + crystallite_partionedLp0(1:3,1:3,g,i,e) = crystallite_Lp0(1:3,1:3,g,i,e) + crystallite_partionedFi0(1:3,1:3,g,i,e) = crystallite_Fi0(1:3,1:3,g,i,e) + crystallite_partionedLi0(1:3,1:3,g,i,e) = crystallite_Li0(1:3,1:3,g,i,e) + crystallite_partionedF0(1:3,1:3,g,i,e) = crystallite_F0(1:3,1:3,g,i,e) + crystallite_partionedS0(1:3,1:3,g,i,e) = crystallite_S0(1:3,1:3,g,i,e) - plasticState (phaseAt(g,i,e))%partionedState0(:,phasememberAt(g,i,e)) = & - plasticState (phaseAt(g,i,e))%state0( :,phasememberAt(g,i,e)) - do mySource = 1, phase_Nsources(phaseAt(g,i,e)) - sourceState(phaseAt(g,i,e))%p(mySource)%partionedState0(:,phasememberAt(g,i,e)) = & - sourceState(phaseAt(g,i,e))%p(mySource)%state0( :,phasememberAt(g,i,e)) enddo - crystallite_partionedFp0(1:3,1:3,g,i,e) = crystallite_Fp0(1:3,1:3,g,i,e) ! ...plastic def grads - crystallite_partionedLp0(1:3,1:3,g,i,e) = crystallite_Lp0(1:3,1:3,g,i,e) ! ...plastic velocity grads - crystallite_partionedFi0(1:3,1:3,g,i,e) = crystallite_Fi0(1:3,1:3,g,i,e) ! ...intermediate def grads - crystallite_partionedLi0(1:3,1:3,g,i,e) = crystallite_Li0(1:3,1:3,g,i,e) ! ...intermediate velocity grads - crystallite_partionedF0(1:3,1:3,g,i,e) = crystallite_F0(1:3,1:3,g,i,e) ! ...def grads - crystallite_partionedS0(1:3,1:3,g,i,e) = crystallite_S0(1:3,1:3,g,i,e) ! ...2nd PK stress - enddo; enddo - forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e)) - materialpoint_subF0(1:3,1:3,i,e) = materialpoint_F0(1:3,1:3,i,e) ! ...def grad + materialpoint_subF0(1:3,1:3,i,e) = materialpoint_F0(1:3,1:3,i,e) materialpoint_subFrac(i,e) = 0.0_pReal materialpoint_subStep(i,e) = 1.0_pReal/subStepSizeHomog ! <> materialpoint_converged(i,e) = .false. ! pretend failed step of twice the required size materialpoint_requested(i,e) = .true. ! everybody requires calculation - endforall - forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), & - homogState(material_homogenizationAt(e))%sizeState > 0) & - homogState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) = & - homogState(material_homogenizationAt(e))%State0( :,mappingHomogenization(1,i,e)) ! ...internal homogenization state - forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), & - thermalState(material_homogenizationAt(e))%sizeState > 0) & - thermalState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) = & - thermalState(material_homogenizationAt(e))%State0( :,mappingHomogenization(1,i,e)) ! ...internal thermal state - forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), & - damageState(material_homogenizationAt(e))%sizeState > 0) & - damageState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) = & - damageState(material_homogenizationAt(e))%State0( :,mappingHomogenization(1,i,e)) ! ...internal damage state + + if (homogState(material_homogenizationAt(e))%sizeState > 0) & + homogState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) = & + homogState(material_homogenizationAt(e))%State0( :,mappingHomogenization(1,i,e)) ! ...internal homogenization state + + if (thermalState(material_homogenizationAt(e))%sizeState > 0) & + thermalState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) = & + thermalState(material_homogenizationAt(e))%State0( :,mappingHomogenization(1,i,e)) ! ...internal thermal state + + if (damageState(material_homogenizationAt(e))%sizeState > 0) & + damageState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) = & + damageState(material_homogenizationAt(e))%State0( :,mappingHomogenization(1,i,e)) ! ...internal damage state + enddo enddo + NiterationHomog = 0 cutBackLooping: do while (.not. terminallyIll .and. & @@ -370,7 +373,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) myNgrains = homogenization_Ngrains(mesh_element(3,e)) IpLooping1: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - converged: if ( materialpoint_converged(i,e) ) then + converged: if (materialpoint_converged(i,e)) then #ifdef DEBUG if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0 & .and. ((e == debug_e .and. i == debug_i) & @@ -391,22 +394,22 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) ! wind forward grain starting point of... crystallite_partionedF0 (1:3,1:3,1:myNgrains,i,e) = & - crystallite_partionedF(1:3,1:3,1:myNgrains,i,e) ! ...def grads + crystallite_partionedF(1:3,1:3,1:myNgrains,i,e) crystallite_partionedFp0 (1:3,1:3,1:myNgrains,i,e) = & - crystallite_Fp (1:3,1:3,1:myNgrains,i,e) ! ...plastic def grads + crystallite_Fp (1:3,1:3,1:myNgrains,i,e) crystallite_partionedLp0 (1:3,1:3,1:myNgrains,i,e) = & - crystallite_Lp (1:3,1:3,1:myNgrains,i,e) ! ...plastic velocity grads + crystallite_Lp (1:3,1:3,1:myNgrains,i,e) crystallite_partionedFi0 (1:3,1:3,1:myNgrains,i,e) = & - crystallite_Fi (1:3,1:3,1:myNgrains,i,e) ! ...intermediate def grads + crystallite_Fi (1:3,1:3,1:myNgrains,i,e) crystallite_partionedLi0 (1:3,1:3,1:myNgrains,i,e) = & - crystallite_Li (1:3,1:3,1:myNgrains,i,e) ! ...intermediate velocity grads + crystallite_Li (1:3,1:3,1:myNgrains,i,e) crystallite_partionedS0 (1:3,1:3,1:myNgrains,i,e) = & - crystallite_S (1:3,1:3,1:myNgrains,i,e) ! ...2nd PK stress + crystallite_S (1:3,1:3,1:myNgrains,i,e) do g = 1,myNgrains plasticState (phaseAt(g,i,e))%partionedState0(:,phasememberAt(g,i,e)) = & @@ -417,23 +420,22 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) enddo enddo - forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), & - homogState(material_homogenizationAt(e))%sizeState > 0) & + if(homogState(material_homogenizationAt(e))%sizeState > 0) & homogState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) = & - homogState(material_homogenizationAt(e))%State (:,mappingHomogenization(1,i,e)) ! ...internal homogenization state - forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), & - thermalState(material_homogenizationAt(e))%sizeState > 0) & + homogState(material_homogenizationAt(e))%State (:,mappingHomogenization(1,i,e)) + if(thermalState(material_homogenizationAt(e))%sizeState > 0) & thermalState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) = & - thermalState(material_homogenizationAt(e))%State (:,mappingHomogenization(1,i,e)) ! ...internal thermal state - forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), & - damageState(material_homogenizationAt(e))%sizeState > 0) & + thermalState(material_homogenizationAt(e))%State (:,mappingHomogenization(1,i,e)) + if(damageState(material_homogenizationAt(e))%sizeState > 0) & damageState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) = & - damageState(material_homogenizationAt(e))%State (:,mappingHomogenization(1,i,e)) ! ...internal damage state - materialpoint_subF0(1:3,1:3,i,e) = materialpoint_subF(1:3,1:3,i,e) ! ...def grad + damageState(material_homogenizationAt(e))%State (:,mappingHomogenization(1,i,e)) + + materialpoint_subF0(1:3,1:3,i,e) = materialpoint_subF(1:3,1:3,i,e) + endif steppingNeeded else converged - if ( (myNgrains == 1 .and. materialpoint_subStep(i,e) <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite + if ( (myNgrains == 1 .and. materialpoint_subStep(i,e) <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite subStepSizeHomog * materialpoint_subStep(i,e) <= subStepMinHomog ) then ! would require too small subStep ! cutback makes no sense !$OMP FLUSH(terminallyIll) @@ -462,16 +464,16 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) ! restore... if (materialpoint_subStep(i,e) < 1.0_pReal) then ! protect against fake cutback from \Delta t = 2 to 1. Maybe that "trick" is not necessary anymore at all? I.e. start with \Delta t = 1 crystallite_Lp(1:3,1:3,1:myNgrains,i,e) = & - crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e) ! ...plastic velocity grads + crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e) crystallite_Li(1:3,1:3,1:myNgrains,i,e) = & - crystallite_partionedLi0(1:3,1:3,1:myNgrains,i,e) ! ...intermediate velocity grads + crystallite_partionedLi0(1:3,1:3,1:myNgrains,i,e) endif ! maybe protecting everything from overwriting (not only L) makes even more sense crystallite_Fp(1:3,1:3,1:myNgrains,i,e) = & - crystallite_partionedFp0(1:3,1:3,1:myNgrains,i,e) ! ...plastic def grads + crystallite_partionedFp0(1:3,1:3,1:myNgrains,i,e) crystallite_Fi(1:3,1:3,1:myNgrains,i,e) = & - crystallite_partionedFi0(1:3,1:3,1:myNgrains,i,e) ! ...intermediate def grads + crystallite_partionedFi0(1:3,1:3,1:myNgrains,i,e) crystallite_S(1:3,1:3,1:myNgrains,i,e) = & - crystallite_partionedS0(1:3,1:3,1:myNgrains,i,e) ! ...2nd PK stress + crystallite_partionedS0(1:3,1:3,1:myNgrains,i,e) do g = 1, myNgrains plasticState (phaseAt(g,i,e))%state( :,phasememberAt(g,i,e)) = & plasticState (phaseAt(g,i,e))%partionedState0(:,phasememberAt(g,i,e)) @@ -480,18 +482,15 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) sourceState(phaseAt(g,i,e))%p(mySource)%partionedState0(:,phasememberAt(g,i,e)) enddo enddo - forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), & - homogState(material_homogenizationAt(e))%sizeState > 0) & + if(homogState(material_homogenizationAt(e))%sizeState > 0) & homogState(material_homogenizationAt(e))%State( :,mappingHomogenization(1,i,e)) = & - homogState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) ! ...internal homogenization state - forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), & - thermalState(material_homogenizationAt(e))%sizeState > 0) & + homogState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) + if(thermalState(material_homogenizationAt(e))%sizeState > 0) & thermalState(material_homogenizationAt(e))%State( :,mappingHomogenization(1,i,e)) = & - thermalState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) ! ...internal thermal state - forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), & - damageState(material_homogenizationAt(e))%sizeState > 0) & + thermalState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) + if(damageState(material_homogenizationAt(e))%sizeState > 0) & damageState(material_homogenizationAt(e))%State( :,mappingHomogenization(1,i,e)) = & - damageState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) ! ...internal damage state + damageState(material_homogenizationAt(e))%subState0(:,mappingHomogenization(1,i,e)) endif endif converged From c7036e39701c3e7b0fcbf4cfc80db47cc76ac92c Mon Sep 17 00:00:00 2001 From: Test User Date: Fri, 17 May 2019 09:06:28 +0200 Subject: [PATCH 19/23] [skip ci] updated version information after successful test of v2.0.3-304-g7b14263c --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 69faed072..2d8c83361 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-301-g789420c9 +v2.0.3-304-g7b14263c From 1de0c7e652fcc6572859976e645374b6a36bd1fb Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 18 May 2019 14:19:31 +0200 Subject: [PATCH 20/23] wrong name --- src/plastic_dislotwin.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/plastic_dislotwin.f90 b/src/plastic_dislotwin.f90 index 7db1f5f7f..05e9a49fd 100644 --- a/src/plastic_dislotwin.f90 +++ b/src/plastic_dislotwin.f90 @@ -1068,7 +1068,7 @@ subroutine plastic_dislotwin_results(instance,group) call results_writeDataset(group,stt%rho_dip,'rho_dip',& 'dislocation dipole density''1/m²') case (dot_gamma_sl_ID) - call results_writeDataset(group,stt%gamma_sl,'dot_gamma_sl',& + call results_writeDataset(group,stt%gamma_sl,'gamma_sl',& 'plastic shear','1') case (Lambda_sl_ID) call results_writeDataset(group,dst%Lambda_sl,'Lambda_sl',& From c4f07a9ad97102f9e8065dcdab06a1283298fbf3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 18 May 2019 21:25:05 +0200 Subject: [PATCH 21/23] need to correct tensor order --- src/results.f90 | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/src/results.f90 b/src/results.f90 index c5582b927..05db831f7 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -299,18 +299,26 @@ end subroutine results_writeVectorDataset_real !-------------------------------------------------------------------------------------------------- subroutine results_writeTensorDataset_real(group,dataset,label,description,SIunit) - character(len=*), intent(in) :: label,group,description - character(len=*), intent(in), optional :: SIunit - real(pReal), intent(inout), dimension(:,:,:) :: dataset + character(len=*), intent(in) :: label,group,description + character(len=*), intent(in), optional :: SIunit + real(pReal), intent(in), dimension(:,:,:) :: dataset + integer :: i integer(HID_T) :: groupHandle - + real(pReal), dimension(:,:,:), allocatable :: dataset_transposed + + + allocate(dataset_transposed,mold=dataset) + do i=1,size(dataset,3) + dataset_transposed(1:3,1:3,i) = transpose(dataset(1:3,1:3,i)) + enddo + groupHandle = results_openGroup(group) #ifdef PETSc - call HDF5_write(groupHandle,dataset,label,.true.) + call HDF5_write(groupHandle,dataset_transposed,label,.true.) #else - call HDF5_write(groupHandle,dataset,label,.false.) + call HDF5_write(groupHandle,dataset_transposed,label,.false.) #endif if (HDF5_objectExists(groupHandle,label)) & From d2e64df6a1813b4b169fdf3fbeb3d401ab7d69a1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 18 May 2019 21:27:27 +0200 Subject: [PATCH 22/23] some tests based on HDF5 output --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 639c6f4a5..3a2f89547 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 639c6f4a5eafc893c83c740c57f417eaaabc45ae +Subproject commit 3a2f89547c264044a7bfab9d33aee78eec495a76 From b35465b591f4d4866f421045d03d1ac7da9a3a99 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 20 May 2019 19:08:56 +0200 Subject: [PATCH 23/23] gamma_slip_ID should be used to write result --- src/plastic_dislotwin.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/plastic_dislotwin.f90 b/src/plastic_dislotwin.f90 index 05e9a49fd..cb4b3cbae 100644 --- a/src/plastic_dislotwin.f90 +++ b/src/plastic_dislotwin.f90 @@ -1067,7 +1067,7 @@ subroutine plastic_dislotwin_results(instance,group) case (rho_dip_ID) call results_writeDataset(group,stt%rho_dip,'rho_dip',& 'dislocation dipole density''1/m²') - case (dot_gamma_sl_ID) + case (gamma_sl_ID) call results_writeDataset(group,stt%gamma_sl,'gamma_sl',& 'plastic shear','1') case (Lambda_sl_ID)