diff --git a/VERSION b/VERSION index 3a0d0e38b..bbcfb4711 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v3.0.0-alpha2-530-g0d0226f70 +v3.0.0-alpha2-602-ge2d4ab427 diff --git a/examples/ConfigFiles/Microstructure_DP_Steel.config b/examples/ConfigFiles/Microstructure_DP_Steel.config deleted file mode 100644 index 6c083e2c7..000000000 --- a/examples/ConfigFiles/Microstructure_DP_Steel.config +++ /dev/null @@ -1,4 +0,0 @@ -[DP_Steel] -crystallite 1 -(constituent) phase 1 texture 1 fraction 0.82 -(constituent) phase 2 texture 2 fraction 0.18 diff --git a/examples/ConfigFiles/Phase_Dislotwin_TWIP-Steel-FeMnC.config b/examples/ConfigFiles/Phase_Dislotwin_TWIP-Steel-FeMnC.config deleted file mode 100644 index 3ca635d73..000000000 --- a/examples/ConfigFiles/Phase_Dislotwin_TWIP-Steel-FeMnC.config +++ /dev/null @@ -1,64 +0,0 @@ -[TWIP_Steel_FeMnC] - -elasticity hooke -plasticity dislotwin - -(output) rho_mob -(output) rho_dip -(output) gamma_sl -(output) lambda_sl -(output) tau_pass -(output) f_tw -(output) lambda_tw -(output) tau_hat_tw -(output) f_tr - - -### Material parameters ### -lattice_structure fcc -C11 175.0e9 # From Music et al. Applied Physics Letters 91, 191904 (2007) -C12 115.0e9 -C44 135.0e9 -grainsize 2.0e-5 # Average grain size [m] -SolidSolutionStrength 1.5e8 # Strength due to elements in solid solution - -### Dislocation glide parameters ### -Nslip 12 -slipburgers 2.56e-10 # Burgers vector of slip system [m] -rhoedgedip0 1.0 # Initial dislocation density [m/m**3] -rhoedge0 1.0e12 # Initial dislocation density [m/m**3] -v0 1.0e-4 # Initial glide velocity [m/s] -Qedge 3.7e-19 # Activation energy for dislocation glide [J] -p_slip 1.0 # p-exponent in glide velocity -q_slip 1.0 # q-exponent in glide velocity - -# hardening of glide -CLambdaSlip 10.0 # Adj. parameter controlling dislocation mean free path -D0 4.0e-5 # Vacancy diffusion prefactor [m**2/s] -Qsd 4.5e-19 # Activation energy for climb [J] -Catomicvolume 1.0 # Adj. parameter controlling the atomic volume [in b^3] -Cedgedipmindistance 1.0 # Adj. parameter controlling the minimum dipole distance [in b] -interactionSlipSlip 0.122 0.122 0.625 0.07 0.137 0.122 # Interaction coefficients (Kubin et al. 2008) - -### Shearband parameters ### -shearbandresistance 180e6 -shearbandvelocity 0e-4 # set to zero to turn shear banding of -QedgePerSbSystem 3.7e-19 # Activation energy for shear banding [J] -p_shearband 1.0 # p-exponent in glide velocity -q_shearband 1.0 # q-exponent in glide velocity - -### Twinning parameters ### -Ntwin 12 -twinburgers 1.47e-10 # Burgers vector of twin system [m] -twinsize 5.0e-8 # Twin stack mean thickness [m] -L0_twin 442.0 # Length of twin nuclei in Burgers vectors -maxtwinfraction 1.0 # Maximum admissible twin volume fraction -xc_twin 1.0e-9 # critical distance for formation of twin nucleus -VcrossSlip 1.67e-29 # cross slip volume -r_twin 10.0 # r-exponent in twin formation probability -Cmfptwin 1.0 # Adj. parameter controlling twin mean free path -Cthresholdtwin 1.0 # Adj. parameter controlling twin threshold stress -interactionSlipTwin 0.0 1.0 1.0 # Dislocation-Twin interaction coefficients -interactionTwinTwin 0.0 1.0 # Twin-Twin interaction coefficients -SFE_0K -0.0396 # stacking fault energy at zero K; TWIP steel: -0.0526; Cu: -0.0396 -dSFE_dT 0.0002 # temperature dependance of stacking fault energy diff --git a/examples/ConfigFiles/Phase_Dislotwin_TWIP-Steel-FeMnC.yaml b/examples/ConfigFiles/Phase_Dislotwin_TWIP-Steel-FeMnC.yaml new file mode 100644 index 000000000..009443284 --- /dev/null +++ b/examples/ConfigFiles/Phase_Dislotwin_TWIP-Steel-FeMnC.yaml @@ -0,0 +1,41 @@ +TWIP_Steel_FeMnC: + lattice: cF + mechanics: + elasticity: {type: hooke, C_11: 175.0e9, C_12: 115.0e9, C_44: 135.0e9} + plasticity: + type: dislotwin + output: [rho_mob, rho_dip, gamma_sl, Lambda_sl, tau_pass, f_tw, Lambda_tw, tau_hat_tw, f_tr] + D: 2.0e-5 + N_sl: [12] + b_sl: [2.56e-10] + rho_mob_0: [1.0e12] + rho_dip_0: [1.0] + v_0: [1.0e4] + Q_s: [3.7e-19] + p_sl: [1.0] + q_sl: [1.0] + tau_0: [1.5e8] + i_sl: [10.0] # Adj. parameter controlling dislocation mean free path + D_0: 4.0e-5 # Vacancy diffusion prefactor / m^2/s + D_a: 1.0 # minimum dipole distance / b + Q_cl: 4.5e-19 # Activation energy for climb / J + h_sl_sl: [0.122, 0.122, 0.625, 0.07, 0.137, 0.122] # Interaction coefficients (Kubin et al. 2008) +# shear band parameters + xi_sb: 180.0e6 + Q_sb: 3.7e-19 + p_sb: 1.0 + q_sb: 1.0 + v_sb: 0.0 # set to 0, to turn it off +# twinning parameters + N_tw: [12] + b_tw: [1.47e-10] # Burgers vector length of twin system / b + t_tw: [5.0e-8] # Twin stack mean thickness / m + L_tw: 442.0 # Length of twin nuclei / b + x_c_tw: 1.0e-9 # critical distance for formation of twin nucleus / m + V_cs: 1.67e-29 # cross slip volume / m^3 + p_tw: [10.0] # r-exponent in twin formation probability + i_tw: 1.0 # Adj. parameter controlling twin mean free path + h_sl_tw: [0.0, 1.0, 1.0] # dislocation-twin interaction coefficients + h_tw_tw: [0.0, 1.0] # twin-twin interaction coefficients + Gamma_sf_0K: -0.0396 # stacking fault energy / J/m^2 at zero K; TWIP steel: -0.0526; Cu: -0.0396 + dGamma_sf_dT: 0.0002 # temperature dependence / J/(m^2 K) of stacking fault energy diff --git a/examples/ConfigFiles/Phase_Dislotwin_Tungsten.config b/examples/ConfigFiles/Phase_Dislotwin_Tungsten.config deleted file mode 100644 index 30c04cb9a..000000000 --- a/examples/ConfigFiles/Phase_Dislotwin_Tungsten.config +++ /dev/null @@ -1,36 +0,0 @@ -[Tungsten] - -elasticity hooke -plasticity dislotwin - -### Material parameters ### -lattice_structure bcc -C11 523.0e9 # From Marinica et al. Journal of Physics: Condensed Matter(2013) -C12 202.0e9 -C44 161.0e9 - -grainsize 2.0e-5 # Average grain size [m] -SolidSolutionStrength 1.5e8 # Strength due to elements in solid solution - -### Dislocation glide parameters ### -#per family -Nslip 12 -slipburgers 2.72e-10 # Burgers vector of slip system [m] -rhoedge0 1.0e12 # Initial edge dislocation density [m/m**3] -rhoedgedip0 1.0 # Initial edged dipole dislocation density [m/m**3] -v0 1.0e-4 # Initial glide velocity [m/s] -Qedge 2.725e-19 # Activation energy for dislocation glide [J] -p_slip 0.78 # p-exponent in glide velocity -q_slip 1.58 # q-exponent in glide velocity -tau_peierls 2.03e9 # peierls stress (for bcc) -dipoleformationfactor 0 # to have hardening due to dipole formation off - -#hardening -CLambdaSlip 10.0 # Adj. parameter controlling dislocation mean free path -D0 4.0e-5 # Vacancy diffusion prefactor [m**2/s] -Qsd 4.5e-19 # Activation energy for climb [J] -Catomicvolume 1.0 # Adj. parameter controlling the atomic volume [in b] -Cedgedipmindistance 1.0 # Adj. parameter controlling the minimum dipole distance [in b] -interaction_slipslip 1 1 1.4 1.4 1.4 1.4 - - diff --git a/examples/ConfigFiles/Phase_Dislotwin_Tungsten.yaml b/examples/ConfigFiles/Phase_Dislotwin_Tungsten.yaml new file mode 100644 index 000000000..36467192b --- /dev/null +++ b/examples/ConfigFiles/Phase_Dislotwin_Tungsten.yaml @@ -0,0 +1,21 @@ +Tungsten: + lattice: cI + mechanics: + elasticity: {type: hooke, C_11: 523.0e9, C_12: 202.0e9, C_44: 161.0e9} # Marinica et al. Journal of Physics: Condensed Matter(2013) + plasticity: + type: dislotwin + D: 2.0e-5 # Average grain size / m + N_sl: [12] + b_sl: [2.72e-10] # Burgers vector length of slip families / m + rho_mob_0: [1.0e12] + rho_dip_0: [1.0] + v_0: [1.0e4] # Initial glide velocity / m/s + Q_s: [2.725e-19] # Activation energy for dislocation glide / J + p_sl: [0.78] # p-exponent in glide velocity + q_sl: [1.58] # q-exponent in glide velocity + tau_0: [1.5e8] # solid solution strength / Pa + i_sl: [10.0] # Adj. parameter controlling dislocation mean free path + D_0: 4.0e-5 # Vacancy diffusion prefactor / m^2/s + D_a: 1.0 # minimum dipole distance / b + Q_cl: 4.5e-19 # Activation energy for climb / J + h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4] diff --git a/examples/ConfigFiles/Phase_Hydrogen.config b/examples/ConfigFiles/Phase_Hydrogen.config deleted file mode 100644 index c9ccfdc86..000000000 --- a/examples/ConfigFiles/Phase_Hydrogen.config +++ /dev/null @@ -1,3 +0,0 @@ -hydrogenflux_diffusion11 1.0 -hydrogenflux_mobility11 1.0 -hydrogenVolume 1e-28 diff --git a/examples/ConfigFiles/Phase_Phenopowerlaw_cpTi.yaml b/examples/ConfigFiles/Phase_Phenopowerlaw_cpTi.yaml index 7931ec267..aa5262454 100644 --- a/examples/ConfigFiles/Phase_Phenopowerlaw_cpTi.yaml +++ b/examples/ConfigFiles/Phase_Phenopowerlaw_cpTi.yaml @@ -1,13 +1,14 @@ # M. Levy, Handbook of Elastic Properties of Solids, Liquids, and Gases (2001) # C. Zambaldi, "Orientation informed nanoindentation of a-titanium: Indentation pileup in hexagonal metals deforming by prismatic slip", J. Mater. Res., Vol. 27, No. 1, Jan 14, 2012 -Ti-alpha: +# Better use values from L. Wang, Z. Zheng, H. Phukan, P. Kenesei, J.-S. Park, J. Lind, R.M. Suter, T.R. Bieler, Direct measurement of critical resolved shear stress of prismatic and basal slip in polycrystalline Ti using high energy X-ray diffraction microscopy, Acta Mater 2017 +cpTi: lattice: hP c/a: 1.587 mechanics: output: [F, P, F_e, F_p, L_p, O] elasticity: {C_11: 160.0e9, C_12: 90.0e9, C_13: 66.0e9, C_33: 181.7e9, C_44: 46.5e9, type: hooke} plasticity: - N_sl: [3, 3, 0, 0, 12] + N_sl: [3, 3, 0, 6, 12] a_sl: 2.0 dot_gamma_0_sl: 0.001 h_0_sl_sl: 200e6 @@ -15,5 +16,5 @@ Ti-alpha: n_sl: 20 output: [gamma_sl] type: phenopowerlaw - xi_0_sl: [349e6, 150e6, 0, 0, 1107e6] - xi_inf_sl: [568e6, 1502e6, 0, 0, 3420e6] + xi_0_sl: [0.15e9, 0.09e9, 0, 0.20e9, 0.25e9] + xi_inf_sl: [0.24e9, 0.5e9, 0, 0.6e9, 0.8e9] diff --git a/python/damask/_config.py b/python/damask/_config.py index edf3a5676..85f0c208c 100644 --- a/python/damask/_config.py +++ b/python/damask/_config.py @@ -1,5 +1,6 @@ import copy from io import StringIO +from collections.abc import Iterable import abc import numpy as np @@ -44,6 +45,42 @@ class Config(dict): copy = __copy__ + def __or__(self,other): + """ + Update configuration with contents of other. + + Parameters + ---------- + other : damask.Config or dict + Key-value pairs that update self. + + """ + duplicate = self.copy() + duplicate.update(other) + return duplicate + + + def __ior__(self,other): + """Update configuration with contents of other.""" + return self.__or__(other) + + + def delete(self,keys): + """ + Remove configuration keys. + + Parameters + ---------- + keys : iterable or scalar + Label of the key(s) to remove. + + """ + duplicate = self.copy() + for k in keys if isinstance(keys, Iterable) and not isinstance(keys, str) else [keys]: + del duplicate[k] + return duplicate + + @classmethod def load(cls,fname): """ @@ -99,30 +136,6 @@ class Config(dict): fhandle.write(yaml.dump(self,Dumper=NiceDumper,**kwargs)) - def add(self,d): - """ - Add dictionary. - - d : dict - Dictionary to append. - """ - duplicate = self.copy() - duplicate.update(d) - return duplicate - - - def delete(self,key): - """ - Delete item. - - key : str or scalar - Label of the key to remove. - """ - duplicate = self.copy() - del duplicate[key] - return duplicate - - @property @abc.abstractmethod def is_complete(self): diff --git a/python/damask/_result.py b/python/damask/_result.py index 6f2494299..4ecd0ba51 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -105,7 +105,7 @@ class Result: self.view('increments',visible_increments) in_between = '' if len(visible_increments) < 3 else \ - ''.join([f'\n{inc}\n ...\n' for inc in visible_increments[1:-2]]) + ''.join([f'\n{inc}\n ...\n' for inc in visible_increments[1:-1]]) return util.srepr(first + in_between + last) @@ -119,10 +119,10 @@ class Result: action : str Select from 'set', 'add', and 'del'. what : str - Attribute to change (must be in self.visible). + Attribute to change (must be from self.visible). datasets : list of str or bool - Name of datasets as list, supports ? and * wildcards. - True is equivalent to [*], False is equivalent to [] + Name of datasets as list; supports ? and * wildcards. + True is equivalent to [*], False is equivalent to []. """ def natural_sort(key): @@ -200,7 +200,7 @@ class Result: self._allow_modification = True def disallow_modification(self): - """Disllow to overwrite existing data (default case).""" + """Disallow to overwrite existing data (default case).""" self._allow_modification = False @@ -272,10 +272,10 @@ class Result: Parameters ---------- what : str - attribute to change (must be from self.visible) + Attribute to change (must be from self.visible). datasets : list of str or bool - name of datasets as list, supports ? and * wildcards. - True is equivalent to [*], False is equivalent to [] + Name of datasets as list; supports ? and * wildcards. + True is equivalent to [*], False is equivalent to []. """ self._manage_view('set',what,datasets) @@ -288,10 +288,10 @@ class Result: Parameters ---------- what : str - attribute to change (must be from self.visible) + Attribute to change (must be from self.visible). datasets : list of str or bool - name of datasets as list, supports ? and * wildcards. - True is equivalent to [*], False is equivalent to [] + Name of datasets as list; supports ? and * wildcards. + True is equivalent to [*], False is equivalent to []. """ self._manage_view('add',what,datasets) @@ -304,10 +304,10 @@ class Result: Parameters ---------- what : str - attribute to change (must be from self.visible) + Attribute to change (must be from self.visible). datasets : list of str or bool - name of datasets as list, supports ? and * wildcards. - True is equivalent to [*], False is equivalent to [] + Name of datasets as list; supports ? and * wildcards. + True is equivalent to [*], False is equivalent to []. """ self._manage_view('del',what,datasets) @@ -315,14 +315,14 @@ class Result: def rename(self,name_old,name_new): """ - Rename datasets. + Rename dataset. Parameters ---------- name_old : str - name of the datasets to be renamed + Name of the dataset to be renamed. name_new : str - new name of the datasets + New name of the dataset. """ if self._allow_modification: @@ -353,13 +353,13 @@ class Result: ---------- datasets : iterable or str constituent : int - Constituent to consider for phase data + Constituent to consider for phase data. tagged : bool - tag Table.column name with '#constituent' - defaults to False + Tag Table.column name with '#constituent'. + Defaults to False. split : bool - split Table by increment and return dictionary of Tables - defaults to True + Split Table by increment and return dictionary of Tables. + Defaults to True. """ sets = datasets if hasattr(datasets,'__iter__') and not isinstance(datasets,str) else \ @@ -371,7 +371,7 @@ class Result: with h5py.File(self.fname,'r') as f: for dataset in sets: for group in self.groups_with_datasets(dataset): - path = os.path.join(group,dataset) + path = '/'.join([group,dataset]) inc,prop,name,cat,item = (path.split('/') + ['']*5)[:5] key = '/'.join([prop,name+tag]) if key not in inGeom: @@ -388,15 +388,15 @@ class Result: np.nan, dtype=np.dtype(f[path])) data[inGeom[key]] = (f[path] if len(shape)>1 else np.expand_dims(f[path],1))[inData[key]] - path = (os.path.join(*([prop,name]+([cat] if cat else [])+([item] if item else []))) if split else path)+tag + path = ('/'.join([prop,name]+([cat] if cat else [])+([item] if item else [])) if split else path)+tag if split: try: - tbl[inc].add(path,data) + tbl[inc] = tbl[inc].add(path,data) except KeyError: tbl[inc] = Table(data.reshape(self.N_materialpoints,-1),{path:data.shape[1:]}) else: try: - tbl.add(path,data) + tbl = tbl.add(path,data) except AttributeError: tbl = Table(data.reshape(self.N_materialpoints,-1),{path:data.shape[1:]}) @@ -415,7 +415,7 @@ class Result: are considered as they contain user-relevant data. Single strings will be treated as list with one entry. - Wild card matching is allowed, but the number of arguments need to fit. + Wild card matching is allowed, but the number of arguments needs to fit. Parameters ---------- diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 86311546e..34a5f81b7 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -475,15 +475,17 @@ class Rotation: Parameters ---------- - vector : bool, optional - Return as actual Rodrigues-Frank vector, i.e. axis - and angle argument are not separated. + compact : bool, optional + Return as actual Rodrigues-Frank vector, + i.e. axis and angle argument are not separated. Returns ------- - rho : numpy.ndarray of shape (...,4) unless vector == True: - numpy.ndarray of shape (...,3) - Rodrigues-Frank vector: [n_1, n_2, n_3, tan(ω/2)], ǀnǀ = 1 and ω ∈ [0,π]. + rho : numpy.ndarray of shape (...,4) containing + [n_1, n_2, n_3, tan(ω/2)], ǀnǀ = 1 and ω ∈ [0,π] + unless compact == True: + numpy.ndarray of shape (...,3) containing + tan(ω/2) [n_1, n_2, n_3], ω ∈ [0,π]. """ ro = Rotation._qu2ro(self.quaternion) diff --git a/python/damask/_table.py b/python/damask/_table.py index 68d6f94bf..ee64ba017 100644 --- a/python/damask/_table.py +++ b/python/damask/_table.py @@ -189,6 +189,11 @@ class Table: label : str Column label. + Returns + ------- + data : numpy.ndarray + Array of column data. + """ if re.match(r'[0-9]*?_',label): idx,key = label.split('_',1) @@ -212,6 +217,11 @@ class Table: info : str, optional Human-readable information about the new data. + Returns + ------- + table : Table + Updated table. + """ dup = self.copy() dup._add_comment(label,data.shape[1:],info) @@ -238,6 +248,11 @@ class Table: info : str, optional Human-readable information about the modified data. + Returns + ------- + table : Table + Updated table. + """ dup = self.copy() dup._add_comment(label,data.shape[1:],info) @@ -261,6 +276,11 @@ class Table: label : str Column label. + Returns + ------- + table : Table + Updated table. + """ dup = self.copy() dup.data.drop(columns=label,inplace=True) @@ -279,6 +299,11 @@ class Table: label_new : str or iterable of str New column label(s). + Returns + ------- + table : Table + Updated table. + """ dup = self.copy() columns = dict(zip([old] if isinstance(old,str) else old, @@ -300,6 +325,11 @@ class Table: ascending : bool or list, optional Set sort order. + Returns + ------- + table : Table + Updated table. + """ dup = self.copy() dup._label_discrete() @@ -320,6 +350,11 @@ class Table: other : Table Table to append. + Returns + ------- + table : Table + Concatenated table. + """ if self.shapes != other.shapes or not self.data.columns.equals(other.data.columns): raise KeyError('Labels or shapes or order do not match') @@ -340,6 +375,11 @@ class Table: other : Table Table to join. + Returns + ------- + table : Table + Joined table. + """ if set(self.shapes) & set(other.shapes) or self.data.shape[0] != other.data.shape[0]: raise KeyError('Dublicated keys or row count mismatch') diff --git a/python/damask/util.py b/python/damask/util.py index cda532bc0..3722fe33f 100644 --- a/python/damask/util.py +++ b/python/damask/util.py @@ -133,6 +133,8 @@ def execute(cmd, stdout = stdout.decode('utf-8').replace('\x08','') stderr = stderr.decode('utf-8').replace('\x08','') if process.returncode != 0: + print(stdout) + print(stderr) raise RuntimeError(f"'{cmd}' failed with returncode {process.returncode}") return stdout, stderr @@ -193,7 +195,7 @@ def scale_to_coprime(v): return m -def project_stereographic(vector,normalize=False): +def project_stereographic(vector,direction='z',normalize=True,keepdims=False): """ Apply stereographic projection to vector. @@ -201,18 +203,37 @@ def project_stereographic(vector,normalize=False): ---------- vector : numpy.ndarray of shape (...,3) Vector coordinates to be projected. + direction : str + Projection direction 'x', 'y', or 'z'. + Defaults to 'z'. normalize : bool - Ensure unit length for vector. Defaults to False. + Ensure unit length of input vector. Defaults to True. + keepdims : bool + Maintain three-dimensional output coordinates. + Default two-dimensional output uses right-handed frame spanned by + the next and next-next axis relative to the projection direction, + e.g. x-y when projecting along z and z-x when projecting along y. Returns ------- - coordinates : numpy.ndarray of shape (...,2) + coordinates : numpy.ndarray of shape (...,2 | 3) Projected coordinates. + Examples + -------- + >>> project_stereographic(np.ones(3)) + [0.3660254, 0.3660254] + >>> project_stereographic(np.ones(3),direction='x',normalize=False,keepdims=True) + [0, 0.5, 0.5] + >>> project_stereographic([0,1,1],direction='y',normalize=True,keepdims=False) + [0.41421356, 0] + """ - v_ = vector/np.linalg.norm(vector,axis=-1,keepdims=True) if normalize else vector - return np.block([v_[...,:2]/(1+np.abs(v_[...,2:3])), - np.zeros_like(v_[...,2:3])]) + shift = 'zyx'.index(direction) + v_ = np.roll(vector/np.linalg.norm(vector,axis=-1,keepdims=True) if normalize else vector, + shift,axis=-1) + return np.roll(np.block([v_[...,:2]/(1+np.abs(v_[...,2:3])),np.zeros_like(v_[...,2:3])]), + -shift if keepdims else 0,axis=-1)[...,:3 if keepdims else 2] def execution_stamp(class_name,function_name=None): @@ -418,7 +439,7 @@ class _ProgressBar: bar = '█' * filled_length + '░' * (self.bar_length - filled_length) delta_time = datetime.datetime.now() - self.start_time remaining_time = (self.total - (iteration+1)) * delta_time / (iteration+1) - remaining_time -= datetime.timedelta(microseconds=remaining_time.microseconds) # remove μs + remaining_time -= datetime.timedelta(microseconds=remaining_time.microseconds) # remove μs sys.stderr.write(f'\r{self.prefix} {bar} {fraction:>4.0%} ETA {remaining_time}') sys.stderr.flush() diff --git a/python/tests/test_Config.py b/python/tests/test_Config.py index 0319fb6de..9324c28c5 100644 --- a/python/tests/test_Config.py +++ b/python/tests/test_Config.py @@ -23,8 +23,17 @@ class TestConfig: assert Config.load(f) == config def test_add_remove(self): + dummy = {'hello':'world','foo':'bar'} config = Config() - assert config.add({'hello':'world'}).delete('hello') == config + config |= dummy + assert config == Config() | dummy + config = config.delete(dummy) + assert config == Config() + assert (config | dummy ).delete( 'hello' ) == config | {'foo':'bar'} + assert (config | dummy ).delete([ 'hello', 'foo' ]) == config + assert (config | Config(dummy)).delete({ 'hello':1,'foo':2 }) == config + assert (config | Config(dummy)).delete(Config({'hello':1 })) == config | {'foo':'bar'} + def test_repr(self,tmp_path): config = Config() diff --git a/python/tests/test_util.py b/python/tests/test_util.py index eb1084b09..397926682 100644 --- a/python/tests/test_util.py +++ b/python/tests/test_util.py @@ -49,17 +49,18 @@ class TestUtil: dist_sampled = np.histogram(centers[selected],bins)[0]/N_samples*np.sum(dist) assert np.sqrt(((dist - dist_sampled) ** 2).mean()) < .025 and selected.shape[0]==N_samples - @pytest.mark.parametrize('point,normalize,answer', + @pytest.mark.parametrize('point,direction,normalize,keepdims,answer', [ - ([1,0,0],False,[1,0,0]), - ([1,0,0],True, [1,0,0]), - ([0,1,1],False,[0,0.5,0]), - ([0,1,1],True, [0,0.41421356,0]), - ([1,1,1],False,[0.5,0.5,0]), - ([1,1,1],True, [0.3660254, 0.3660254, 0]), + ([1,0,0],'z',False,True, [1,0,0]), + ([1,0,0],'z',True, False,[1,0]), + ([0,1,1],'z',False,True, [0,0.5,0]), + ([0,1,1],'y',True, False,[0.41421356,0]), + ([1,1,0],'x',False,False,[0.5,0]), + ([1,1,1],'y',True, True, [0.3660254, 0,0.3660254]), ]) - def test_project_stereographic(self,point,normalize,answer): - assert np.allclose(util.project_stereographic(np.array(point),normalize=normalize),answer) + def test_project_stereographic(self,point,direction,normalize,keepdims,answer): + assert np.allclose(util.project_stereographic(np.array(point),direction=direction, + normalize=normalize,keepdims=keepdims),answer) @pytest.mark.parametrize('fro,to,mode,answer', [ 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/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 diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index d250e2f53..0aa462e7e 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -37,9 +37,10 @@ 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 real(pReal), dimension(9) :: temp_valueVector !< temporarily from loadcase file when reading in tensors (initialize to 0.0) @@ -53,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 = 300.0_pReal, & time = 0.0_pReal, & !< elapsed time time0 = 0.0_pReal, & !< begin of interval timeinc = 1.0_pReal, & !< current time interval @@ -78,8 +80,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 +99,13 @@ program DAMASK_grid quit class (tNode), pointer :: & num_grid, & - debug_grid, & ! pointer to grid debug options config_load, & load_steps, & load_step, & + solver, & + initial_conditions, & + ic_thermal, & + thermal, & step_bc, & step_mech, & step_discretization, & @@ -112,17 +116,11 @@ 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' -!-------------------------------------------------------------------------------------------------- -! 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 @@ -133,19 +131,22 @@ 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,32 +161,36 @@ 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 (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)) + + field = 1 + ID(field) = FIELD_MECH_ID ! mechanical active by default + thermalActive: if (solver%get_asString('thermal',defaultVal = 'n/a') == 'spectral') then + field = field + 1 + ID(field) = FIELD_THERMAL_ID + endif thermalActive + damageActive: if (solver%get_asString('damage',defaultVal = 'n/a') == 'spectral') then + field = field + 1 + ID(field) = FIELD_DAMAGE_ID + endif damageActive + !-------------------------------------------------------------------------------------------------- -! 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 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') @@ -220,19 +225,17 @@ 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 - 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 @@ -286,13 +289,13 @@ 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, ext_msg = loadcase_string) ! exit with error message + if (errorID > 0) call IO_error(error_ID = errorID, el = l) endif reportAndCheck enddo @@ -301,12 +304,14 @@ 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 case(FIELD_THERMAL_ID) - call grid_thermal_spectral_init + 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 @@ -377,7 +382,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, & @@ -397,7 +402,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) diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index afbae027f..8d3f913fa 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -160,7 +160,7 @@ function grid_damage_spectral_solution(timeinc) result(solution) real(pReal), intent(in) :: & timeinc !< increment in time for current solution - integer :: i, j, k, cell + integer :: i, j, k, ce type(tSolutionState) :: solution PetscInt :: devNull PetscReal :: phi_min, phi_max, stagNorm, solnNorm @@ -193,10 +193,10 @@ function grid_damage_spectral_solution(timeinc) result(solution) !-------------------------------------------------------------------------------------------------- ! updating damage state - cell = 0 + ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 - call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),1,cell) + ce = ce + 1 + call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),ce) enddo; enddo; enddo call VecMin(solution_vec,devNull,phi_min,ierr); CHKERRQ(ierr) @@ -216,7 +216,7 @@ end function grid_damage_spectral_solution subroutine grid_damage_spectral_forward(cutBack) logical, intent(in) :: cutBack - integer :: i, j, k, cell + integer :: i, j, k, ce DM :: dm_local PetscScalar, dimension(:,:,:), pointer :: x_scal PetscErrorCode :: ierr @@ -226,14 +226,14 @@ subroutine grid_damage_spectral_forward(cutBack) phi_stagInc = phi_lastInc !-------------------------------------------------------------------------------------------------- ! reverting damage field state - cell = 0 + ce = 0 call SNESGetDM(damage_snes,dm_local,ierr); CHKERRQ(ierr) call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with x_scal(xstart:xend,ystart:yend,zstart:zend) = phi_current call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 - call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),1,cell) + ce = ce + 1 + call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),ce) enddo; enddo; enddo else phi_lastInc = phi_current @@ -258,7 +258,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) f_scal PetscObject :: dummy PetscErrorCode :: ierr - integer :: i, j, k, cell + integer :: i, j, k, ce real(pReal) :: phiDot, dPhiDot_dPhi, mobility phi_current = x_scal @@ -269,20 +269,20 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) call utilities_FFTscalarForward call utilities_fourierScalarGradient !< calculate gradient of damage field call utilities_FFTvectorBackward - cell = 0 + ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 - vectorField_real(1:3,i,j,k) = matmul(damage_nonlocal_getDiffusion(1,cell) - K_ref, & + ce = ce + 1 + vectorField_real(1:3,i,j,k) = matmul(damage_nonlocal_getDiffusion(ce) - K_ref, & vectorField_real(1:3,i,j,k)) enddo; enddo; enddo call utilities_FFTvectorForward call utilities_fourierVectorDivergence !< calculate damage divergence in fourier field call utilities_FFTscalarBackward - cell = 0 + ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 - call damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi_current(i,j,k), 1, cell) - mobility = damage_nonlocal_getMobility(1,cell) + ce = ce + 1 + call damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi_current(i,j,k),ce) + mobility = damage_nonlocal_getMobility(ce) scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + phiDot) & + mobility*(phi_lastInc(i,j,k) - phi_current(i,j,k)) & + mu_ref*phi_current(i,j,k) @@ -310,15 +310,15 @@ end subroutine formResidual !-------------------------------------------------------------------------------------------------- subroutine updateReference - integer :: i,j,k,cell,ierr + integer :: i,j,k,ce,ierr - cell = 0 + ce = 0 K_ref = 0.0_pReal mu_ref = 0.0_pReal do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 - K_ref = K_ref + damage_nonlocal_getDiffusion(1,cell) - mu_ref = mu_ref + damage_nonlocal_getMobility(1,cell) + ce = ce + 1 + K_ref = K_ref + damage_nonlocal_getDiffusion(ce) + mu_ref = mu_ref + damage_nonlocal_getMobility(ce) enddo; enddo; enddo K_ref = K_ref*wgt call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) 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' diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 866914f5b..5ac043a67 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 @@ -131,9 +133,10 @@ subroutine grid_thermal_spectral_init 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 @@ -268,7 +271,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 - vectorField_real(1:3,i,j,k) = matmul(thermal_conduction_getConductivity(1,ce) - K_ref, & + vectorField_real(1:3,i,j,k) = matmul(thermal_conduction_getConductivity(ce) - K_ref, & vectorField_real(1:3,i,j,k)) enddo; enddo; enddo call utilities_FFTvectorForward @@ -277,7 +280,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 - call thermal_conduction_getSource(Tdot, 1,ce) + call thermal_conduction_getSource(Tdot,1,ce) scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + Tdot) & + thermal_conduction_getMassDensity (ce)* & thermal_conduction_getSpecificHeat(ce)*(T_lastInc(i,j,k) - & @@ -310,7 +313,7 @@ subroutine updateReference mu_ref = 0.0_pReal do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 - K_ref = K_ref + thermal_conduction_getConductivity(1,ce) + K_ref = K_ref + thermal_conduction_getConductivity(ce) mu_ref = mu_ref + thermal_conduction_getMassDensity(ce)* thermal_conduction_getSpecificHeat(ce) enddo; enddo; enddo K_ref = K_ref*wgt diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 1de8bb93a..f6ca104ae 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -35,19 +35,11 @@ module homogenization homogState, & damageState_h - integer, dimension(:), allocatable, public, protected :: & - homogenization_typeInstance, & !< instance of particular type of each homogenization - thermal_typeInstance, & !< instance of particular type of each thermal transport - damage_typeInstance !< instance of particular type of each nonlocal damage - - 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 :: & + integer(kind(DAMAGE_none_ID)), dimension(:), allocatable :: & 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 @@ -95,12 +87,11 @@ module homogenization module subroutine damage_init end subroutine damage_init - module subroutine mechanical_partition(subF,ip,el) + module subroutine mechanical_partition(subF,ce) real(pReal), intent(in), dimension(3,3) :: & subF integer, intent(in) :: & - ip, & !< integration point - el !< element number + ce end subroutine mechanical_partition module subroutine thermal_partition(ce) @@ -115,37 +106,31 @@ module homogenization integer, intent(in) :: ip,el end subroutine thermal_homogenize - module subroutine mechanical_homogenize(dt,ip,el) + module subroutine mechanical_homogenize(dt,ce) real(pReal), intent(in) :: dt integer, intent(in) :: & - ip, & !< integration point - el !< element number + ce !< cell end subroutine mechanical_homogenize - module subroutine mechanical_results(group_base,h) + module subroutine mechanical_results(group_base,ho) character(len=*), intent(in) :: group_base - integer, intent(in) :: h + integer, intent(in) :: ho end subroutine mechanical_results - module function mechanical_updateState(subdt,subF,ip,el) result(doneAndHappy) + module function mechanical_updateState(subdt,subF,ce) result(doneAndHappy) real(pReal), intent(in) :: & - subdt !< current time step + subdt !< current time step real(pReal), intent(in), dimension(3,3) :: & subF integer, intent(in) :: & - ip, & !< integration point - el !< element number + ce !< cell logical, dimension(2) :: doneAndHappy end function mechanical_updateState - module function thermal_conduction_getConductivity(ip,el) result(K) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number + 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) @@ -173,48 +158,37 @@ module homogenization real(pReal) :: T end function homogenization_thermal_T - module subroutine thermal_conduction_getSource(Tdot, ip,el) + module subroutine thermal_conduction_getSource(Tdot, ip, el) integer, intent(in) :: & - ip, & !< integration point number - el !< element number + ip, & + el real(pReal), intent(out) :: Tdot end subroutine thermal_conduction_getSource - module function damage_nonlocal_getMobility(ip,el) result(M) - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - 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, ip, el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number + module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ce) + integer, intent(in) :: ce real(pReal), intent(in) :: & phi real(pReal) :: & phiDot, dPhiDot_dPhi end subroutine damage_nonlocal_getSourceAndItsTangent - - module subroutine damage_nonlocal_putNonLocalDamage(phi,ip,el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number + 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(homog,group) - - integer, intent(in) :: homog + module subroutine damage_nonlocal_results(ho,group) + integer, intent(in) :: ho character(len=*), intent(in) :: group - end subroutine damage_nonlocal_results + end interface public :: & @@ -257,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 @@ -318,7 +289,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE if (.not. doneAndHappy(1)) then - call mechanical_partition(homogenization_F(1:3,1:3,ce),ip,el) + call mechanical_partition(homogenization_F(1:3,1:3,ce),ce) converged = .true. do co = 1, myNgrains converged = converged .and. crystallite_stress(dt,co,ip,el) @@ -327,7 +298,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE if (.not. converged) then doneAndHappy = [.true.,.false.] else - doneAndHappy = mechanical_updateState(dt,homogenization_F(1:3,1:3,ce),ip,el) + doneAndHappy = mechanical_updateState(dt,homogenization_F(1:3,1:3,ce),ce) converged = all(doneAndHappy) endif endif @@ -341,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 @@ -355,21 +326,22 @@ 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 - !$OMP DO PRIVATE(ho) + !$OMP DO PRIVATE(ho,ce) elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2) ho = material_homogenizationAt(el) IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2) + ce = (el-1)*discretization_nIPs + ip do co = 1, homogenization_Nconstituents(ho) call crystallite_orientations(co,ip,el) enddo - call mechanical_homogenize(dt,ip,el) + call mechanical_homogenize(dt,ce) enddo IpLooping3 enddo elementLooping3 !$OMP END DO @@ -527,26 +499,25 @@ end subroutine damage_nonlocal_init !-------------------------------------------------------------------------------------------------- !> @brief returns homogenized non local damage diffusion tensor in reference configuration !-------------------------------------------------------------------------------------------------- -function damage_nonlocal_getDiffusion(ip,el) +function damage_nonlocal_getDiffusion(ce) - integer, intent(in) :: & - ip, & !< integration point number - el !< element number + integer, intent(in) :: ce real(pReal), dimension(3,3) :: & damage_nonlocal_getDiffusion integer :: & - homog, & - grain + ho, & + co - homog = material_homogenizationAt(el) + ho = material_homogenizationAt2(ce) damage_nonlocal_getDiffusion = 0.0_pReal - do grain = 1, homogenization_Nconstituents(homog) + + do co = 1, homogenization_Nconstituents(ho) damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + & - crystallite_push33ToRef(grain,ip,el,lattice_D(1:3,1:3,material_phaseAt(grain,el))) + crystallite_push33ToRef(co,ce,lattice_D(1:3,1:3,material_phaseAt2(co,ce))) enddo damage_nonlocal_getDiffusion = & - num_damage%charLength**2*damage_nonlocal_getDiffusion/real(homogenization_Nconstituents(homog),pReal) + num_damage%charLength**2*damage_nonlocal_getDiffusion/real(homogenization_Nconstituents(ho),pReal) end function damage_nonlocal_getDiffusion @@ -568,13 +539,9 @@ 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(homogenization_typeInstance(size(material_name_homogenization)), source=0) - allocate(thermal_typeInstance(size(material_name_homogenization)), source=0) - allocate(damage_typeInstance(size(material_name_homogenization)), source=0) - 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) do h=1, size(material_name_homogenization) homog => material_homogenization%get(h) @@ -590,28 +557,20 @@ subroutine material_parseHomogenization call IO_error(500,ext_msg=homogMech%get_asString('type')) end select - homogenization_typeInstance(h) = count(homogenization_type==homogenization_type(h)) - - if(homog%contains('thermal')) then + 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('isothermal') - thermal_type(h) = THERMAL_isothermal_ID - case('conduction') + case('pass') thermal_type(h) = THERMAL_conduction_ID case default call IO_error(500,ext_msg=homogThermal%get_asString('type')) end select endif - if(homog%contains('damage')) then + 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')) @@ -619,12 +578,6 @@ subroutine material_parseHomogenization endif enddo - do h=1, size(material_name_homogenization) - homogenization_typeInstance(h) = count(homogenization_type(1:h) == homogenization_type(h)) - thermal_typeInstance(h) = count(thermal_type (1:h) == thermal_type (h)) - damage_typeInstance(h) = count(damage_type (1:h) == damage_type (h)) - enddo - end subroutine material_parseHomogenization diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90 index 502174309..993ea84f6 100644 --- a/src/homogenization_damage.f90 +++ b/src/homogenization_damage.f90 @@ -85,22 +85,20 @@ end subroutine damage_partition !-------------------------------------------------------------------------------------------------- !> @brief Returns homogenized nonlocal damage mobility !-------------------------------------------------------------------------------------------------- -module function damage_nonlocal_getMobility(ip,el) result(M) +module function damage_nonlocal_getMobility(ce) result(M) - integer, intent(in) :: & - ip, & !< integration point number - el !< element number + integer, intent(in) :: ce integer :: & co real(pReal) :: M M = 0.0_pReal - do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - M = M + lattice_M(material_phaseAt(co,el)) + do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) + M = M + lattice_M(material_phaseAt2(co,ce)) enddo - M = M/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) + M = M/real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) end function damage_nonlocal_getMobility @@ -108,11 +106,9 @@ end function damage_nonlocal_getMobility !-------------------------------------------------------------------------------------------------- !> @brief calculates homogenized damage driving forces !-------------------------------------------------------------------------------------------------- -module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el) +module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ce) - integer, intent(in) :: & - ip, & !< integration point number - el !< element number + integer, intent(in) :: ce real(pReal), intent(in) :: & phi real(pReal) :: & @@ -121,9 +117,9 @@ module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, p phiDot = 0.0_pReal dPhiDot_dPhi = 0.0_pReal - call phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) - phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) - dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) + call phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ce) + phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) + dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) end subroutine damage_nonlocal_getSourceAndItsTangent @@ -131,20 +127,18 @@ end subroutine damage_nonlocal_getSourceAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief updated nonlocal damage field with solution from damage phase field PDE !-------------------------------------------------------------------------------------------------- -module subroutine damage_nonlocal_putNonLocalDamage(phi,ip,el) +module subroutine damage_nonlocal_putNonLocalDamage(phi,ce) - integer, intent(in) :: & - ip, & !< integration point number - el !< element number + integer, intent(in) :: ce real(pReal), intent(in) :: & phi integer :: & - homog, & - offset + ho, & + me - homog = material_homogenizationAt(el) - offset = material_homogenizationMemberAt(ip,el) - damagestate_h(homog)%state(1,offset) = phi + ho = material_homogenizationAt2(ce) + me = material_homogenizationMemberAt2(ce) + damagestate_h(ho)%state(1,me) = phi end subroutine damage_nonlocal_putNonLocalDamage @@ -152,18 +146,18 @@ end subroutine damage_nonlocal_putNonLocalDamage !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -module subroutine damage_nonlocal_results(homog,group) +module subroutine damage_nonlocal_results(ho,group) - integer, intent(in) :: homog + integer, intent(in) :: ho character(len=*), intent(in) :: group integer :: o - associate(prm => param(damage_typeInstance(homog))) + associate(prm => param(ho)) outputsLoop: do o = 1,size(prm%output) select case(prm%output(o)) case ('phi') - call results_writeDataset(group,damagestate_h(homog)%state(1,:),prm%output(o),& + call results_writeDataset(group,damagestate_h(ho)%state(1,:),prm%output(o),& 'damage indicator','-') end select enddo outputsLoop diff --git a/src/homogenization_mechanical.f90 b/src/homogenization_mechanical.f90 index a19a61f72..f14ddf4db 100644 --- a/src/homogenization_mechanical.f90 +++ b/src/homogenization_mechanical.f90 @@ -24,35 +24,34 @@ submodule(homogenization) mechanical real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point end subroutine mechanical_isostrain_partitionDeformation - module subroutine mechanical_RGC_partitionDeformation(F,avgF,instance,of) + module subroutine mechanical_RGC_partitionDeformation(F,avgF,ce) 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 + ce end subroutine mechanical_RGC_partitionDeformation - module subroutine mechanical_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) + module subroutine mechanical_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho) 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 + integer, intent(in) :: ho end subroutine mechanical_isostrain_averageStressAndItsTangent - module subroutine mechanical_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) + module subroutine mechanical_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho) 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 + integer, intent(in) :: ho end subroutine mechanical_RGC_averageStressAndItsTangent - module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAndHappy) + module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) logical, dimension(2) :: doneAndHappy real(pReal), dimension(:,:,:), intent(in) :: & P,& !< partitioned stresses @@ -61,13 +60,12 @@ submodule(homogenization) mechanical 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 + ce !< cell end function mechanical_RGC_updateState - module subroutine mechanical_RGC_results(instance,group) - integer, intent(in) :: instance !< homogenization instance + module subroutine mechanical_RGC_results(ho,group) + integer, intent(in) :: ho !< homogenization type character(len=*), intent(in) :: group !< group name in HDF5 file end subroutine mechanical_RGC_results @@ -104,19 +102,18 @@ end subroutine mechanical_init !-------------------------------------------------------------------------------------------------- !> @brief Partition F onto the individual constituents. !-------------------------------------------------------------------------------------------------- -module subroutine mechanical_partition(subF,ip,el) +module subroutine mechanical_partition(subF,ce) real(pReal), intent(in), dimension(3,3) :: & subF integer, intent(in) :: & - ip, & !< integration point - el !< element number + ce integer :: co - real(pReal), dimension (3,3,homogenization_Nconstituents(material_homogenizationAt(el))) :: Fs + real(pReal), dimension (3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) :: Fs - chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) + chosenHomogenization: select case(homogenization_type(material_homogenizationAt2(ce))) case (HOMOGENIZATION_NONE_ID) chosenHomogenization Fs(1:3,1:3,1) = subF @@ -125,12 +122,12 @@ module subroutine mechanical_partition(subF,ip,el) call mechanical_isostrain_partitionDeformation(Fs,subF) case (HOMOGENIZATION_RGC_ID) chosenHomogenization - call mechanical_RGC_partitionDeformation(Fs,subF,ip,el) + call mechanical_RGC_partitionDeformation(Fs,subF,ce) end select chosenHomogenization - do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) - call phase_mechanical_setF(Fs(1:3,1:3,co),co,ip,el) + do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce)) + call phase_mechanical_setF(Fs(1:3,1:3,co),co,ce) enddo @@ -140,46 +137,43 @@ end subroutine mechanical_partition !-------------------------------------------------------------------------------------------------- !> @brief Average P and dPdF from the individual constituents. !-------------------------------------------------------------------------------------------------- -module subroutine mechanical_homogenize(dt,ip,el) +module subroutine mechanical_homogenize(dt,ce) real(pReal), intent(in) :: dt - integer, intent(in) :: & - ip, & !< integration point - el !< element number + integer, intent(in) :: ce - integer :: co,ce - real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el))) - real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationAt(el))) + integer :: co + real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) + real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) - ce = (el-1)* discretization_nIPs + ip - chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) + chosenHomogenization: select case(homogenization_type(material_homogenizationAt2(ce))) case (HOMOGENIZATION_NONE_ID) chosenHomogenization - homogenization_P(1:3,1:3,ce) = phase_mechanical_getP(1,ip,el) - homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(dt,1,ip,el) + homogenization_P(1:3,1:3,ce) = phase_mechanical_getP(1,ce) + homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(dt,1,ce) case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization - do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ip,el) - Ps(:,:,co) = phase_mechanical_getP(co,ip,el) + do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) + dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ce) + Ps(:,:,co) = phase_mechanical_getP(co,ce) enddo call mechanical_isostrain_averageStressAndItsTangent(& homogenization_P(1:3,1:3,ce), & homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& Ps,dPdFs, & - homogenization_typeInstance(material_homogenizationAt(el))) + material_homogenizationAt2(ce)) case (HOMOGENIZATION_RGC_ID) chosenHomogenization - do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ip,el) - Ps(:,:,co) = phase_mechanical_getP(co,ip,el) + do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) + dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ce) + Ps(:,:,co) = phase_mechanical_getP(co,ce) enddo call mechanical_RGC_averageStressAndItsTangent(& homogenization_P(1:3,1:3,ce), & homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& Ps,dPdFs, & - homogenization_typeInstance(material_homogenizationAt(el))) + material_homogenizationAt2(ce)) end select chosenHomogenization @@ -190,30 +184,29 @@ end subroutine mechanical_homogenize !> @brief update the internal state of the homogenization scheme and tell whether "done" and !> "happy" with result !-------------------------------------------------------------------------------------------------- -module function mechanical_updateState(subdt,subF,ip,el) result(doneAndHappy) +module function mechanical_updateState(subdt,subF,ce) result(doneAndHappy) real(pReal), intent(in) :: & subdt !< current time step real(pReal), intent(in), dimension(3,3) :: & subF integer, intent(in) :: & - ip, & !< integration point - el !< element number + ce logical, dimension(2) :: doneAndHappy integer :: co - real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el))) - real(pReal) :: Fs(3,3,homogenization_Nconstituents(material_homogenizationAt(el))) - real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationAt(el))) + real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) + real(pReal) :: Fs(3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) + real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) - if (homogenization_type(material_homogenizationAt(el)) == HOMOGENIZATION_RGC_ID) then - do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(subdt,co,ip,el) - Fs(:,:,co) = phase_mechanical_getF(co,ip,el) - Ps(:,:,co) = phase_mechanical_getP(co,ip,el) + if (homogenization_type(material_homogenizationAt2(ce)) == HOMOGENIZATION_RGC_ID) then + do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) + dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(subdt,co,ce) + Fs(:,:,co) = phase_mechanical_getF(co,ce) + Ps(:,:,co) = phase_mechanical_getP(co,ce) enddo - doneAndHappy = mechanical_RGC_updateState(Ps,Fs,subF,subdt,dPdFs,ip,el) + doneAndHappy = mechanical_RGC_updateState(Ps,Fs,subF,subdt,dPdFs,ce) else doneAndHappy = .true. endif @@ -224,20 +217,20 @@ end function mechanical_updateState !-------------------------------------------------------------------------------------------------- !> @brief Write results to file. !-------------------------------------------------------------------------------------------------- -module subroutine mechanical_results(group_base,h) +module subroutine mechanical_results(group_base,ho) character(len=*), intent(in) :: group_base - integer, intent(in) :: h + integer, intent(in) :: ho character(len=:), allocatable :: group group = trim(group_base)//'/mech' call results_closeGroup(results_addGroup(group)) - select case(homogenization_type(h)) + select case(homogenization_type(ho)) case(HOMOGENIZATION_rgc_ID) - call mechanical_RGC_results(homogenization_typeInstance(h),group) + call mechanical_RGC_results(ho,group) end select diff --git a/src/homogenization_mechanical_RGC.f90 b/src/homogenization_mechanical_RGC.f90 index f1bcda380..a68e9b772 100644 --- a/src/homogenization_mechanical_RGC.f90 +++ b/src/homogenization_mechanical_RGC.f90 @@ -77,8 +77,7 @@ module subroutine mechanical_RGC_init(num_homogMech) num_homogMech !< pointer to mechanical homogenization numerics data integer :: & - Ninstances, & - h, & + ho, & Nmaterialpoints, & sizeState, nIntFaceTot @@ -90,8 +89,7 @@ module subroutine mechanical_RGC_init(num_homogMech) print'(/,a)', ' <<<+- homogenization:mechanical:RGC init -+>>>' - Ninstances = count(homogenization_type == HOMOGENIZATION_RGC_ID) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) + print'(a,i2)', ' # instances: ',count(homogenization_type == HOMOGENIZATION_RGC_ID); flush(IO_STDOUT) print*, 'Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009' print*, 'https://doi.org/10.1007/s12289-009-0619-1'//IO_EOL @@ -101,10 +99,11 @@ module subroutine mechanical_RGC_init(num_homogMech) - allocate(param(Ninstances)) - allocate(state(Ninstances)) - allocate(state0(Ninstances)) - allocate(dependentState(Ninstances)) + material_homogenization => config_material%get('homogenization') + allocate(param(material_homogenization%length)) + allocate(state(material_homogenization%length)) + allocate(state0(material_homogenization%length)) + allocate(dependentState(material_homogenization%length)) num_RGC => num_homogMech%get('RGC',defaultVal=emptyDict) @@ -137,15 +136,14 @@ module subroutine mechanical_RGC_init(num_homogMech) if (num%volDiscrPow <= 0.0_pReal) call IO_error(301,ext_msg='volDiscrPw_RGC') - material_homogenization => config_material%get('homogenization') - do h = 1, size(homogenization_type) - if (homogenization_type(h) /= HOMOGENIZATION_RGC_ID) cycle - homog => material_homogenization%get(h) + do ho = 1, size(homogenization_type) + if (homogenization_type(ho) /= HOMOGENIZATION_RGC_ID) cycle + homog => material_homogenization%get(ho) homogMech => homog%get('mechanics') - associate(prm => param(homogenization_typeInstance(h)), & - stt => state(homogenization_typeInstance(h)), & - st0 => state0(homogenization_typeInstance(h)), & - dst => dependentState(homogenization_typeInstance(h))) + associate(prm => param(ho), & + stt => state(ho), & + st0 => state0(ho), & + dst => dependentState(ho)) #if defined (__GFORTRAN__) prm%output = output_asStrings(homogMech) @@ -154,7 +152,7 @@ module subroutine mechanical_RGC_init(num_homogMech) #endif prm%N_constituents = homogMech%get_asInts('cluster_size',requiredSize=3) - if (homogenization_Nconstituents(h) /= product(prm%N_constituents)) & + if (homogenization_Nconstituents(ho) /= product(prm%N_constituents)) & call IO_error(211,ext_msg='N_constituents (mechanical_RGC)') prm%xi_alpha = homogMech%get_asFloat('xi_alpha') @@ -163,18 +161,18 @@ module subroutine mechanical_RGC_init(num_homogMech) prm%D_alpha = homogMech%get_asFloats('D_alpha', requiredSize=3) prm%a_g = homogMech%get_asFloats('a_g', requiredSize=3) - Nmaterialpoints = count(material_homogenizationAt == h) + Nmaterialpoints = count(material_homogenizationAt == ho) nIntFaceTot = 3*( (prm%N_constituents(1)-1)*prm%N_constituents(2)*prm%N_constituents(3) & + prm%N_constituents(1)*(prm%N_constituents(2)-1)*prm%N_constituents(3) & + prm%N_constituents(1)*prm%N_constituents(2)*(prm%N_constituents(3)-1)) sizeState = nIntFaceTot - homogState(h)%sizeState = sizeState - allocate(homogState(h)%state0 (sizeState,Nmaterialpoints), source=0.0_pReal) - allocate(homogState(h)%state (sizeState,Nmaterialpoints), source=0.0_pReal) + homogState(ho)%sizeState = sizeState + allocate(homogState(ho)%state0 (sizeState,Nmaterialpoints), source=0.0_pReal) + allocate(homogState(ho)%state (sizeState,Nmaterialpoints), source=0.0_pReal) - stt%relaxationVector => homogState(h)%state(1:nIntFaceTot,:) - st0%relaxationVector => homogState(h)%state0(1:nIntFaceTot,:) + stt%relaxationVector => homogState(ho)%state(1:nIntFaceTot,:) + st0%relaxationVector => homogState(ho)%state0(1:nIntFaceTot,:) allocate(dst%volumeDiscrepancy( Nmaterialpoints), source=0.0_pReal) allocate(dst%relaxationRate_avg( Nmaterialpoints), source=0.0_pReal) @@ -183,7 +181,7 @@ module subroutine mechanical_RGC_init(num_homogMech) !-------------------------------------------------------------------------------------------------- ! assigning cluster orientations - dependentState(homogenization_typeInstance(h))%orientation = spread(eu2om(prm%a_g*inRad),3,Nmaterialpoints) + dependentState(ho)%orientation = spread(eu2om(prm%a_g*inRad),3,Nmaterialpoints) !dst%orientation = spread(eu2om(prm%a_g*inRad),3,Nmaterialpoints) ifort version 18.0.1 crashes (for whatever reason) end associate @@ -196,22 +194,23 @@ end subroutine mechanical_RGC_init !-------------------------------------------------------------------------------------------------- !> @brief partitions the deformation gradient onto the constituents !-------------------------------------------------------------------------------------------------- -module subroutine mechanical_RGC_partitionDeformation(F,avgF,instance,of) +module subroutine mechanical_RGC_partitionDeformation(F,avgF,ce) real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned F per grain real(pReal), dimension (3,3), intent(in) :: avgF !< averaged F integer, intent(in) :: & - instance, & - of + ce real(pReal), dimension(3) :: aVect,nVect integer, dimension(4) :: intFace integer, dimension(3) :: iGrain3 - integer :: iGrain,iFace,i,j - - associate(prm => param(instance)) + integer :: iGrain,iFace,i,j,ho,me + associate(prm => param(material_homogenizationAt2(ce))) + + ho = material_homogenizationAt2(ce) + me = material_homogenizationMemberAt2(ce) !-------------------------------------------------------------------------------------------------- ! compute the deformation gradient of individual grains due to relaxations F = 0.0_pReal @@ -219,8 +218,8 @@ module subroutine mechanical_RGC_partitionDeformation(F,avgF,instance,of) iGrain3 = grain1to3(iGrain,prm%N_constituents) do iFace = 1,6 intFace = getInterface(iFace,iGrain3) ! identifying 6 interfaces of each grain - aVect = relaxationVector(intFace,instance,of) ! get the relaxation vectors for each interface from global relaxation vector array - nVect = interfaceNormal(intFace,instance,of) + aVect = relaxationVector(intFace,ho,me) ! get the relaxation vectors for each interface from global relaxation vector array + nVect = interfaceNormal(intFace,ho,me) forall (i=1:3,j=1:3) & F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! calculating deformation relaxations due to interface relaxation enddo @@ -236,7 +235,7 @@ end subroutine mechanical_RGC_partitionDeformation !> @brief update the internal state of the homogenization scheme and tell whether "done" and ! "happy" with result !-------------------------------------------------------------------------------------------------- -module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAndHappy) +module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) logical, dimension(2) :: doneAndHappy real(pReal), dimension(:,:,:), intent(in) :: & P,& !< partitioned stresses @@ -245,12 +244,11 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn 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 + ce !< cell 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 + integer :: ho,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ipert,iGrain,nGrain, me real(pReal), dimension(3,3,size(P,3)) :: R,pF,pR,D,pD real(pReal), dimension(3,size(P,3)) :: NN,devNull real(pReal), dimension(3) :: normP,normN,mornP,mornN @@ -264,10 +262,10 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn return endif zeroTimeStep - instance = homogenization_typeInstance(material_homogenizationAt(el)) - of = material_homogenizationMemberAt(ip,el) + ho = material_homogenizationAt2(ce) - associate(stt => state(instance), st0 => state0(instance), dst => dependentState(instance), prm => param(instance)) + me = material_homogenizationMemberAt2(ce) + associate(stt => state(ho), st0 => state0(ho), dst => dependentState(ho), prm => param(ho)) !-------------------------------------------------------------------------------------------------- ! get the dimension of the cluster (grains and interfaces) @@ -281,36 +279,36 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn ! allocate the size of the global relaxation arrays/jacobian matrices depending on the size of the cluster allocate(resid(3*nIntFaceTot), source=0.0_pReal) allocate(tract(nIntFaceTot,3), source=0.0_pReal) - relax = stt%relaxationVector(:,of) - drelax = stt%relaxationVector(:,of) - st0%relaxationVector(:,of) + relax = stt%relaxationVector(:,me) + drelax = stt%relaxationVector(:,me) - st0%relaxationVector(:,me) !-------------------------------------------------------------------------------------------------- ! computing interface mismatch and stress penalty tensor for all interfaces of all grains - call stressPenalty(R,NN,avgF,F,ip,el,instance,of) + call stressPenalty(R,NN,avgF,F,ho,me) !-------------------------------------------------------------------------------------------------- ! calculating volume discrepancy and stress penalty related to overall volume discrepancy - call volumePenalty(D,dst%volumeDiscrepancy(of),avgF,F,nGrain,instance,of) + call volumePenalty(D,dst%volumeDiscrepancy(me),avgF,F,nGrain) !------------------------------------------------------------------------------------------------ ! computing the residual stress from the balance of traction at all (interior) interfaces do iNum = 1,nIntFaceTot - faceID = interface1to4(iNum,param(instance)%N_constituents) ! identifying the interface ID in local coordinate system (4-dimensional index) + faceID = interface1to4(iNum,param(ho)%N_constituents) ! identifying the interface ID in local coordinate system (4-dimensional index) !-------------------------------------------------------------------------------------------------- ! identify the left/bottom/back grain (-|N) iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate system (3-dimensional index) - iGrN = grain3to1(iGr3N,param(instance)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) + iGrN = grain3to1(iGr3N,param(ho)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) intFaceN = getInterface(2*faceID(1),iGr3N) - normN = interfaceNormal(intFaceN,instance,of) + normN = interfaceNormal(intFaceN,ho,me) !-------------------------------------------------------------------------------------------------- ! identify the right/up/front grain (+|P) iGr3P = iGr3N iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate system (3-dimensional index) - iGrP = grain3to1(iGr3P,param(instance)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) + iGrP = grain3to1(iGr3P,param(ho)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) intFaceP = getInterface(2*faceID(1)-1,iGr3P) - normP = interfaceNormal(intFaceP,instance,of) + normP = interfaceNormal(intFaceP,ho,me) !-------------------------------------------------------------------------------------------------- ! compute the residual of traction at the interface (in local system, 4-dimensional index) @@ -338,9 +336,9 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn if (residMax < num%rtol*stresMax .or. residMax < num%atol) then doneAndHappy = .true. - dst%mismatch(1:3,of) = sum(NN,2)/real(nGrain,pReal) - dst%relaxationRate_avg(of) = sum(abs(drelax))/dt/real(3*nIntFaceTot,pReal) - dst%relaxationRate_max(of) = maxval(abs(drelax))/dt + dst%mismatch(1:3,me) = sum(NN,2)/real(nGrain,pReal) + dst%relaxationRate_avg(me) = sum(abs(drelax))/dt/real(3*nIntFaceTot,pReal) + dst%relaxationRate_max(me) = maxval(abs(drelax))/dt return @@ -359,18 +357,18 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn ! ... of the constitutive stress tangent, assembled from dPdF or material constitutive model "smatrix" allocate(smatrix(3*nIntFaceTot,3*nIntFaceTot), source=0.0_pReal) do iNum = 1,nIntFaceTot - faceID = interface1to4(iNum,param(instance)%N_constituents) ! assembling of local dPdF into global Jacobian matrix + faceID = interface1to4(iNum,param(ho)%N_constituents) ! assembling of local dPdF into global Jacobian matrix !-------------------------------------------------------------------------------------------------- ! identify the left/bottom/back grain (-|N) iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate sytem - iGrN = grain3to1(iGr3N,param(instance)%N_constituents) ! translate into global grain ID + iGrN = grain3to1(iGr3N,param(ho)%N_constituents) ! translate into global grain ID intFaceN = getInterface(2*faceID(1),iGr3N) ! identifying the connecting interface in local coordinate system - normN = interfaceNormal(intFaceN,instance,of) + normN = interfaceNormal(intFaceN,ho,me) do iFace = 1,6 intFaceN = getInterface(iFace,iGr3N) ! identifying all interfaces that influence relaxation of the above interface - mornN = interfaceNormal(intFaceN,instance,of) - iMun = interface4to1(intFaceN,param(instance)%N_constituents) ! translate the interfaces ID into local 4-dimensional index + mornN = interfaceNormal(intFaceN,ho,me) + iMun = interface4to1(intFaceN,param(ho)%N_constituents) ! translate the interfaces ID into local 4-dimensional index if (iMun > 0) then ! get the corresponding tangent do i=1,3; do j=1,3; do k=1,3; do l=1,3 smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) & @@ -385,13 +383,13 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn ! identify the right/up/front grain (+|P) iGr3P = iGr3N iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate sytem - iGrP = grain3to1(iGr3P,param(instance)%N_constituents) ! translate into global grain ID + iGrP = grain3to1(iGr3P,param(ho)%N_constituents) ! translate into global grain ID intFaceP = getInterface(2*faceID(1)-1,iGr3P) ! identifying the connecting interface in local coordinate system - normP = interfaceNormal(intFaceP,instance,of) + normP = interfaceNormal(intFaceP,ho,me) do iFace = 1,6 intFaceP = getInterface(iFace,iGr3P) ! identifying all interfaces that influence relaxation of the above interface - mornP = interfaceNormal(intFaceP,instance,of) - iMun = interface4to1(intFaceP,param(instance)%N_constituents) ! translate the interfaces ID into local 4-dimensional index + mornP = interfaceNormal(intFaceP,ho,me) + iMun = interface4to1(intFaceP,param(ho)%N_constituents) ! translate the interfaces ID into local 4-dimensional index if (iMun > 0) then ! get the corresponding tangent do i=1,3; do j=1,3; do k=1,3; do l=1,3 smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) & @@ -411,31 +409,31 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn do ipert = 1,3*nIntFaceTot p_relax = relax p_relax(ipert) = relax(ipert) + num%pPert ! perturb the relaxation vector - stt%relaxationVector(:,of) = p_relax - call grainDeformation(pF,avgF,instance,of) ! rain deformation from perturbed state - call stressPenalty(pR,DevNull, avgF,pF,ip,el,instance,of) ! stress penalty due to interface mismatch from perturbed state - call volumePenalty(pD,devNull(1,1), avgF,pF,nGrain,instance,of) ! stress penalty due to volume discrepancy from perturbed state + stt%relaxationVector(:,me) = p_relax + call grainDeformation(pF,avgF,ho,me) ! rain deformation from perturbed state + call stressPenalty(pR,DevNull, avgF,pF,ho,me) ! stress penalty due to interface mismatch from perturbed state + call volumePenalty(pD,devNull(1,1), avgF,pF,nGrain) ! stress penalty due to volume discrepancy from perturbed state !-------------------------------------------------------------------------------------------------- ! computing the global stress residual array from the perturbed state p_resid = 0.0_pReal do iNum = 1,nIntFaceTot - faceID = interface1to4(iNum,param(instance)%N_constituents) ! identifying the interface ID in local coordinate system (4-dimensional index) + faceID = interface1to4(iNum,param(ho)%N_constituents) ! identifying the interface ID in local coordinate system (4-dimensional index) !-------------------------------------------------------------------------------------------------- ! identify the left/bottom/back grain (-|N) iGr3N = faceID(2:4) ! identify the grain ID in local coordinate system (3-dimensional index) - iGrN = grain3to1(iGr3N,param(instance)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) + iGrN = grain3to1(iGr3N,param(ho)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) intFaceN = getInterface(2*faceID(1),iGr3N) ! identify the interface ID of the grain - normN = interfaceNormal(intFaceN,instance,of) + normN = interfaceNormal(intFaceN,ho,me) !-------------------------------------------------------------------------------------------------- ! identify the right/up/front grain (+|P) iGr3P = iGr3N iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identify the grain ID in local coordinate system (3-dimensional index) - iGrP = grain3to1(iGr3P,param(instance)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) + iGrP = grain3to1(iGr3P,param(ho)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) intFaceP = getInterface(2*faceID(1)-1,iGr3P) ! identify the interface ID of the grain - normP = interfaceNormal(intFaceP,instance,of) + normP = interfaceNormal(intFaceP,ho,me) !-------------------------------------------------------------------------------------------------- ! compute the residual stress (contribution of mismatch and volume penalties) from perturbed state @@ -475,11 +473,11 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn do i = 1,3*nIntFaceTot;do j = 1,3*nIntFaceTot drelax(i) = drelax(i) - jnverse(i,j)*resid(j) ! Calculate the correction for the state variable enddo; enddo - stt%relaxationVector(:,of) = relax + drelax ! Updateing the state variable for the next iteration + stt%relaxationVector(:,me) = relax + drelax ! Updateing the state variable for the next iteration if (any(abs(drelax) > num%maxdRelax)) then ! Forcing cutback when the incremental change of relaxation vector becomes too large doneAndHappy = [.true.,.false.] !$OMP CRITICAL (write2out) - print'(a,i3,a,i3,a)',' RGC_updateState: ip ',ip,' | el ',el,' enforces cutback' + print'(a,i3,a,i3,a)',' RGC_updateState: enforces cutback' print'(a,e15.8)',' due to large relaxation change = ',maxval(abs(drelax)) flush(IO_STDOUT) !$OMP END CRITICAL (write2out) @@ -491,14 +489,14 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn !------------------------------------------------------------------------------------------------ !> @brief calculate stress-like penalty due to deformation mismatch !------------------------------------------------------------------------------------------------ - subroutine stressPenalty(rPen,nMis,avgF,fDef,ip,el,instance,of) + subroutine stressPenalty(rPen,nMis,avgF,fDef,ho,me) real(pReal), dimension (:,:,:), intent(out) :: rPen !< stress-like penalty real(pReal), dimension (:,:), intent(out) :: nMis !< total amount of mismatch real(pReal), dimension (:,:,:), intent(in) :: fDef !< deformation gradients real(pReal), dimension (3,3), intent(in) :: avgF !< initial effective stretch tensor - integer, intent(in) :: ip,el,instance,of + integer, intent(in) :: ho, me integer, dimension (4) :: intFace integer, dimension (3) :: iGrain3,iGNghb3,nGDim @@ -510,7 +508,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn nDefToler = 1.0e-10_pReal, & b = 2.5e-10_pReal ! Length of Burgers vector - nGDim = param(instance)%N_constituents + nGDim = param(ho)%N_constituents rPen = 0.0_pReal nMis = 0.0_pReal @@ -518,27 +516,26 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn ! get the correction factor the modulus of penalty stress representing the evolution of area of ! the interfaces due to deformations - surfCorr = surfaceCorrection(avgF,instance,of) - - associate(prm => param(instance)) + surfCorr = surfaceCorrection(avgF,ho,me) + associate(prm => param(ho)) !----------------------------------------------------------------------------------------------- ! computing the mismatch and penalty stress tensor of all grains grainLoop: do iGrain = 1,product(prm%N_constituents) - muGrain = equivalentMu(iGrain,ip,el) + muGrain = equivalentMu(iGrain,ce) iGrain3 = grain1to3(iGrain,prm%N_constituents) ! get the grain ID in local 3-dimensional index (x,y,z)-position interfaceLoop: do iFace = 1,6 intFace = getInterface(iFace,iGrain3) ! get the 4-dimensional index of the interface in local numbering system of the grain - nVect = interfaceNormal(intFace,instance,of) + nVect = interfaceNormal(intFace,ho,me) iGNghb3 = iGrain3 ! identify the neighboring grain across the interface iGNghb3(abs(intFace(1))) = iGNghb3(abs(intFace(1))) & + int(real(intFace(1),pReal)/real(abs(intFace(1)),pReal)) where(iGNghb3 < 1) iGNghb3 = nGDim where(iGNghb3 >nGDim) iGNghb3 = 1 iGNghb = grain3to1(iGNghb3,prm%N_constituents) ! get the ID of the neighboring grain - muGNghb = equivalentMu(iGNghb,ip,el) + muGNghb = equivalentMu(iGNghb,ce) gDef = 0.5_pReal*(fDef(1:3,1:3,iGNghb) - fDef(1:3,1:3,iGrain)) ! difference/jump in deformation gradeint across the neighbor !------------------------------------------------------------------------------------------- @@ -577,7 +574,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn !------------------------------------------------------------------------------------------------ !> @brief calculate stress-like penalty due to volume discrepancy !------------------------------------------------------------------------------------------------ - subroutine volumePenalty(vPen,vDiscrep,fAvg,fDef,nGrain,instance,of) + subroutine volumePenalty(vPen,vDiscrep,fAvg,fDef,nGrain) real(pReal), dimension (:,:,:), intent(out) :: vPen ! stress-like penalty due to volume real(pReal), intent(out) :: vDiscrep ! total volume discrepancy @@ -585,9 +582,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn real(pReal), dimension (:,:,:), intent(in) :: fDef ! deformation gradients real(pReal), dimension (3,3), intent(in) :: fAvg ! overall deformation gradient integer, intent(in) :: & - Ngrain, & - instance, & - of + Ngrain real(pReal), dimension(size(vPen,3)) :: gVol integer :: i @@ -617,14 +612,14 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn !> @brief compute the correction factor accouted for surface evolution (area change) due to ! deformation !-------------------------------------------------------------------------------------------------- - function surfaceCorrection(avgF,instance,of) + function surfaceCorrection(avgF,ho,me) real(pReal), dimension(3) :: surfaceCorrection real(pReal), dimension(3,3), intent(in) :: avgF !< average F integer, intent(in) :: & - instance, & - of + ho, & + me real(pReal), dimension(3,3) :: invC real(pReal), dimension(3) :: nVect real(pReal) :: detF @@ -635,7 +630,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn surfaceCorrection = 0.0_pReal do iBase = 1,3 - nVect = interfaceNormal([iBase,1,1,1],instance,of) + nVect = interfaceNormal([iBase,1,1,1],ho,me) do i = 1,3; do j = 1,3 surfaceCorrection(iBase) = surfaceCorrection(iBase) + invC(i,j)*nVect(i)*nVect(j) ! compute the component of (the inverse of) the stretch in the direction of the normal enddo; enddo @@ -648,17 +643,16 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn !------------------------------------------------------------------------------------------------- !> @brief compute the equivalent shear and bulk moduli from the elasticity tensor !------------------------------------------------------------------------------------------------- - real(pReal) function equivalentMu(grainID,ip,el) + real(pReal) function equivalentMu(grainID,ce) integer, intent(in) :: & grainID,& - ip, & !< integration point number - el !< element number + ce real(pReal), dimension(6,6) :: C - C = phase_homogenizedC(material_phaseAt(grainID,el),material_phaseMemberAt(grainID,ip,el)) + C = phase_homogenizedC(material_phaseAt2(grainID,ce),material_phaseMemberAt2(grainID,ce)) equivalentMu = lattice_equivalent_mu(C,'voigt') end function equivalentMu @@ -668,14 +662,14 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn !> @brief calculating the grain deformation gradient (the same with ! homogenization_RGC_partitionDeformation, but used only for perturbation scheme) !------------------------------------------------------------------------------------------------- - subroutine grainDeformation(F, avgF, instance, of) + subroutine grainDeformation(F, avgF, ho, me) real(pReal), dimension(:,:,:), intent(out) :: F !< partitioned F per grain real(pReal), dimension(:,:), intent(in) :: avgF !< averaged F integer, intent(in) :: & - instance, & - of + ho, & + me real(pReal), dimension(3) :: aVect,nVect integer, dimension(4) :: intFace @@ -685,15 +679,15 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn !----------------------------------------------------------------------------------------------- ! compute the deformation gradient of individual grains due to relaxations - associate(prm => param(instance)) + associate (prm => param(ho)) F = 0.0_pReal do iGrain = 1,product(prm%N_constituents) iGrain3 = grain1to3(iGrain,prm%N_constituents) do iFace = 1,6 intFace = getInterface(iFace,iGrain3) - aVect = relaxationVector(intFace,instance,of) - nVect = interfaceNormal(intFace,instance,of) + aVect = relaxationVector(intFace,ho,me) + nVect = interfaceNormal(intFace,ho,me) forall (i=1:3,j=1:3) & F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations enddo @@ -710,17 +704,17 @@ end function mechanical_RGC_updateState !-------------------------------------------------------------------------------------------------- !> @brief derive average stress and stiffness from constituent quantities !-------------------------------------------------------------------------------------------------- -module subroutine mechanical_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) +module subroutine mechanical_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho) 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 + integer, intent(in) :: ho - avgP = sum(P,3) /real(product(param(instance)%N_constituents),pReal) - dAvgPdAvgF = sum(dPdF,5)/real(product(param(instance)%N_constituents),pReal) + avgP = sum(P,3) /real(product(param(ho)%N_constituents),pReal) + dAvgPdAvgF = sum(dPdF,5)/real(product(param(ho)%N_constituents),pReal) end subroutine mechanical_RGC_averageStressAndItsTangent @@ -728,14 +722,14 @@ end subroutine mechanical_RGC_averageStressAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -module subroutine mechanical_RGC_results(instance,group) +module subroutine mechanical_RGC_results(ho,group) - integer, intent(in) :: instance + integer, intent(in) :: ho character(len=*), intent(in) :: group integer :: o - associate(stt => state(instance), dst => dependentState(instance), prm => param(instance)) + associate(stt => state(ho), dst => dependentState(ho), prm => param(ho)) outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) case('M') @@ -760,11 +754,11 @@ end subroutine mechanical_RGC_results !-------------------------------------------------------------------------------------------------- !> @brief collect relaxation vectors of an interface !-------------------------------------------------------------------------------------------------- -pure function relaxationVector(intFace,instance,of) +pure function relaxationVector(intFace,ho,me) real(pReal), dimension (3) :: relaxationVector - integer, intent(in) :: instance,of + integer, intent(in) :: ho,me integer, dimension(4), intent(in) :: intFace !< set of interface ID in 4D array (normal and position) integer :: iNum @@ -772,29 +766,35 @@ pure function relaxationVector(intFace,instance,of) !-------------------------------------------------------------------------------------------------- ! collect the interface relaxation vector from the global state array - iNum = interface4to1(intFace,param(instance)%N_constituents) ! identify the position of the interface in global state array + associate (prm => param(ho), & + stt => state(ho)) + + iNum = interface4to1(intFace,prm%N_constituents) ! identify the position of the interface in global state array if (iNum > 0) then - relaxationVector = state(instance)%relaxationVector((3*iNum-2):(3*iNum),of) + relaxationVector = stt%relaxationVector((3*iNum-2):(3*iNum),me) else relaxationVector = 0.0_pReal endif + end associate + end function relaxationVector !-------------------------------------------------------------------------------------------------- !> @brief identify the normal of an interface !-------------------------------------------------------------------------------------------------- -pure function interfaceNormal(intFace,instance,of) +pure function interfaceNormal(intFace,ho,me) real(pReal), dimension(3) :: interfaceNormal integer, dimension(4), intent(in) :: intFace !< interface ID in 4D array (normal and position) integer, intent(in) :: & - instance, & - of + ho, & + me integer :: nPos + associate (dst => dependentState(ho)) !-------------------------------------------------------------------------------------------------- ! get the normal of the interface, identified from the value of intFace(1) @@ -802,7 +802,9 @@ pure function interfaceNormal(intFace,instance,of) nPos = abs(intFace(1)) ! identify the position of the interface in global state array interfaceNormal(nPos) = real(intFace(1)/abs(intFace(1)),pReal) ! get the normal vector w.r.t. cluster axis - interfaceNormal = matmul(dependentState(instance)%orientation(1:3,1:3,of),interfaceNormal) ! map the normal vector into sample coordinate system (basis) + interfaceNormal = matmul(dst%orientation(1:3,1:3,me),interfaceNormal) ! map the normal vector into sample coordinate system (basis) + + end associate end function interfaceNormal diff --git a/src/homogenization_mechanical_isostrain.f90 b/src/homogenization_mechanical_isostrain.f90 index b492499d8..4773e3081 100644 --- a/src/homogenization_mechanical_isostrain.f90 +++ b/src/homogenization_mechanical_isostrain.f90 @@ -29,7 +29,6 @@ contains module subroutine mechanical_isostrain_init integer :: & - Ninstances, & h, & Nmaterialpoints class(tNode), pointer :: & @@ -39,17 +38,16 @@ module subroutine mechanical_isostrain_init print'(/,a)', ' <<<+- homogenization:mechanical:isostrain init -+>>>' - Ninstances = count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) - - allocate(param(Ninstances)) ! one container of parameters per instance + print'(a,i2)', ' # instances: ',count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID); flush(IO_STDOUT) material_homogenization => config_material%get('homogenization') + allocate(param(material_homogenization%length)) ! one container of parameters per homog + do h = 1, size(homogenization_type) if (homogenization_type(h) /= HOMOGENIZATION_ISOSTRAIN_ID) cycle homog => material_homogenization%get(h) homogMech => homog%get('mechanics') - associate(prm => param(homogenization_typeInstance(h))) + associate(prm => param(h)) prm%N_constituents = homogenization_Nconstituents(h) select case(homogMech%get_asString('mapping',defaultVal = 'sum')) @@ -90,16 +88,16 @@ end subroutine mechanical_isostrain_partitionDeformation !-------------------------------------------------------------------------------------------------- !> @brief derive average stress and stiffness from constituent quantities !-------------------------------------------------------------------------------------------------- -module subroutine mechanical_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) +module subroutine mechanical_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho) 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 + integer, intent(in) :: ho - associate(prm => param(instance)) + associate(prm => param(ho)) select case (prm%mapping) case (parallel_ID) diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90 index 67027f52a..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)) @@ -99,24 +99,21 @@ end subroutine thermal_homogenize !-------------------------------------------------------------------------------------------------- !> @brief return homogenized thermal conductivity in reference configuration !-------------------------------------------------------------------------------------------------- -module function thermal_conduction_getConductivity(ip,el) result(K) +module function thermal_conduction_getConductivity(ce) result(K) - integer, intent(in) :: & - ip, & !< integration point number - el !< element number + integer, intent(in) :: ce real(pReal), dimension(3,3) :: K integer :: & co - K = 0.0_pReal - do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - K = K + crystallite_push33ToRef(co,ip,el,lattice_K(:,:,material_phaseAt(co,el))) + do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) + K = K + crystallite_push33ToRef(co,ce,lattice_K(:,:,material_phaseAt2(co,ce))) enddo - K = K / real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) + K = K / real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) end function thermal_conduction_getConductivity @@ -219,11 +216,11 @@ end function homogenization_thermal_T !-------------------------------------------------------------------------------------------------- !> @brief return heat generation rate !-------------------------------------------------------------------------------------------------- -module subroutine thermal_conduction_getSource(Tdot, ip,el) +module subroutine thermal_conduction_getSource(Tdot, ip, el) integer, intent(in) :: & - ip, & !< integration point number - el !< element number + ip, & + el real(pReal), intent(out) :: & Tdot diff --git a/src/material.f90 b/src/material.f90 index 3656f0f0a..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 @@ -39,9 +39,6 @@ module material integer, dimension(:,:,:), allocatable, public, protected :: & ! (constituent,IP,elem) material_phaseMemberAt !< position of the element within its phase instance - type(Rotation), dimension(:,:,:), allocatable, public, protected :: & - material_orientation0 !< initial orientation of each grain,IP,element - public :: & material_init @@ -125,8 +122,6 @@ subroutine material_parseMaterial allocate(material_phaseAt2(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0) allocate(material_phaseMemberAt2(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0) - allocate(material_orientation0(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems)) - do el = 1, discretization_Nelems material => materials%get(discretization_materialAt(el)) constituents => material%get('constituents') @@ -135,7 +130,7 @@ subroutine material_parseMaterial do ip = 1, discretization_nIPs ce = (el-1)*discretization_nIPs + ip counterHomogenization(material_homogenizationAt(el)) = counterHomogenization(material_homogenizationAt(el)) + 1 - material_homogenizationMemberAt(ip,el) = counterHomogenization(material_homogenizationAt(el)) + material_homogenizationMemberAt(ip,el) = counterHomogenization(material_homogenizationAt(el)) material_homogenizationAt2(ce) = material_homogenizationAt(el) material_homogenizationMemberAt2(ce) = material_homogenizationMemberAt(ip,el) enddo @@ -153,7 +148,6 @@ subroutine material_parseMaterial material_phaseAt2(co,ce) = material_phaseAt(co,el) material_phaseMemberAt2(co,ce) = material_phaseMemberAt(co,ip,el) - call material_orientation0(co,ip,el)%fromQuaternion(constituent%get_asFloats('O',requiredSize=4)) ! should be done in crystallite enddo enddo diff --git a/src/math.f90 b/src/math.f90 index 6b89a9923..6ef942677 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -1133,6 +1133,7 @@ real(pReal) pure function math_areaTriangle(v1,v2,v3) real(pReal), dimension (3), intent(in) :: v1,v2,v3 + math_areaTriangle = 0.5_pReal * norm2(math_cross(v1-v2,v1-v3)) end function math_areaTriangle @@ -1147,11 +1148,13 @@ real(pReal) pure elemental function math_clip(a, left, right) real(pReal), intent(in) :: a real(pReal), intent(in), optional :: left, right + math_clip = a if (present(left)) math_clip = max(left,math_clip) if (present(right)) math_clip = min(right,math_clip) - if (present(left) .and. present(right)) & - math_clip = merge (IEEE_value(1.0_pReal,IEEE_quiet_NaN),math_clip, left>right) + if (present(left) .and. present(right)) then + if(left>right) error stop 'left > right' + endif end function math_clip @@ -1182,6 +1185,7 @@ subroutine selfTest integer :: d logical :: e + if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal,3.0_pReal,3.0_pReal,3.0_pReal] - & math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1,2,3,0])) > tol_math_check)) & error stop 'math_expand [1,2,3] by [1,2,3,0] => [1,2,2,3,3,3]' diff --git a/src/phase.f90 b/src/phase.f90 index a36ddc157..add26ae0b 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -19,6 +19,9 @@ module phase implicit none private + type(Rotation), dimension(:,:,:), allocatable :: & + material_orientation0 !< initial orientation of each grain,IP,element + type(rotation), dimension(:,:,:), allocatable :: & crystallite_orientation !< current orientation @@ -77,8 +80,8 @@ module phase interface ! == cleaned:begin ================================================================================= - module subroutine mechanical_init(phases) - class(tNode), pointer :: phases + module subroutine mechanical_init(materials,phases) + class(tNode), pointer :: materials,phases end subroutine mechanical_init module subroutine damage_init @@ -117,12 +120,11 @@ module phase end subroutine mechanical_restore - module function phase_mechanical_dPdF(dt,co,ip,el) result(dPdF) + module function phase_mechanical_dPdF(dt,co,ce) result(dPdF) real(pReal), intent(in) :: dt integer, intent(in) :: & co, & !< counter in constituent loop - ip, & !< counter in integration point loop - el !< counter in element loop + ce real(pReal), dimension(3,3,3,3) :: dPdF end function phase_mechanical_dPdF @@ -147,8 +149,8 @@ module phase real(pReal), dimension(3,3) :: L_p end function mechanical_L_p - module function phase_mechanical_getF(co,ip,el) result(F) - integer, intent(in) :: co, ip, el + module function phase_mechanical_getF(co,ce) result(F) + integer, intent(in) :: co, ce real(pReal), dimension(3,3) :: F end function phase_mechanical_getF @@ -157,8 +159,8 @@ module phase real(pReal), dimension(3,3) :: F_e end function mechanical_F_e - module function phase_mechanical_getP(co,ip,el) result(P) - integer, intent(in) :: co, ip, el + module function phase_mechanical_getP(co,ce) result(P) + integer, intent(in) :: co, ce real(pReal), dimension(3,3) :: P end function phase_mechanical_getP @@ -183,9 +185,9 @@ module phase end function damage_phi - module subroutine phase_mechanical_setF(F,co,ip,el) + module subroutine phase_mechanical_setF(F,co,ce) real(pReal), dimension(3,3), intent(in) :: F - integer, intent(in) :: co, ip, el + integer, intent(in) :: co, ce end subroutine phase_mechanical_setF module subroutine phase_thermal_setField(T,dot_T, co,ce) @@ -229,10 +231,8 @@ module phase end function phase_homogenizedC - module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) - integer, intent(in) :: & - ip, & !< integration point number - el !< element number + module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ce) + integer, intent(in) :: ce real(pReal), intent(in) :: & phi !< damage parameter real(pReal), intent(inout) :: & @@ -342,6 +342,7 @@ subroutine phase_init so !< counter in source loop class (tNode), pointer :: & debug_constitutive, & + materials, & phases @@ -356,9 +357,10 @@ subroutine phase_init debugConstitutive%grain = config_debug%get_asInt('grain',defaultVal = 1) + materials => config_material%get('material') phases => config_material%get('phase') - call mechanical_init(phases) + call mechanical_init(materials,phases) call damage_init call thermal_init(phases) @@ -624,19 +626,20 @@ end subroutine crystallite_orientations !-------------------------------------------------------------------------------------------------- !> @brief Map 2nd order tensor to reference config !-------------------------------------------------------------------------------------------------- -function crystallite_push33ToRef(co,ip,el, tensor33) +function crystallite_push33ToRef(co,ce, tensor33) real(pReal), dimension(3,3), intent(in) :: tensor33 integer, intent(in):: & - el, & - ip, & - co + co, & + ce real(pReal), dimension(3,3) :: crystallite_push33ToRef real(pReal), dimension(3,3) :: T + integer :: ph, me - - T = matmul(material_orientation0(co,ip,el)%asMatrix(),transpose(math_inv33(phase_mechanical_getF(co,ip,el)))) ! ToDo: initial orientation correct? + ph = material_phaseAt2(co,ce) + me = material_phaseMemberAt2(co,ce) + T = matmul(material_orientation0(co,ph,me)%asMatrix(),transpose(math_inv33(phase_mechanical_getF(co,ce)))) ! ToDo: initial orientation correct? crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T)) diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 5e4e43017..c85075288 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) @@ -179,11 +179,9 @@ end subroutine damage_init !---------------------------------------------------------------------------------------------- !< @brief returns local part of nonlocal damage driving force !---------------------------------------------------------------------------------------------- -module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) +module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ce) - integer, intent(in) :: & - ip, & !< integration point number - el !< element number + integer, intent(in) :: ce real(pReal), intent(in) :: & phi !< damage parameter real(pReal), intent(inout) :: & @@ -201,9 +199,9 @@ module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, phiDot = 0.0_pReal dPhiDot_dPhi = 0.0_pReal - do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - ph = material_phaseAt(co,el) - me = material_phasememberAt(co,ip,el) + do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) + ph = material_phaseAt2(co,ce) + me = material_phasememberAt2(co,ce) select case(phase_source(ph)) case (DAMAGE_ISOBRITTLE_ID) 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 57d4ea5ec..e642b22ef 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -177,21 +177,26 @@ contains !> @brief Initialize mechanical field related constitutive models !> @details Initialize elasticity, plasticity and stiffness degradation models. !-------------------------------------------------------------------------------------------------- -module subroutine mechanical_init(phases) +module subroutine mechanical_init(materials,phases) class(tNode), pointer :: & + materials, & phases integer :: & el, & ip, & co, & + ce, & ph, & me, & stiffDegradationCtr, & - Nconstituents + Nmembers class(tNode), pointer :: & num_crystallite, & + material, & + constituents, & + constituent, & phase, & mech, & elastic, & @@ -221,23 +226,25 @@ module subroutine mechanical_init(phases) allocate(phase_mechanical_P(phases%length)) allocate(phase_mechanical_S0(phases%length)) - do ph = 1, phases%length - Nconstituents = count(material_phaseAt == ph) * discretization_nIPs + allocate(material_orientation0(homogenization_maxNconstituents,phases%length,maxVal(material_phaseMemberAt))) - 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)) + do ph = 1, phases%length + Nmembers = count(material_phaseAt == ph) * discretization_nIPs + + 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') @@ -271,14 +278,19 @@ module subroutine mechanical_init(phases) enddo endif - !$OMP PARALLEL DO PRIVATE(ph,me) + 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)) + constituents => material%get('constituents') + constituent => constituents%get(co) ph = material_phaseAt(co,el) me = material_phaseMemberAt(co,ip,el) - phase_mechanical_Fp0(ph)%data(1:3,1:3,me) = material_orientation0(co,ip,el)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005) + call material_orientation0(co,ph,me)%fromQuaternion(constituent%get_asFloats('O',requiredSize=4)) + + phase_mechanical_Fp0(ph)%data(1:3,1:3,me) = material_orientation0(co,ph,me)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005) phase_mechanical_Fp0(ph)%data(1:3,1:3,me) = phase_mechanical_Fp0(ph)%data(1:3,1:3,me) & / math_det33(phase_mechanical_Fp0(ph)%data(1:3,1:3,me))**(1.0_pReal/3.0_pReal) phase_mechanical_Fi0(ph)%data(1:3,1:3,me) = math_I3 @@ -292,7 +304,6 @@ module subroutine mechanical_init(phases) enddo enddo; enddo - !$OMP END PARALLEL DO ! initialize plasticity @@ -342,12 +353,11 @@ end subroutine mechanical_init !> the elastic and intermediate deformation gradients using Hooke's law !-------------------------------------------------------------------------------------------------- subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & - Fe, Fi, co, ip, el) + Fe, Fi, ph, me) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element + ph, & + me real(pReal), intent(in), dimension(3,3) :: & Fe, & !< elastic deformation gradient Fi !< intermediate deformation gradient @@ -360,17 +370,15 @@ subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & real(pReal), dimension(3,3) :: E real(pReal), dimension(3,3,3,3) :: C integer :: & - ho, & !< homogenization - d, & !< counter in degradation loop - i, j, ph, me + d, & !< counter in degradation loop + i, j - ho = material_homogenizationAt(el) - C = math_66toSym3333(phase_homogenizedC(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el))) + C = math_66toSym3333(phase_homogenizedC(ph,me)) - DegradationLoop: do d = 1, phase_NstiffnessDegradations(material_phaseAt(co,el)) - degradationType: select case(phase_stiffnessDegradation(d,material_phaseAt(co,el))) + DegradationLoop: do d = 1, phase_NstiffnessDegradations(ph) + degradationType: select case(phase_stiffnessDegradation(d,ph)) case (STIFFNESS_DEGRADATION_damage_ID) degradationType - C = C * phase_damage_get_phi(co,ip,el)**2 + C = C * damage_phi(ph,me)**2 end select degradationType enddo DegradationLoop @@ -528,7 +536,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken) B = math_I3 - Delta_t*Lpguess Fe = matmul(matmul(A,B), invFi_new) call phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & - Fe, Fi_new, co, ip, el) + Fe, Fi_new, ph, me) call plastic_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, & S, Fi_new, ph,me) @@ -1250,13 +1258,12 @@ end subroutine mechanical_restore !-------------------------------------------------------------------------------------------------- !> @brief Calculate tangent (dPdF). !-------------------------------------------------------------------------------------------------- -module function phase_mechanical_dPdF(dt,co,ip,el) result(dPdF) +module function phase_mechanical_dPdF(dt,co,ce) result(dPdF) real(pReal), intent(in) :: dt integer, intent(in) :: & co, & !< counter in constituent loop - ip, & !< counter in integration point loop - el !< counter in element loop + ce real(pReal), dimension(3,3,3,3) :: dPdF integer :: & @@ -1281,12 +1288,12 @@ module function phase_mechanical_dPdF(dt,co,ip,el) result(dPdF) logical :: error - ph = material_phaseAt(co,el) - me = material_phaseMemberAt(co,ip,el) + ph = material_phaseAt2(co,ce) + me = material_phaseMemberAt2(co,ce) call phase_hooke_SandItsTangents(devNull,dSdFe,dSdFi, & phase_mechanical_Fe(ph)%data(1:3,1:3,me), & - phase_mechanical_Fi(ph)%data(1:3,1:3,me),co,ip,el) + phase_mechanical_Fi(ph)%data(1:3,1:3,me),ph,me) call phase_LiAndItsTangents(devNull,dLidS,dLidFi, & phase_mechanical_S(ph)%data(1:3,1:3,me), & phase_mechanical_Fi(ph)%data(1:3,1:3,me), & @@ -1311,7 +1318,7 @@ module function phase_mechanical_dPdF(dt,co,ip,el) result(dPdF) enddo; enddo call math_invert(temp_99,error,math_3333to99(lhs_3333)) if (error) then - call IO_warning(warning_ID=600,el=el,ip=ip,g=co, & + call IO_warning(warning_ID=600, & ext_msg='inversion error in analytic tangent calculation') dFidS = 0.0_pReal else @@ -1341,7 +1348,7 @@ module function phase_mechanical_dPdF(dt,co,ip,el) result(dPdF) call math_invert(temp_99,error,math_eye(9)+math_3333to99(lhs_3333)) if (error) then - call IO_warning(warning_ID=600,el=el,ip=ip,g=co, & + call IO_warning(warning_ID=600, & ext_msg='inversion error in analytic tangent calculation') dSdF = rhs_3333 else @@ -1440,13 +1447,13 @@ end function mechanical_L_p !---------------------------------------------------------------------------------------------- !< @brief Get deformation gradient (for use by homogenization) !---------------------------------------------------------------------------------------------- -module function phase_mechanical_getF(co,ip,el) result(F) +module function phase_mechanical_getF(co,ce) result(F) - integer, intent(in) :: co, ip, el + integer, intent(in) :: co, ce real(pReal), dimension(3,3) :: F - F = phase_mechanical_F(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)) + F = phase_mechanical_F(material_phaseAt2(co,ce))%data(1:3,1:3,material_phaseMemberAt2(co,ce)) end function phase_mechanical_getF @@ -1469,25 +1476,25 @@ end function mechanical_F_e !---------------------------------------------------------------------------------------------- !< @brief Get second Piola-Kichhoff stress (for use by homogenization) !---------------------------------------------------------------------------------------------- -module function phase_mechanical_getP(co,ip,el) result(P) +module function phase_mechanical_getP(co,ce) result(P) - integer, intent(in) :: co, ip, el + integer, intent(in) :: co, ce real(pReal), dimension(3,3) :: P - P = phase_mechanical_P(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)) + P = phase_mechanical_P(material_phaseAt2(co,ce))%data(1:3,1:3,material_phaseMemberAt2(co,ce)) end function phase_mechanical_getP ! setter for homogenization -module subroutine phase_mechanical_setF(F,co,ip,el) +module subroutine phase_mechanical_setF(F,co,ce) real(pReal), dimension(3,3), intent(in) :: F - integer, intent(in) :: co, ip, el + integer, intent(in) :: co, ce - phase_mechanical_F(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)) = F + phase_mechanical_F(material_phaseAt2(co,ce))%data(1:3,1:3,material_phaseMemberAt2(co,ce)) = F end subroutine phase_mechanical_setF 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 0c52176ce..5500ae731 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -24,7 +24,6 @@ submodule(phase:plastic) dislotwin q_sb = 1.0_pReal, & !< q-exponent in shear band velocity D_a = 1.0_pReal, & !< adjustment parameter to calculate minimum dipole distance i_tw = 1.0_pReal, & !< adjustment parameter to calculate MFP for twinning - tau_0 = 1.0_pReal, & !< strength due to elements in solid solution L_tw = 1.0_pReal, & !< Length of twin nuclei in Burgers vectors L_tr = 1.0_pReal, & !< Length of trans nuclei in Burgers vectors x_c_tw = 1.0_pReal, & !< critical distance for formation of twin nucleus @@ -53,6 +52,7 @@ submodule(phase:plastic) dislotwin q, & !< q-exponent in glide velocity r, & !< r-exponent in twin nucleation rate s, & !< s-exponent in trans nucleation rate + tau_0, & !< strength due to elements in solid solution gamma_char, & !< characteristic shear for twins B !< drag coefficient real(pReal), allocatable, dimension(:,:) :: & @@ -81,7 +81,7 @@ submodule(phase:plastic) dislotwin logical :: & ExtendedDislocations, & !< consider split into partials for climb calculation fccTwinTransNucleation, & !< twinning and transformation models are for fcc - dipoleFormation !< flag indicating consideration of dipole formation + omitDipoles !< flag controlling consideration of dipole formation end type !< container type for internal constitutive parameters type :: tDislotwinState @@ -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 :: & @@ -213,10 +213,10 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%i_sl = pl%get_asFloats('i_sl', requiredSize=size(N_sl)) prm%p = pl%get_asFloats('p_sl', requiredSize=size(N_sl)) prm%q = pl%get_asFloats('q_sl', requiredSize=size(N_sl)) + prm%tau_0 = pl%get_asFloats('tau_0', requiredSize=size(N_sl)) prm%B = pl%get_asFloats('B', requiredSize=size(N_sl), & defaultVal=[(0.0_pReal, i=1,size(N_sl))]) - prm%tau_0 = pl%get_asFloat('tau_0') prm%D_a = pl%get_asFloat('D_a') prm%D_0 = pl%get_asFloat('D_0') prm%Q_cl = pl%get_asFloat('Q_cl') @@ -226,7 +226,7 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%dGamma_sf_dT = pl%get_asFloat('dGamma_sf_dT') endif - prm%dipoleformation = .not. pl%get_asBool('no_dipole_formation',defaultVal = .false.) + prm%omitDipoles = pl%get_asBool('omit_dipoles',defaultVal = .false.) ! multiplication factor according to crystal structure (nearest neighbors bcc vs fcc/hex) ! details: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 @@ -242,6 +242,7 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%i_sl = math_expand(prm%i_sl, N_sl) prm%p = math_expand(prm%p, N_sl) prm%q = math_expand(prm%q, N_sl) + prm%tau_0 = math_expand(prm%tau_0, N_sl) prm%B = math_expand(prm%B, N_sl) ! sanity checks @@ -405,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' @@ -427,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) @@ -443,28 +444,28 @@ module function plastic_dislotwin_init() result(myPlasticity) endIndex = endIndex + prm%sum_N_tw stt%f_tw=>plasticState(ph)%state(startIndex:endIndex,:) dot%f_tw=>plasticState(ph)%dotState(startIndex:endIndex,:) - plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('f_twin',defaultVal=1.0e-7_pReal) - if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' f_twin' + plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_f_tw',defaultVal=1.0e-7_pReal) + if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_f_tw' startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tr stt%f_tr=>plasticState(ph)%state(startIndex:endIndex,:) dot%f_tr=>plasticState(ph)%dotState(startIndex:endIndex,:) - plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('f_trans',defaultVal=1.0e-6_pReal) - if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' f_trans' + 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 @@ -535,9 +536,9 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me) real(pReal), dimension(param(ph)%sum_N_sl) :: & dot_gamma_sl,ddot_gamma_dtau_slip real(pReal), dimension(param(ph)%sum_N_tw) :: & - dot_gamma_twin,ddot_gamma_dtau_twin + dot_gamma_tw,ddot_gamma_dtau_tw real(pReal), dimension(param(ph)%sum_N_tr) :: & - dot_gamma_tr,ddot_gamma_dtau_trans + dot_gamma_tr,ddot_gamma_dtau_tr real(pReal):: dot_gamma_sb real(pReal), dimension(3,3) :: eigVectors, P_sb real(pReal), dimension(3) :: eigValues @@ -578,20 +579,20 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me) + ddot_gamma_dtau_slip(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i) enddo slipContribution - call kinetics_twin(Mp,T,dot_gamma_sl,ph,me,dot_gamma_twin,ddot_gamma_dtau_twin) + call kinetics_twin(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tw,ddot_gamma_dtau_tw) twinContibution: do i = 1, prm%sum_N_tw - Lp = Lp + dot_gamma_twin(i)*prm%P_tw(1:3,1:3,i) + Lp = Lp + dot_gamma_tw(i)*prm%P_tw(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + ddot_gamma_dtau_twin(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) + + ddot_gamma_dtau_tw(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) enddo twinContibution - call kinetics_trans(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tr,ddot_gamma_dtau_trans) + call kinetics_trans(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tr,ddot_gamma_dtau_tr) transContibution: do i = 1, prm%sum_N_tr Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + ddot_gamma_dtau_trans(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i) + + ddot_gamma_dtau_tr(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i) enddo transContibution Lp = Lp * f_unrotated @@ -646,7 +647,6 @@ module subroutine dislotwin_dotState(Mp,T,ph,me) f_unrotated, & rho_dip_distance, & v_cl, & !< climb velocity - Gamma, & !< stacking fault energy tau, & sigma_cl, & !< climb stress b_d !< ratio of Burgers vector to stacking fault width @@ -656,7 +656,7 @@ module subroutine dislotwin_dotState(Mp,T,ph,me) rho_dip_distance_min, & dot_gamma_sl real(pReal), dimension(param(ph)%sum_N_tw) :: & - dot_gamma_twin + dot_gamma_tw real(pReal), dimension(param(ph)%sum_N_tr) :: & dot_gamma_tr @@ -675,7 +675,7 @@ module subroutine dislotwin_dotState(Mp,T,ph,me) slipState: do i = 1, prm%sum_N_sl tau = math_tensordot(Mp,prm%P_sl(1:3,1:3,i)) - significantSlipStress: if (dEq0(tau)) then + significantSlipStress: if (dEq0(tau) .or. prm%omitDipoles) then dot_rho_dip_formation(i) = 0.0_pReal dot_rho_dip_climb(i) = 0.0_pReal else significantSlipStress @@ -683,24 +683,18 @@ module subroutine dislotwin_dotState(Mp,T,ph,me) rho_dip_distance = math_clip(rho_dip_distance, right = dst%Lambda_sl(i,me)) rho_dip_distance = math_clip(rho_dip_distance, left = rho_dip_distance_min(i)) - if (prm%dipoleFormation) then - dot_rho_dip_formation(i) = 2.0_pReal*(rho_dip_distance-rho_dip_distance_min(i))/prm%b_sl(i) & - * stt%rho_mob(i,me)*abs(dot_gamma_sl(i)) - else - dot_rho_dip_formation(i) = 0.0_pReal - endif + dot_rho_dip_formation(i) = 2.0_pReal*(rho_dip_distance-rho_dip_distance_min(i))/prm%b_sl(i) & + * stt%rho_mob(i,me)*abs(dot_gamma_sl(i)) if (dEq(rho_dip_distance,rho_dip_distance_min(i))) then dot_rho_dip_climb(i) = 0.0_pReal else - !@details: Refer: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 + ! Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 sigma_cl = dot_product(prm%n0_sl(1:3,i),matmul(Mp,prm%n0_sl(1:3,i))) - if (prm%ExtendedDislocations) then - Gamma = prm%Gamma_sf_0K + prm%dGamma_sf_dT * T - b_d = 24.0_pReal*PI*(1.0_pReal - prm%nu)/(2.0_pReal + prm%nu)* Gamma/(prm%mu*prm%b_sl(i)) - else - b_d = 1.0_pReal - endif + b_d = merge(24.0_pReal*PI*(1.0_pReal - prm%nu)/(2.0_pReal + prm%nu) & + * (prm%Gamma_sf_0K + prm%dGamma_sf_dT * T) / (prm%mu*prm%b_sl(i)), & + 1.0_pReal, & + prm%ExtendedDislocations) v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Q_cl/(kB*T)) & * (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal) @@ -718,8 +712,8 @@ module subroutine dislotwin_dotState(Mp,T,ph,me) - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,me)*abs(dot_gamma_sl) & - dot_rho_dip_climb - call kinetics_twin(Mp,T,dot_gamma_sl,ph,me,dot_gamma_twin) - dot%f_tw(:,me) = f_unrotated*dot_gamma_twin/prm%gamma_char + call kinetics_twin(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tw) + dot%f_tw(:,me) = f_unrotated*dot_gamma_tw/prm%gamma_char call kinetics_trans(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tr) dot%f_tr(:,me) = f_unrotated*dot_gamma_tr @@ -741,11 +735,9 @@ module subroutine dislotwin_dependentState(T,ph,me) T real(pReal) :: & - sumf_twin,Gamma,sumf_trans + sumf_tw,Gamma,sumf_tr real(pReal), dimension(param(ph)%sum_N_sl) :: & - inv_lambda_sl_sl, & !< 1/mean free distance between 2 forest dislocations seen by a moving dislocation - inv_lambda_sl_tw, & !< 1/mean free distance between 2 twin stacks from different systems seen by a moving dislocation - inv_lambda_sl_tr !< 1/mean free distance between 2 martensite lamellar from different systems seen by a moving dislocation + inv_lambda_sl real(pReal), dimension(param(ph)%sum_N_tw) :: & inv_lambda_tw_tw, & !< 1/mean free distance between 2 twin stacks from different systems seen by a growing twin f_over_t_tw @@ -760,38 +752,27 @@ module subroutine dislotwin_dependentState(T,ph,me) stt => state(ph),& dst => dependentState(ph)) - sumf_twin = sum(stt%f_tw(1:prm%sum_N_tw,me)) - sumf_trans = sum(stt%f_tr(1:prm%sum_N_tr,me)) + sumf_tw = sum(stt%f_tw(1:prm%sum_N_tw,me)) + sumf_tr = sum(stt%f_tr(1:prm%sum_N_tr,me)) Gamma = prm%Gamma_sf_0K + prm%dGamma_sf_dT * T !* rescaled volume fraction for topology f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,me)/prm%t_tw ! this is per system ... - f_over_t_tr = sumf_trans/prm%t_tr ! but this not + f_over_t_tr = sumf_tr/prm%t_tr ! but this not ! ToDo ...Physically correct, but naming could be adjusted - inv_lambda_sl_sl = sqrt(matmul(prm%forestProjection, & - stt%rho_mob(:,me)+stt%rho_dip(:,me)))/prm%i_sl - + inv_lambda_sl = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,me)+stt%rho_dip(:,me)))/prm%i_sl if (prm%sum_N_tw > 0 .and. prm%sum_N_sl > 0) & - inv_lambda_sl_tw = matmul(prm%h_sl_tw,f_over_t_tw)/(1.0_pReal-sumf_twin) - - inv_lambda_tw_tw = matmul(prm%h_tw_tw,f_over_t_tw)/(1.0_pReal-sumf_twin) - + inv_lambda_sl = inv_lambda_sl + matmul(prm%h_sl_tw,f_over_t_tw)/(1.0_pReal-sumf_tw) if (prm%sum_N_tr > 0 .and. prm%sum_N_sl > 0) & - inv_lambda_sl_tr = matmul(prm%h_sl_tr,f_over_t_tr)/(1.0_pReal-sumf_trans) - - inv_lambda_tr_tr = matmul(prm%h_tr_tr,f_over_t_tr)/(1.0_pReal-sumf_trans) - - if ((prm%sum_N_tw > 0) .or. (prm%sum_N_tr > 0)) then ! ToDo: better logic needed here - dst%Lambda_sl(:,me) = prm%D & - / (1.0_pReal+prm%D*(inv_lambda_sl_sl + inv_lambda_sl_tw + inv_lambda_sl_tr)) - else - dst%Lambda_sl(:,me) = prm%D & - / (1.0_pReal+prm%D*inv_lambda_sl_sl) !!!!!! correct? - endif + inv_lambda_sl = inv_lambda_sl + matmul(prm%h_sl_tr,f_over_t_tr)/(1.0_pReal-sumf_tr) + dst%Lambda_sl(:,me) = prm%D / (1.0_pReal+prm%D*inv_lambda_sl) + inv_lambda_tw_tw = matmul(prm%h_tw_tw,f_over_t_tw)/(1.0_pReal-sumf_tw) dst%Lambda_tw(:,me) = prm%i_tw*prm%D/(1.0_pReal+prm%D*inv_lambda_tw_tw) + + inv_lambda_tr_tr = matmul(prm%h_tr_tr,f_over_t_tr)/(1.0_pReal-sumf_tr) dst%Lambda_tr(:,me) = prm%i_tr*prm%D/(1.0_pReal+prm%D*inv_lambda_tr_tr) !* threshold stress for dislocation motion @@ -957,7 +938,7 @@ end subroutine kinetics_slip ! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,me,& - dot_gamma_twin,ddot_gamma_dtau_twin) + dot_gamma_tw,ddot_gamma_dtau_tw) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -970,9 +951,9 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,me,& dot_gamma_sl real(pReal), dimension(param(ph)%sum_N_tw), intent(out) :: & - dot_gamma_twin + dot_gamma_tw real(pReal), dimension(param(ph)%sum_N_tw), optional, intent(out) :: & - ddot_gamma_dtau_twin + ddot_gamma_dtau_tw real, dimension(param(ph)%sum_N_tw) :: & tau, & @@ -1004,16 +985,16 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,me,& significantStress: where(tau > tol_math_check) StressRatio_r = (dst%tau_hat_tw(:,me)/tau)**prm%r - dot_gamma_twin = prm%gamma_char * dst%V_tw(:,me) * Ndot0*exp(-StressRatio_r) - ddot_gamma_dtau = (dot_gamma_twin*prm%r/tau)*StressRatio_r + dot_gamma_tw = prm%gamma_char * dst%V_tw(:,me) * Ndot0*exp(-StressRatio_r) + ddot_gamma_dtau = (dot_gamma_tw*prm%r/tau)*StressRatio_r else where significantStress - dot_gamma_twin = 0.0_pReal + dot_gamma_tw = 0.0_pReal ddot_gamma_dtau = 0.0_pReal end where significantStress end associate - if(present(ddot_gamma_dtau_twin)) ddot_gamma_dtau_twin = ddot_gamma_dtau + if(present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw = ddot_gamma_dtau end subroutine kinetics_twin @@ -1026,7 +1007,7 @@ end subroutine kinetics_twin ! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,me,& - dot_gamma_tr,ddot_gamma_dtau_trans) + dot_gamma_tr,ddot_gamma_dtau_tr) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -1041,7 +1022,7 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,me,& real(pReal), dimension(param(ph)%sum_N_tr), intent(out) :: & dot_gamma_tr real(pReal), dimension(param(ph)%sum_N_tr), optional, intent(out) :: & - ddot_gamma_dtau_trans + ddot_gamma_dtau_tr real, dimension(param(ph)%sum_N_tr) :: & tau, & @@ -1081,7 +1062,7 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,me,& end associate - if(present(ddot_gamma_dtau_trans)) ddot_gamma_dtau_trans = ddot_gamma_dtau + if(present(ddot_gamma_dtau_tr)) ddot_gamma_dtau_tr = ddot_gamma_dtau end subroutine kinetics_trans 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 fbfcfa5af..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 @@ -1403,8 +1403,10 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) integer :: & n, & ! neighbor index + me, & neighbor_e, & ! element index of my neighbor neighbor_i, & ! integration point index of my neighbor + neighbor_me, & neighbor_phase, & ns, & ! number of active slip systems s1, & ! slip system index (me) @@ -1422,6 +1424,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) associate(prm => param(ph)) ns = prm%sum_N_sl + me = material_phaseMemberAt(1,i,e) !*** start out fully compatible my_compatibility = 0.0_pReal forall(s1 = 1:ns) my_compatibility(:,s1,s1,:) = 1.0_pReal @@ -1429,7 +1432,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) neighbors: do n = 1,nIPneighbors neighbor_e = IPneighborhood(1,n,i,e) neighbor_i = IPneighborhood(2,n,i,e) - + neighbor_me = material_phaseMemberAt(1,neighbor_i,neighbor_e) neighbor_phase = material_phaseAt(1,neighbor_e) if (neighbor_e <= 0 .or. neighbor_i <= 0) then @@ -1447,8 +1450,8 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) elseif (prm%chi_GB >= 0.0_pReal) then !* GRAIN BOUNDARY ! !* fixed transmissivity for adjacent ips with different texture (only if explicitly given in material.config) - if (any(dNeq(material_orientation0(1,i,e)%asQuaternion(), & - material_orientation0(1,neighbor_i,neighbor_e)%asQuaternion())) .and. & + if (any(dNeq(material_orientation0(1,ph,me)%asQuaternion(), & + material_orientation0(1,neighbor_phase,neighbor_me)%asQuaternion())) .and. & (.not. phase_localPlasticity(neighbor_phase))) & forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB) else @@ -1576,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, & @@ -1599,7 +1602,7 @@ subroutine stateInit(ini,phase,Nconstituents) totalVolume, & densityBinning, & minimumIpVolume - real(pReal), dimension(Nconstituents) :: & + real(pReal), dimension(Nmembers) :: & volume @@ -1619,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)) @@ -1806,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 @@ -1819,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 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