From c5d0c7e52ea63eb1e0d15b04897cd15171b6b2af Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Jun 2020 10:43:07 +0200 Subject: [PATCH 01/43] easier to read, more flexible --- python/damask/__init__.py | 4 ++-- python/damask/_result.py | 10 ++++------ python/damask/_vtk.py | 12 ++++++------ 3 files changed, 12 insertions(+), 14 deletions(-) diff --git a/python/damask/__init__.py b/python/damask/__init__.py index 0f343a581..6d65f7cd1 100644 --- a/python/damask/__init__.py +++ b/python/damask/__init__.py @@ -1,9 +1,9 @@ """Tools for pre and post processing of DAMASK simulations.""" -import os as _os +from pathlib import Path as _Path import re as _re name = 'damask' -with open(_os.path.join(_os.path.dirname(__file__),'VERSION')) as _f: +with open(_Path(__file__).parent/_Path('VERSION')) as _f: version = _re.sub(r'^v','',_f.readline().strip()) # make classes directly accessible as damask.Class diff --git a/python/damask/_result.py b/python/damask/_result.py index 73a038c4c..b6ac3782e 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -6,6 +6,7 @@ import os import datetime import xml.etree.ElementTree as ET import xml.dom.minidom +from pathlib import Path from functools import partial import h5py @@ -88,7 +89,7 @@ class Result: 'con_physics': self.con_physics, 'mat_physics': self.mat_physics } - self.fname = os.path.abspath(fname) + self.fname = Path(fname).absolute() self._allow_overwrite = False @@ -1163,7 +1164,7 @@ class Result: 'Dimensions': '{} {} {} {}'.format(*self.grid,np.prod(shape))} data_items[-1].text='{}:{}'.format(os.path.split(self.fname)[1],name) - with open(os.path.splitext(self.fname)[0]+'.xdmf','w') as f: + with open(self.fname.with_suffix('.xdmf').name,'w') as f: f.write(xml.dom.minidom.parseString(ET.tostring(xdmf).decode()).toprettyxml()) @@ -1239,7 +1240,4 @@ class Result: u = self.read_dataset(self.get_dataset_location('u_n' if mode.lower() == 'cell' else 'u_p')) v.add(u,'u') - file_out = '{}_inc{}'.format(os.path.splitext(os.path.basename(self.fname))[0], - inc[3:].zfill(N_digits)) - - v.write(file_out) + v.write('{}_inc{}'.format(self.fname.stem,inc[3:].zfill(N_digits))) diff --git a/python/damask/_vtk.py b/python/damask/_vtk.py index 58b43dc1f..72e224bda 100644 --- a/python/damask/_vtk.py +++ b/python/damask/_vtk.py @@ -1,4 +1,4 @@ -import os +from pathlib import Path import pandas as pd import numpy as np @@ -126,8 +126,8 @@ class VTK: vtkUnstructuredGrid, and vtkPolyData. """ - ext = os.path.splitext(fname)[1] - if ext == '.vtk': + ext = Path(fname).suffix + if ext == '.vtk' or dataset_type: reader = vtk.vtkGenericDataObjectReader() reader.SetFileName(fname) reader.Update() @@ -176,10 +176,10 @@ class VTK: writer = vtk.vtkXMLPolyDataWriter() default_ext = writer.GetDefaultFileExtension() - name, ext = os.path.splitext(fname) + ext = Path(fname).suffix if ext and ext != '.'+default_ext: - raise ValueError('Given extension {} is not .{}'.format(ext,default_ext)) - writer.SetFileName('{}.{}'.format(name,default_ext)) + raise ValueError('Given extension {} does not match default .{}'.format(ext,default_ext)) + writer.SetFileName(str(Path(fname).with_suffix('.'+default_ext))) writer.SetCompressorTypeToZLib() writer.SetDataModeToBinary() writer.SetInputData(self.geom) From c67fbacfc7fb828bf9ac8e462c3b0e3e87510b71 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Jun 2020 11:03:31 +0200 Subject: [PATCH 02/43] higher test coverage - invalid operations - legacy output --- python/tests/test_VTK.py | 34 +++++++++++++++++++++++++--------- 1 file changed, 25 insertions(+), 9 deletions(-) diff --git a/python/tests/test_VTK.py b/python/tests/test_VTK.py index 627ea4007..8795e7161 100644 --- a/python/tests/test_VTK.py +++ b/python/tests/test_VTK.py @@ -17,18 +17,24 @@ class TestVTK: size = np.random.random(3) + 1.0 origin = np.random.random(3) v = VTK.from_rectilinearGrid(grid,size,origin) - s = v.__repr__() + string = v.__repr__() v.write(os.path.join(tmp_path,'rectilinearGrid')) - v = VTK.from_file(os.path.join(tmp_path,'rectilinearGrid.vtr')) - assert(s == v.__repr__()) + vtr = VTK.from_file(os.path.join(tmp_path,'rectilinearGrid.vtr')) + with open(os.path.join(tmp_path,'rectilinearGrid.vtk'),'w') as f: + f.write(string) + vtk = VTK.from_file(os.path.join(tmp_path,'rectilinearGrid.vtk'),'VTK_rectilinearGrid') + assert(string == vtr.__repr__() == vtk.__repr__()) def test_polyData(self,tmp_path): points = np.random.rand(3,100) v = VTK.from_polyData(points) - s = v.__repr__() + string = v.__repr__() v.write(os.path.join(tmp_path,'polyData')) - v = VTK.from_file(os.path.join(tmp_path,'polyData.vtp')) - assert(s == v.__repr__()) + vtp = VTK.from_file(os.path.join(tmp_path,'polyData.vtp')) + with open(os.path.join(tmp_path,'polyData.vtk'),'w') as f: + f.write(string) + vtk = VTK.from_file(os.path.join(tmp_path,'polyData.vtk'),'polyData') + assert(string == vtp.__repr__() == vtk.__repr__()) @pytest.mark.parametrize('cell_type,n',[ ('VTK_hexahedron',8), @@ -41,7 +47,17 @@ class TestVTK: nodes = np.random.rand(n,3) connectivity = np.random.choice(np.arange(n),n,False).reshape(-1,n) v = VTK.from_unstructuredGrid(nodes,connectivity,cell_type) - s = v.__repr__() + string = v.__repr__() v.write(os.path.join(tmp_path,'unstructuredGrid')) - v = VTK.from_file(os.path.join(tmp_path,'unstructuredGrid.vtu')) - assert(s == v.__repr__()) + vtu = VTK.from_file(os.path.join(tmp_path,'unstructuredGrid.vtu')) + with open(os.path.join(tmp_path,'unstructuredGrid.vtk'),'w') as f: + f.write(string) + vtk = VTK.from_file(os.path.join(tmp_path,'unstructuredGrid.vtk'),'unstructuredgrid') + assert(string == vtu.__repr__() == vtk.__repr__()) + + @pytest.mark.parametrize('name,dataset_type',[('this_file_does_not_exist.vtk',None), + ('this_file_does_not_exist.vtk','vtk'), + ('this_file_does_not_exist.vtx', None)]) + def test_invalid_dataset_type(self,dataset_type,name): + with pytest.raises(TypeError): + VTK.from_file('this_file_does_not_exist.vtk',dataset_type) From 3be1a33820c5ab276c4e487b0276b689732849c0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Jun 2020 13:32:47 +0200 Subject: [PATCH 03/43] easier to read --- installation/symlink_Processing.py | 56 +++++++-------------------- processing/pre/mentat_pbcOnBoxMesh.py | 2 +- processing/pre/mentat_spectralBox.py | 2 +- python/damask/_environment.py | 25 ++++++------ python/damask/solver/_marc.py | 27 ++++++------- 5 files changed, 44 insertions(+), 68 deletions(-) diff --git a/installation/symlink_Processing.py b/installation/symlink_Processing.py index 5972a7f6a..9ee6db99d 100755 --- a/installation/symlink_Processing.py +++ b/installation/symlink_Processing.py @@ -1,56 +1,30 @@ #!/usr/bin/env python3 # Makes postprocessing routines accessible from everywhere. -import os import sys +from pathlib import Path import damask -damaskEnv = damask.Environment() -baseDir = damaskEnv.relPath('processing/') -binDir = damaskEnv.relPath('bin/') +env = damask.Environment() +bin_dir = env.root_dir/Path('bin') -if not os.path.isdir(binDir): - os.mkdir(binDir) +if not bin_dir.exists(): + bin_dir.mkdir() -#define ToDo list -processing_subDirs = ['pre', - 'post', - ] sys.stdout.write('\nsymbolic linking...\n') +for sub_dir in ['pre','post']: + the_dir = env.root_dir/Path('processing')/Path(sub_dir) -for subDir in processing_subDirs: - theDir = os.path.abspath(os.path.join(baseDir,subDir)) - - sys.stdout.write('\n'+binDir+' ->\n'+theDir+damask.util.deemph(' ...')+'\n') - - for theFile in os.listdir(theDir): - theName,theExt = os.path.splitext(theFile) - if theExt in ['.py']: - - src = os.path.abspath(os.path.join(theDir,theFile)) - sym_link = os.path.abspath(os.path.join(binDir,theName)) - - if os.path.lexists(sym_link): - os.remove(sym_link) - output = theName+damask.util.deemph(theExt) - else: - output = damask.util.emph(theName)+damask.util.deemph(theExt) - - sys.stdout.write(damask.util.deemph('... ')+output+'\n') - os.symlink(src,sym_link) + for the_file in the_dir.glob('*.py'): + src = the_dir/the_file + dst = bin_dir/Path(the_file.with_suffix('').name) + dst.unlink(True) + dst.symlink_to(src) sys.stdout.write('\npruning broken links...\n') - -brokenLinks = 0 - -for filename in os.listdir(binDir): - path = os.path.join(binDir,filename) - if os.path.islink(path) and not os.path.exists(path): - sys.stdout.write(' '+damask.util.delete(path)+'\n') - os.remove(path) - brokenLinks += 1 - -sys.stdout.write(('none.' if brokenLinks == 0 else '')+'\n') +for filename in bin_dir.glob('*'): + if not filename.is_file(): + filename.unlink diff --git a/processing/pre/mentat_pbcOnBoxMesh.py b/processing/pre/mentat_pbcOnBoxMesh.py index c171c6ccd..defcbb09b 100755 --- a/processing/pre/mentat_pbcOnBoxMesh.py +++ b/processing/pre/mentat_pbcOnBoxMesh.py @@ -6,7 +6,7 @@ import numpy as np from optparse import OptionParser import damask -sys.path.append(damask.solver.Marc().libraryPath()) +sys.path.append(str(damask.solver.Marc().library_path)) scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) diff --git a/processing/pre/mentat_spectralBox.py b/processing/pre/mentat_spectralBox.py index a61bef57a..e6d138952 100755 --- a/processing/pre/mentat_spectralBox.py +++ b/processing/pre/mentat_spectralBox.py @@ -9,7 +9,7 @@ import damask scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) -sys.path.append(damask.solver.Marc().libraryPath()) +sys.path.append(str(damask.solver.Marc().library_path)) #------------------------------------------------------------------------------------------------- def outMentat(cmd,locals): diff --git a/python/damask/_environment.py b/python/damask/_environment.py index df690314c..7d12050aa 100644 --- a/python/damask/_environment.py +++ b/python/damask/_environment.py @@ -1,11 +1,11 @@ import os +from pathlib import Path import tkinter class Environment: def __init__(self): """Read and provide values of DAMASK configuration.""" - self.options = self._get_options() try: tk = tkinter.Tk() self.screen_width = tk.winfo_screenwidth() @@ -15,17 +15,8 @@ class Environment: self.screen_width = 1024 self.screen_height = 768 - def relPath(self,relative = '.'): - """Return absolute path from path relative to DAMASK root.""" - return os.path.join(self.rootDir(),relative) - - - def rootDir(self): - """Return DAMASK root path.""" - return os.path.normpath(os.path.join(os.path.realpath(__file__),'../../../')) - - - def _get_options(self): + @property + def options(self): options = {} for item in ['DAMASK_NUM_THREADS', 'MSC_ROOT', @@ -34,3 +25,13 @@ class Environment: options[item] = os.environ[item] if item in os.environ else None return options + + @property + def root_dir(self): + """Return DAMASK root path.""" + return Path(__file__).parents[2] + + + # for compatibility + def rootDir(self): + return str(self.root_dir) diff --git a/python/damask/solver/_marc.py b/python/damask/solver/_marc.py index 22e9842f5..9d8b30406 100644 --- a/python/damask/solver/_marc.py +++ b/python/damask/solver/_marc.py @@ -2,6 +2,7 @@ import os import subprocess import shlex import string +from pathlib import Path from .._environment import Environment @@ -25,22 +26,22 @@ class Marc: self.version = -1 -#-------------------------- - def libraryPath(self): + @property + def library_path(self): path_MSC = Environment().options['MSC_ROOT'] - path_lib = '{}/mentat{}/shlib/linux64'.format(path_MSC,self.version) + path_lib = Path('{}/mentat{}/shlib/linux64'.format(path_MSC,self.version)) - return path_lib if os.path.exists(path_lib) else '' + return path_lib if path_lib.is_file() else None -#-------------------------- - def toolsPath(self): + @property + def tools_path(self): path_MSC = Environment().options['MSC_ROOT'] path_tools = '{}/marc{}/tools'.format(path_MSC,self.version) - return path_tools if os.path.exists(path_tools) else '' + return path_tools if path_tools.is_file() else None #-------------------------- @@ -53,21 +54,21 @@ class Marc: ): - damaskEnv = Environment() + env = Environment() - user = os.path.join(damaskEnv.relPath('src'),'DAMASK_marc{}.{}'.format(self.version,'f90' if compile else 'marc')) - if not os.path.isfile(user): + user = env.root_dir/Path('src')/Path('DAMASK_marc{}'.format(version)).with_suffix('.f90' if compile else '.marc') + if not user.is_file(): raise FileNotFoundError("DAMASK4Marc ({}) '{}' not found".format(('source' if compile else 'binary'),user)) # Define options [see Marc Installation and Operation Guide, pp 23] script = 'run_damask_{}mp'.format(optimization) - cmd = os.path.join(self.toolsPath(),script) + \ + cmd = str(self.tools_path/Path(script)) + \ ' -jid ' + model + '_' + job + \ ' -nprocd 1 -autorst 0 -ci n -cr n -dcoup 0 -b no -v no' - if compile: cmd += ' -u ' + user + ' -save y' - else: cmd += ' -prog ' + os.path.splitext(user)[0] + if compile: cmd += ' -u ' + str(user) + ' -save y' + else: cmd += ' -prog ' + str(user.with_suffix('')) print('job submission {} compilation: {}'.format('with' if compile else 'without',user)) if logfile: log = open(logfile, 'w') From 33a8aa4db72ee8bb16ba872599eb620dd6dba8ef Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Jun 2020 13:46:53 +0200 Subject: [PATCH 04/43] don't catch all exceptions and give meaningful meassages --- processing/pre/mentat_pbcOnBoxMesh.py | 23 ++++++++++------------- 1 file changed, 10 insertions(+), 13 deletions(-) diff --git a/processing/pre/mentat_pbcOnBoxMesh.py b/processing/pre/mentat_pbcOnBoxMesh.py index defcbb09b..adf0b1aff 100755 --- a/processing/pre/mentat_pbcOnBoxMesh.py +++ b/processing/pre/mentat_pbcOnBoxMesh.py @@ -197,14 +197,15 @@ def add_servoLinks(mfd_data,active=[True,True,True]): # directions on which to if mfd_data[i]['uid'] == 1705: del mfd_data[i] mfd_data.insert(i, links) - + + #-------------------------------------------------------------------------------------------------- # MAIN #-------------------------------------------------------------------------------------------------- + parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """ Set up servo linking to achieve periodic boundary conditions for a regular hexahedral mesh. Use *py_connection to operate on model presently opened in MSC.Mentat. - """, version = scriptID) parser.add_option('-p', '--port', @@ -229,10 +230,7 @@ if remote and filenames != []: if filenames == []: filenames = [None] if remote: - try: import py_mentat - except: - damask.util.croak('no valid Mentat release found.') - sys.exit(-1) + import py_mentat damask.util.report(scriptName, 'waiting to connect...') filenames = [os.path.join(tempfile._get_default_tempdir(), next(tempfile._get_candidate_names()) + '.mfd')] @@ -240,14 +238,14 @@ if remote: py_mentat.py_connect('',options.port) py_mentat.py_send('*set_save_formatted on') py_mentat.py_send('*save_as_model "{}" yes'.format(filenames[0])) - py_mentat.py_get_int("nnodes()") # hopefully blocks until file is written - except: - damask.util.croak('failed. try setting Tools/Python/"Run as Separate Process" & "Initiate".') - sys.exit() + py_mentat.py_get_int("nnodes()") + except py_mentat.InputError as err: + damask.util.croak('{}. Try Tools/Python/"Run as Separate Process" & "Initiate".'.format(err)) + sys.exit(-1) damask.util.croak( 'connected...') for name in filenames: - while remote and not os.path.exists(name): time.sleep(0.5) # wait for Mentat to write MFD file + while remote and not os.path.exists(name): time.sleep(0.5) with open( name,'r') if name is not None else sys.stdin as fileIn: damask.util.report(scriptName, name) mfd = parseMFD(fileIn) @@ -257,5 +255,4 @@ for name in filenames: fileOut.write(asMFD(mfd)) if remote: - try: py_mentat.py_send('*open_model "{}"'.format(filenames[0])) - except: damask.util.croak('lost connection on sending open command for "{}".'.format(filenames[0])) + py_mentat.py_send('*open_model "{}"'.format(filenames[0])) From 9b04a45bbd23a6e46307c0d30dbc8f0dfbc5f39d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Jun 2020 13:53:00 +0200 Subject: [PATCH 05/43] bugfix (wrong variable name) --- python/damask/solver/_marc.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/python/damask/solver/_marc.py b/python/damask/solver/_marc.py index 9d8b30406..8bee9efc0 100644 --- a/python/damask/solver/_marc.py +++ b/python/damask/solver/_marc.py @@ -1,4 +1,3 @@ -import os import subprocess import shlex import string @@ -56,7 +55,7 @@ class Marc: env = Environment() - user = env.root_dir/Path('src')/Path('DAMASK_marc{}'.format(version)).with_suffix('.f90' if compile else '.marc') + user = env.root_dir/Path('src/DAMASK_marc{}'.format(self.version)).with_suffix('.f90' if compile else '.marc') if not user.is_file(): raise FileNotFoundError("DAMASK4Marc ({}) '{}' not found".format(('source' if compile else 'binary'),user)) From 18849c3a69b411bbe5b178509f1979a2c805c959 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Jun 2020 14:03:39 +0200 Subject: [PATCH 06/43] support for older python --- installation/symlink_Processing.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/installation/symlink_Processing.py b/installation/symlink_Processing.py index 9ee6db99d..441e9a9ea 100755 --- a/installation/symlink_Processing.py +++ b/installation/symlink_Processing.py @@ -20,7 +20,7 @@ for sub_dir in ['pre','post']: for the_file in the_dir.glob('*.py'): src = the_dir/the_file dst = bin_dir/Path(the_file.with_suffix('').name) - dst.unlink(True) + if dst.is_file(): dst.unlink() # dst.unlink(True) for Python >3.8 dst.symlink_to(src) From 5da1aa49bc1778c8f38e6f555d9f1891fc977198 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Jun 2020 15:17:00 +0200 Subject: [PATCH 07/43] string should be a Path object --- python/damask/solver/_marc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/solver/_marc.py b/python/damask/solver/_marc.py index 8bee9efc0..ce32b17e5 100644 --- a/python/damask/solver/_marc.py +++ b/python/damask/solver/_marc.py @@ -38,7 +38,7 @@ class Marc: def tools_path(self): path_MSC = Environment().options['MSC_ROOT'] - path_tools = '{}/marc{}/tools'.format(path_MSC,self.version) + path_tools = Path('{}/marc{}/tools'.format(path_MSC,self.version)) return path_tools if path_tools.is_file() else None From bb5485927ed74814d4f6c1151e757e7d32a1564f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Jun 2020 15:18:46 +0200 Subject: [PATCH 08/43] names like 2020.2 cannot be converted to int --- python/damask/solver/_marc.py | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/python/damask/solver/_marc.py b/python/damask/solver/_marc.py index ce32b17e5..0484ba7e6 100644 --- a/python/damask/solver/_marc.py +++ b/python/damask/solver/_marc.py @@ -19,11 +19,7 @@ class Marc: """ self.solver = 'Marc' - try: - self.version = int(version) - except TypeError: - self.version = -1 - + self.version = version @property def library_path(self): From bda555fd1c55fcf46aea4288cd523321bb12b84c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Jun 2020 16:47:03 +0200 Subject: [PATCH 09/43] we are looking for a path, not for a file --- python/damask/solver/_marc.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/damask/solver/_marc.py b/python/damask/solver/_marc.py index 0484ba7e6..dfe45a46d 100644 --- a/python/damask/solver/_marc.py +++ b/python/damask/solver/_marc.py @@ -27,7 +27,7 @@ class Marc: path_MSC = Environment().options['MSC_ROOT'] path_lib = Path('{}/mentat{}/shlib/linux64'.format(path_MSC,self.version)) - return path_lib if path_lib.is_file() else None + return path_lib if path_lib.is_dir() else None @property @@ -36,7 +36,7 @@ class Marc: path_MSC = Environment().options['MSC_ROOT'] path_tools = Path('{}/marc{}/tools'.format(path_MSC,self.version)) - return path_tools if path_tools.is_file() else None + return path_tools if path_tools.is_dir() else None #-------------------------- From a9c61ede69e21c094f304828afbe47b31bdfd29d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Jun 2020 20:24:18 +0200 Subject: [PATCH 10/43] bugfix: should also work if DAMASK_NUM_THREADS is not set --- python/damask/_result.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/python/damask/_result.py b/python/damask/_result.py index b6ac3782e..f68dfbea0 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -1030,14 +1030,17 @@ class Result: Parameters ---------- func : function - Callback function that calculates a new dataset from one or more datasets per HDF5 group. + Callback function that calculates a new dataset from one or + more datasets per HDF5 group. datasets : dictionary - Details of the datasets to be used: label (in HDF5 file) and arg (argument to which the data is parsed in func). + Details of the datasets to be used: label (in HDF5 file) and + arg (argument to which the data is parsed in func). args : dictionary, optional Arguments parsed to func. """ - pool = multiprocessing.Pool(int(Environment().options['DAMASK_NUM_THREADS'])) + num_threads = Environment().options['DAMASK_NUM_THREADS'] + pool = multiprocessing.Pool(int(num_threads) if num_threads is not None else None) lock = multiprocessing.Manager().Lock() groups = self.groups_with_datasets(datasets.values()) From a9e0e93213314d18300935a8e6fad6b4caa9d945 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 5 Jun 2020 13:38:40 +0200 Subject: [PATCH 11/43] need to handle case of zero length file when splitting --- src/IO.f90 | 4 ++-- src/grid/DAMASK_grid.f90 | 3 ++- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/IO.f90 b/src/IO.f90 index e6e3d4a60..8d56bdd70 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -11,7 +11,6 @@ module IO implicit none private character(len=*), parameter, public :: & - IO_EOF = '#EOF#', & !< end of file string IO_WHITESPACE = achar(44)//achar(32)//achar(9)//achar(10)//achar(13) !< whitespace characters character, parameter, public :: & IO_EOL = new_line('DAMASK'), & !< end of line character @@ -86,6 +85,7 @@ function IO_readlines(fileName) result(fileContent) do l=1, len(rawData) if (rawData(l:l) == IO_EOL) N_lines = N_lines+1 enddo + if (l>1) then; if(rawData(l-1:l-1) /= IO_EOL) N_lines = N_lines+1; endif ! no EOL@EOF, need exception for empty file allocate(fileContent(N_lines)) !-------------------------------------------------------------------------------------------------- @@ -94,7 +94,7 @@ function IO_readlines(fileName) result(fileContent) startPos = 1 l = 1 do while (l <= N_lines) - endPos = startPos + scan(rawData(startPos:),IO_EOL) - 2 + endPos = merge(startPos + scan(rawData(startPos:),IO_EOL) - 2,len(rawData),l /= N_lines) if (endPos - startPos > pStringLen-1) then line = rawData(startPos:startPos+pStringLen-1) if (.not. warned) then diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index 7b3265740..8e57a56af 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -140,7 +140,8 @@ program DAMASK_grid !-------------------------------------------------------------------------------------------------- ! reading information from load case file and to sanity checks - fileContent = IO_read_ASCII(trim(loadCaseFile)) + fileContent = IO_readlines(trim(loadCaseFile)) + if(size(fileContent) == 0) call IO_error(307,ext_msg='No load case specified') allocate (loadCases(0)) ! array of load cases do currentLoadCase = 1, size(fileContent) From 9cd9ee71c5e9034e53262e76193b73d4fd34f4bc Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 5 Jun 2020 14:44:31 +0200 Subject: [PATCH 12/43] off-by-one issue fixed --- src/IO.f90 | 4 ++-- src/mesh/discretization_mesh.f90 | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/IO.f90 b/src/IO.f90 index 8d56bdd70..608d70fba 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -94,7 +94,7 @@ function IO_readlines(fileName) result(fileContent) startPos = 1 l = 1 do while (l <= N_lines) - endPos = merge(startPos + scan(rawData(startPos:),IO_EOL) - 2,len(rawData),l /= N_lines) + endPos = merge(startPos + scan(rawData(startPos:),IO_EOL) - 2,len(rawData)-1,l /= N_lines) if (endPos - startPos > pStringLen-1) then line = rawData(startPos:startPos+pStringLen-1) if (.not. warned) then @@ -106,7 +106,7 @@ function IO_readlines(fileName) result(fileContent) endif startPos = endPos + 2 ! jump to next line start - fileContent(l) = line + fileContent(l) = trim(line)//'' l = l + 1 enddo diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index 0880de115..80db37e8d 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -119,7 +119,7 @@ subroutine discretization_mesh_init(restart) call MPI_Bcast(mesh_boundaries,mesh_Nboundaries,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr) if (worldrank == 0) then - fileContent = IO_read_ASCII(geometryFile) + fileContent = IO_readlines(geometryFile) l = 0 do l = l + 1 From 5c544a6e4e9ce61eb682638c3f308ffb0f32a11d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 5 Jun 2020 20:58:36 +0200 Subject: [PATCH 13/43] bugfix - IO_read sanities files such that they end with EOL (unless 0 byte) - IO_readline simply counts EOL to determine number of lines --- src/IO.f90 | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/IO.f90 b/src/IO.f90 index 608d70fba..362f4a8d5 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -85,7 +85,6 @@ function IO_readlines(fileName) result(fileContent) do l=1, len(rawData) if (rawData(l:l) == IO_EOL) N_lines = N_lines+1 enddo - if (l>1) then; if(rawData(l-1:l-1) /= IO_EOL) N_lines = N_lines+1; endif ! no EOL@EOF, need exception for empty file allocate(fileContent(N_lines)) !-------------------------------------------------------------------------------------------------- @@ -94,7 +93,7 @@ function IO_readlines(fileName) result(fileContent) startPos = 1 l = 1 do while (l <= N_lines) - endPos = merge(startPos + scan(rawData(startPos:),IO_EOL) - 2,len(rawData)-1,l /= N_lines) + endPos = startPos + scan(rawData(startPos:),IO_EOL) - 2 if (endPos - startPos > pStringLen-1) then line = rawData(startPos:startPos+pStringLen-1) if (.not. warned) then @@ -114,7 +113,8 @@ end function IO_readlines !-------------------------------------------------------------------------------------------------- -!> @brief reads an entire ASCII file into a string +!> @brief read ASCII file into a string +!> @details ensures that the string ends with a new line (expected UNIX behavior) !-------------------------------------------------------------------------------------------------- function IO_read(fileName) result(fileContent) @@ -130,10 +130,14 @@ function IO_read(fileName) result(fileContent) status='old', position='rewind', action='read',iostat=myStat) if(myStat /= 0) call IO_error(100,ext_msg=trim(fileName)) allocate(character(len=fileLength)::fileContent) + if(fileLength==0) return + read(fileUnit,iostat=myStat) fileContent - if(myStat > 0) call IO_error(102,ext_msg=trim(fileName)) ! <0 for ifort (https://software.intel.com/en-us/comment/1960081) + if(myStat /= 0) call IO_error(102,ext_msg=trim(fileName)) close(fileUnit) + if(fileContent(fileLength:fileLength) /= IO_EOL) fileContent = fileContent//IO_EOL ! ensure EOL@EOF + end function IO_read From 0e11d6f0498fef6ba200dd3d54112a44dbdf9bc3 Mon Sep 17 00:00:00 2001 From: Test User Date: Sun, 14 Jun 2020 15:51:48 +0200 Subject: [PATCH 14/43] [skip ci] updated version information after successful test of v2.0.3-2654-g5c544a6e --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 08fd9e952..4a0b49cc6 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-2650-g512e54a7 +v2.0.3-2654-g5c544a6e From 753fbb70fd7401745f9a3259b5c124affd5e0df6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 16 Jun 2020 13:27:21 +0200 Subject: [PATCH 15/43] cleaning --- src/discretization.f90 | 2 +- src/mesh/discretization_mesh.f90 | 146 ++++++++++++++----------------- 2 files changed, 67 insertions(+), 81 deletions(-) diff --git a/src/discretization.f90 b/src/discretization.f90 index e9a987912..16f76b158 100644 --- a/src/discretization.f90 +++ b/src/discretization.f90 @@ -49,7 +49,7 @@ subroutine discretization_init(homogenizationAt,microstructureAt,& IPcoords0, & NodeCoords0 integer, optional, intent(in) :: & - sharedNodesBegin + sharedNodesBegin ! index of first node shared among different processes (MPI) write(6,'(/,a)') ' <<<+- discretization init -+>>>'; flush(6) diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index 0880de115..329205f77 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -8,53 +8,50 @@ module discretization_mesh #include #include #include - use PETScdmplex - use PETScdmda - use PETScis + use PETScdmplex + use PETScdmda + use PETScis - use DAMASK_interface - use IO - use debug - use discretization - use numerics - use FEsolving - use FEM_quadrature - use prec - - implicit none - private - - integer, public, protected :: & - mesh_Nboundaries, & - mesh_NcpElemsGlobal + use DAMASK_interface + use IO + use debug + use discretization + use numerics + use FEsolving + use FEM_quadrature + use prec - integer :: & - mesh_NcpElems !< total number of CP elements in mesh + implicit none + private + + integer, public, protected :: & + mesh_Nboundaries, & + mesh_NcpElemsGlobal + + integer :: & + mesh_NcpElems !< total number of CP elements in mesh !!!! BEGIN DEPRECATED !!!!! - integer, public, protected :: & - mesh_maxNips !< max number of IPs in any CP element + integer, public, protected :: & + mesh_maxNips !< max number of IPs in any CP element !!!! BEGIN DEPRECATED !!!!! - integer, dimension(:,:), allocatable :: & - mesh_element !DEPRECATED + real(pReal), dimension(:,:), allocatable :: & + mesh_ipVolume, & !< volume associated with IP (initially!) + mesh_node0 !< node x,y,z coordinates (initially!) - real(pReal), dimension(:,:), allocatable :: & - mesh_ipVolume, & !< volume associated with IP (initially!) - mesh_node0 !< node x,y,z coordinates (initially!) - - real(pReal), dimension(:,:,:), allocatable :: & - mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!) + real(pReal), dimension(:,:,:), allocatable :: & + mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!) - DM, public :: geomMesh - - PetscInt, dimension(:), allocatable, public, protected :: & - mesh_boundaries + DM, public :: geomMesh - public :: & - discretization_mesh_init, & - mesh_FEM_build_ipVolumes, & - mesh_FEM_build_ipCoordinates + PetscInt, dimension(:), allocatable, public, protected :: & + mesh_boundaries + + public :: & + discretization_mesh_init, & + mesh_FEM_build_ipVolumes, & + mesh_FEM_build_ipCoordinates contains @@ -67,39 +64,32 @@ subroutine discretization_mesh_init(restart) logical, intent(in) :: restart - integer, dimension(1), parameter:: FE_geomtype = [1] !< geometry type of particular element type - integer, dimension(1) :: FE_Nips !< number of IPs in a specific type of element - integer, allocatable, dimension(:) :: chunkPos integer :: dimPlex, & mesh_Nnodes, & !< total number of nodes in mesh j, l - integer, parameter :: & - mesh_ElemType=1 !< Element type of the mesh (only support homogeneous meshes) PetscSF :: sf DM :: globalMesh PetscInt :: nFaceSets PetscInt, pointer, dimension(:) :: pFaceSets character(len=pStringLen), dimension(:), allocatable :: fileContent - IS :: faceSetIS + IS :: faceSetIS PetscErrorCode :: ierr + integer, dimension(:,:), allocatable :: & + mesh_element + - write(6,'(/,a)') ' <<<+- mesh init -+>>>' - ! read in file call DMPlexCreateFromFile(PETSC_COMM_WORLD,geometryFile,PETSC_TRUE,globalMesh,ierr) CHKERRQ(ierr) - ! get spatial dimension (2 or 3?) call DMGetDimension(globalMesh,dimPlex,ierr) CHKERRQ(ierr) - write(6,*) 'dimension',dimPlex;flush(6) call DMGetStratumSize(globalMesh,'depth',dimPlex,mesh_NcpElemsGlobal,ierr) CHKERRQ(ierr) ! get number of IDs in face sets (for boundary conditions?) call DMGetLabelSize(globalMesh,'Face Sets',mesh_Nboundaries,ierr) CHKERRQ(ierr) - write(6,*) 'number of "Face Sets"',mesh_Nboundaries;flush(6) call MPI_Bcast(mesh_Nboundaries,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr) call MPI_Bcast(mesh_NcpElemsGlobal,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr) call MPI_Bcast(dimPlex,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr) @@ -142,48 +132,44 @@ subroutine discretization_mesh_init(restart) enddo call DMClone(globalMesh,geomMesh,ierr) CHKERRQ(ierr) - else + else call DMPlexDistribute(globalMesh,0,sf,geomMesh,ierr) CHKERRQ(ierr) - endif + endif call DMDestroy(globalMesh,ierr); CHKERRQ(ierr) - + call DMGetStratumSize(geomMesh,'depth',dimPlex,mesh_NcpElems,ierr) CHKERRQ(ierr) call DMGetStratumSize(geomMesh,'depth',0,mesh_Nnodes,ierr) CHKERRQ(ierr) - FE_Nips(FE_geomtype(1)) = FEM_nQuadrature(dimPlex,integrationOrder) - mesh_maxNips = FE_Nips(1) - - write(6,*) 'mesh_maxNips',mesh_maxNips + mesh_maxNips = FEM_nQuadrature(dimPlex,integrationOrder) + call mesh_FEM_build_ipCoordinates(dimPlex,FEM_quadrature_points(dimPlex,integrationOrder)%p) call mesh_FEM_build_ipVolumes(dimPlex) - + allocate (mesh_element (4,mesh_NcpElems)); mesh_element = 0 do j = 1, mesh_NcpElems - mesh_element( 1,j) = -1 ! DEPRECATED - mesh_element( 2,j) = mesh_elemType ! elem type + mesh_element( 1,j) = -1 + mesh_element( 2,j) = -1 mesh_element( 3,j) = 1 ! homogenization call DMGetLabelValue(geomMesh,'material',j-1,mesh_element(4,j),ierr) CHKERRQ(ierr) - end do + end do + + if (debug_e < 1 .or. debug_e > mesh_NcpElems) call IO_error(602,ext_msg='element') + if (debug_i < 1 .or. debug_i > mesh_maxNips) call IO_error(602,ext_msg='IP') + + FEsolving_execElem = [1,mesh_NcpElems] ! parallel loop bounds set to comprise all DAMASK elements + FEsolving_execIP = [1,mesh_maxNips] - if (debug_e < 1 .or. debug_e > mesh_NcpElems) & - call IO_error(602,ext_msg='element') ! selected element does not exist - if (debug_i < 1 .or. debug_i > FE_Nips(FE_geomtype(mesh_element(2,debug_e)))) & - call IO_error(602,ext_msg='IP') ! selected element does not have requested IP - - FEsolving_execElem = [ 1,mesh_NcpElems ] ! parallel loop bounds set to comprise all DAMASK elements - FEsolving_execIP = [1,FE_Nips(FE_geomtype(mesh_element(2,1)))] - allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal) call discretization_init(mesh_element(3,:),mesh_element(4,:),& reshape(mesh_ipCoordinates,[3,mesh_maxNips*mesh_NcpElems]), & mesh_node0) - + end subroutine discretization_mesh_init @@ -191,23 +177,23 @@ end subroutine discretization_mesh_init !> @brief Calculates IP volume. Allocates global array 'mesh_ipVolume' !-------------------------------------------------------------------------------------------------- subroutine mesh_FEM_build_ipVolumes(dimPlex) - + PetscInt,intent(in):: dimPlex PetscReal :: vol PetscReal, pointer,dimension(:) :: pCent, pNorm PetscInt :: cellStart, cellEnd, cell PetscErrorCode :: ierr - + allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems),source=0.0_pReal) - + call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr) allocate(pCent(dimPlex)) allocate(pNorm(dimPlex)) do cell = cellStart, cellEnd-1 - call DMPlexComputeCellGeometryFVM(geomMesh,cell,vol,pCent,pNorm,ierr) + call DMPlexComputeCellGeometryFVM(geomMesh,cell,vol,pCent,pNorm,ierr) CHKERRQ(ierr) mesh_ipVolume(:,cell+1) = vol/real(mesh_maxNips,pReal) - enddo + enddo end subroutine mesh_FEM_build_ipVolumes @@ -219,20 +205,20 @@ subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints) PetscInt, intent(in) :: dimPlex PetscReal, intent(in) :: qPoints(mesh_maxNips*dimPlex) - + PetscReal, pointer,dimension(:) :: pV0, pCellJ, pInvcellJ PetscReal :: detJ PetscInt :: cellStart, cellEnd, cell, qPt, dirI, dirJ, qOffset PetscErrorCode :: ierr - - + + allocate(mesh_ipCoordinates(3,mesh_maxNips,mesh_NcpElems),source=0.0_pReal) - + allocate(pV0(dimPlex)) allocatE(pCellJ(dimPlex**2)) allocatE(pinvCellJ(dimPlex**2)) call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr) - do cell = cellStart, cellEnd-1 !< loop over all elements + do cell = cellStart, cellEnd-1 !< loop over all elements call DMPlexComputeCellGeometryAffineFEM(geomMesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr) CHKERRQ(ierr) qOffset = 0 @@ -246,7 +232,7 @@ subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints) enddo qOffset = qOffset + dimPlex enddo - enddo + enddo end subroutine mesh_FEM_build_ipCoordinates From e0d4ee44a340a7adf327732c3fb9e921ee917b2a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 16 Jun 2020 13:59:59 +0200 Subject: [PATCH 16/43] better name --- src/mesh/discretization_mesh.f90 | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index 329205f77..f462f7196 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -75,8 +75,9 @@ subroutine discretization_mesh_init(restart) character(len=pStringLen), dimension(:), allocatable :: fileContent IS :: faceSetIS PetscErrorCode :: ierr - integer, dimension(:,:), allocatable :: & - mesh_element + integer, dimension(:), allocatable :: & + homogenizationAt, & + microstructureAt write(6,'(/,a)') ' <<<+- mesh init -+>>>' @@ -149,12 +150,10 @@ subroutine discretization_mesh_init(restart) call mesh_FEM_build_ipCoordinates(dimPlex,FEM_quadrature_points(dimPlex,integrationOrder)%p) call mesh_FEM_build_ipVolumes(dimPlex) - allocate (mesh_element (4,mesh_NcpElems)); mesh_element = 0 + allocate(microstructureAt(mesh_NcpElems)) + allocate(homogenizationAt(mesh_NcpElems),source=1) do j = 1, mesh_NcpElems - mesh_element( 1,j) = -1 - mesh_element( 2,j) = -1 - mesh_element( 3,j) = 1 ! homogenization - call DMGetLabelValue(geomMesh,'material',j-1,mesh_element(4,j),ierr) + call DMGetLabelValue(geomMesh,'material',j-1,microstructureAt(j),ierr) CHKERRQ(ierr) end do @@ -166,7 +165,7 @@ subroutine discretization_mesh_init(restart) allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal) - call discretization_init(mesh_element(3,:),mesh_element(4,:),& + call discretization_init(microstructureAt,homogenizationAt,& reshape(mesh_ipCoordinates,[3,mesh_maxNips*mesh_NcpElems]), & mesh_node0) From d0d9245707a665952042db5a2d2090af8a574939 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Jun 2020 15:02:39 +0200 Subject: [PATCH 17/43] clearer intention --- python/tests/test_Result.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py index d7946e5e0..38c0d84cc 100644 --- a/python/tests/test_Result.py +++ b/python/tests/test_Result.py @@ -302,8 +302,8 @@ class TestResult: @pytest.mark.parametrize('allowed',['off','on']) def test_rename(self,default,allowed): - F = default.read_dataset(default.get_dataset_location('F')) if allowed == 'on': + F = default.read_dataset(default.get_dataset_location('F')) default.allow_modification() default.rename('F','new_name') assert np.all(F == default.read_dataset(default.get_dataset_location('new_name'))) From 829896390ca3089d2f9e2e1d7d9d3b2659f2878e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 16 Jun 2020 23:37:24 +0200 Subject: [PATCH 18/43] hopefully not needed any more --- src/crystallite.f90 | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 929fec862..e85f9e517 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -287,11 +287,9 @@ end subroutine crystallite_init !-------------------------------------------------------------------------------------------------- !> @brief calculate stress (P) !-------------------------------------------------------------------------------------------------- -function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) +function crystallite_stress() logical, dimension(discretization_nIP,discretization_nElem) :: crystallite_stress - real(pReal), intent(in), optional :: & - dummyArgumentToPreventInternalCompilerErrorWithGCC real(pReal) :: & formerSubStep integer :: & From 62384b583677347de327158d431511f5ee6347da Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 16 Jun 2020 23:43:15 +0200 Subject: [PATCH 19/43] bugfix: invalid description/unit --- src/constitutive_plastic_dislotwin.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/constitutive_plastic_dislotwin.f90 b/src/constitutive_plastic_dislotwin.f90 index 7c7d24ab8..dc40d269e 100644 --- a/src/constitutive_plastic_dislotwin.f90 +++ b/src/constitutive_plastic_dislotwin.f90 @@ -830,7 +830,7 @@ module subroutine plastic_dislotwin_results(instance,group) 'mobile dislocation density','1/m²') case('rho_dip') if(prm%sum_N_sl>0) call results_writeDataset(group,stt%rho_dip,'rho_dip',& - 'dislocation dipole density''1/m²') + 'dislocation dipole density','1/m²') case('gamma_sl') if(prm%sum_N_sl>0) call results_writeDataset(group,stt%gamma_sl,'gamma_sl',& 'plastic shear','1') From 8fa68101b8e8c13774c4b28be67b3f49dcefd65f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 17 Jun 2020 22:28:29 +0200 Subject: [PATCH 20/43] not needed --- examples/ConfigFiles/Phase_DisloUCLA_Tungsten.config | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/ConfigFiles/Phase_DisloUCLA_Tungsten.config b/examples/ConfigFiles/Phase_DisloUCLA_Tungsten.config index 6920aee2b..2cdd5c58f 100644 --- a/examples/ConfigFiles/Phase_DisloUCLA_Tungsten.config +++ b/examples/ConfigFiles/Phase_DisloUCLA_Tungsten.config @@ -27,7 +27,7 @@ SolidSolutionStrength 0.0 # Strength due to elements in solid solution ### Dislocation glide parameters ### #per family -Nslip 12 0 +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] From 01f21a9e3baa5aca8fe8eec0421bdd76197e5675 Mon Sep 17 00:00:00 2001 From: Test User Date: Thu, 18 Jun 2020 18:04:51 +0200 Subject: [PATCH 21/43] [skip ci] updated version information after successful test of v2.0.3-2664-ge959aaab --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 4a0b49cc6..b25957f75 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-2654-g5c544a6e +v2.0.3-2664-ge959aaab From 7c4afe06c9cfe9db116d0ab75b949a50b44202a3 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Fri, 19 Jun 2020 02:32:33 +0200 Subject: [PATCH 22/43] less generic variables --- PRIVATE | 2 +- src/grid/DAMASK_grid.f90 | 11 ++++----- src/grid/grid_damage_spectral.f90 | 12 +++++----- src/grid/grid_mech_FEM.f90 | 22 ++++++++---------- src/grid/grid_mech_spectral_basic.f90 | 24 ++++++++------------ src/grid/grid_mech_spectral_polarisation.f90 | 24 +++++++------------- src/grid/grid_thermal_spectral.f90 | 12 +++++----- src/grid/spectral_utilities.f90 | 9 +++----- src/mesh/DAMASK_mesh.f90 | 8 +++---- src/mesh/FEM_utilities.f90 | 7 ++---- src/mesh/mesh_mech_FEM.f90 | 13 +++++------ 11 files changed, 60 insertions(+), 84 deletions(-) diff --git a/PRIVATE b/PRIVATE index 6f818e871..f391a36c2 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 6f818e871d6415fa4ebe2a5781920bad614d02e8 +Subproject commit f391a36c2e829f2cb4db2297271941168f7fd7d8 diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index c2ac2ea58..e151222c2 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -92,8 +92,7 @@ program DAMASK_grid external :: & quit class (tNode), pointer :: & - num_grid, & - num_generic + num_grid !-------------------------------------------------------------------------------------------------- ! init DAMASK (all modules) @@ -114,9 +113,9 @@ program DAMASK_grid !------------------------------------------------------------------------------------------------- ! reading field paramters from numerics file and do sanity checks - num_generic => numerics_root%get('generic', defaultVal=emptyDict) - stagItMax = num_generic%get_asInt('maxStaggeredIter',defaultVal=10) - maxCutBack = num_generic%get_asInt('maxCutBack',defaultVal=3) + num_grid => numerics_root%get('grid', defaultVal=emptyDict) + stagItMax = num_grid%get_asInt('maxStaggeredIter',defaultVal=10) + maxCutBack = num_grid%get_asInt('maxCutBack',defaultVal=3) if (stagItMax < 0) call IO_error(301,ext_msg='maxStaggeredIter') if (maxCutBack < 0) call IO_error(301,ext_msg='maxCutBack') @@ -124,8 +123,6 @@ program DAMASK_grid !-------------------------------------------------------------------------------------------------- ! assign mechanics solver depending on selected type - num_grid => numerics_root%get('grid',defaultVal=emptyDict) - select case (trim(num_grid%get_asString('solver', defaultVal = 'Basic'))) case ('Basic') mech_init => grid_mech_spectral_basic_init diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index ffc3b4cc2..659b111a3 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -61,7 +61,7 @@ subroutine grid_damage_spectral_init Vec :: uBound, lBound PetscErrorCode :: ierr class(tNode), pointer :: & - num_generic + num_grid character(len=pStringLen) :: & snes_type, & petsc_options @@ -73,8 +73,8 @@ subroutine grid_damage_spectral_init !------------------------------------------------------------------------------------------------- ! read numerical parameter - num_generic => numerics_root%get('generic',defaultVal=emptyDict) - petsc_options = num_generic%get_asString('petsc_options',defaultVal='') + num_grid => numerics_root%get('grid',defaultVal=emptyDict) + petsc_options = num_grid%get_asString('petsc_options',defaultVal='') !-------------------------------------------------------------------------------------------------- ! set default and user defined options for PETSc @@ -150,7 +150,7 @@ function grid_damage_spectral_solution(timeinc,timeinc_old) result(solution) itmax !< maximum number of iterations type(tSolutionState) :: solution class(tNode), pointer :: & - num_generic + num_grid PetscInt :: devNull PetscReal :: phi_min, phi_max, stagNorm, solnNorm @@ -159,8 +159,8 @@ function grid_damage_spectral_solution(timeinc,timeinc_old) result(solution) !------------------------------------------------------------------- ! reading numerical parameter and do sanity check - num_generic => numerics_root%get('generic',defaultVal=emptyDict) - itmax = num_generic%get_asInt('itmax',defaultVal=250) + num_grid => numerics_root%get('grid',defaultVal=emptyDict) + itmax = num_grid%get_asInt('itmax',defaultVal=250) if (itmax <= 1) call IO_error(301,ext_msg='itmax') solution%converged =.false. diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 03709040f..bdd76b423 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -100,7 +100,7 @@ subroutine grid_mech_FEM_init fileName, & petsc_options class(tNode), pointer :: & - num_generic + num_grid real(pReal), dimension(3,3,3,3) :: devNull PetscScalar, pointer, dimension(:,:,:,:) :: & u_current,u_lastInc @@ -109,8 +109,8 @@ subroutine grid_mech_FEM_init !------------------------------------------------------------------------------------------------- ! read numerical parameter - num_generic => numerics_root%get('generic',defaultVal=emptyDict) - petsc_options = num_generic%get_asString('petsc_options',defaultVal='') + num_grid => numerics_root%get('grid',defaultVal=emptyDict) + petsc_options = num_grid%get_asString('petsc_options',defaultVal='') !-------------------------------------------------------------------------------------------------- ! set default and user defined options for PETSc @@ -431,8 +431,7 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,i eps_stress_atol, & !< absolute tolerance for fullfillment of stress BC eps_stress_rtol !< relative tolerance for fullfillment of stress BC class(tNode), pointer :: & - num_grid, & - num_generic + num_grid !----------------------------------------------------------------------------------- ! reading numerical parameters and do sanity check @@ -442,9 +441,8 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,i eps_stress_atol = num_grid%get_asFloat('eps_stress_atol',defaultVal=1.0e3_pReal) eps_stress_rtol = num_grid%get_asFloat('eps_stress_rtol',defaultVal=0.01_pReal) - num_generic => numerics_root%get('generic',defaultVal=emptyDict) - itmin = num_generic%get_asInt('itmin',defaultVal=1) - itmax = num_generic%get_asInt('itmax',defaultVal=250) + itmin = num_grid%get_asInt('itmin',defaultVal=1) + itmax = num_grid%get_asInt('itmax',defaultVal=250) if (eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol') if (eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol') @@ -505,13 +503,13 @@ subroutine formResidual(da_local,x_local, & itmin, & itmax class(tNode), pointer :: & - num_generic + num_grid !---------------------------------------------------------------------- ! read numerical paramteters and do sanity checks - num_generic => numerics_root%get('generic',defaultVal=emptyDict) - itmin = num_generic%get_asInt('itmin',defaultVal=1) - itmax = num_generic%get_asInt('itmax',defaultVal=250) + num_grid => numerics_root%get('grid',defaultVal=emptyDict) + itmin = num_grid%get_asInt('itmin',defaultVal=1) + itmax = num_grid%get_asInt('itmax',defaultVal=250) if (itmax <= 1) call IO_error(301,ext_msg='itmax') if (itmin > itmax .or. itmin < 1) call IO_error(301,ext_msg='itmin') diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index e754b0d64..61ac18779 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -89,8 +89,7 @@ subroutine grid_mech_spectral_basic_init real(pReal), dimension(3,3) :: & temp33_Real = 0.0_pReal class (tNode), pointer :: & - num_grid, & - num_generic + num_grid PetscErrorCode :: ierr PetscScalar, pointer, dimension(:,:,:,:) :: & @@ -112,10 +111,8 @@ subroutine grid_mech_spectral_basic_init !------------------------------------------------------------------------------------------------- ! read numerical parameters - num_generic => numerics_root%get('generic',defaultVal=emptyDict) - petsc_options = num_generic%get_asString('petsc_options',defaultVal='') - num_grid => numerics_root%get('grid',defaultVal=emptyDict) + petsc_options = num_grid%get_asString('petsc_options',defaultVal='') num%update_gamma = num_grid%get_asBool('update_gamma',defaultVal=.false.) !-------------------------------------------------------------------------------------------------- @@ -403,8 +400,8 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm eps_stress_atol, & !< absolute tolerance for fullfillment of stress BC eps_stress_rtol !< relative tolerance for fullfillment of stress BC class(tNode), pointer :: & - num_grid, & - num_generic + num_grid + !----------------------------------------------------------------------------------- ! reading numerical parameters and do sanity check num_grid => numerics_root%get('grid',defaultVal=emptyDict) @@ -413,9 +410,8 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm eps_stress_atol = num_grid%get_asFloat('eps_stress_atol',defaultVal=1.0e3_pReal) eps_stress_rtol = num_grid%get_asFloat('eps_stress_rtol',defaultVal=0.01_pReal) - num_generic => numerics_root%get('generic',defaultVal=emptyDict) - itmin = num_generic%get_asInt('itmin',defaultVal=1) - itmax = num_generic%get_asInt('itmax',defaultVal=250) + itmin = num_grid%get_asInt('itmin',defaultVal=1) + itmax = num_grid%get_asInt('itmax',defaultVal=250) if (eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol') if (eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol') @@ -474,13 +470,13 @@ subroutine formResidual(in, F, & itmin, & itmax class(tNode), pointer :: & - num_generic + num_grid !---------------------------------------------------------------------- ! read numerical paramteter and do sanity checks - num_generic => numerics_root%get('generic',defaultVal=emptyDict) - itmin = num_generic%get_asInt('itmin',defaultVal=1) - itmax = num_generic%get_asInt('itmax',defaultVal=250) + num_grid => numerics_root%get('grid',defaultVal=emptyDict) + itmin = num_grid%get_asInt('itmin',defaultVal=1) + itmax = num_grid%get_asInt('itmax',defaultVal=250) if (itmax <= 1) call IO_error(301,ext_msg='itmax') if (itmin > itmax .or. itmin < 1) call IO_error(301,ext_msg='itmin') diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 33b84c340..a32216615 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -96,8 +96,7 @@ subroutine grid_mech_spectral_polarisation_init real(pReal), dimension(3,3) :: & temp33_Real = 0.0_pReal class (tNode), pointer :: & - num_grid, & - num_generic + num_grid PetscErrorCode :: ierr PetscScalar, pointer, dimension(:,:,:,:) :: & @@ -118,10 +117,8 @@ subroutine grid_mech_spectral_polarisation_init !------------------------------------------------------------------------------------------------- ! read numerical parameters - num_generic => numerics_root%get('generic',defaultVal=emptyDict) - petsc_options = num_generic%get_asString('petsc_options',defaultVal='') - num_grid => numerics_root%get('grid',defaultVal=emptyDict) + petsc_options = num_grid%get_asString('petsc_options',defaultVal='') num%update_gamma = num_grid%get_asBool('update_gamma',defaultVal=.false.) !-------------------------------------------------------------------------------------------------- @@ -453,8 +450,7 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm eps_stress_atol, & !< absolute tolerance for fullfillment of stress BC eps_stress_rtol !< relative tolerance for fullfillment of stress BC class(tNode), pointer :: & - num_grid, & - num_generic + num_grid !----------------------------------------------------------------------------------- ! reading numerical parameters and do sanity check @@ -466,9 +462,8 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm eps_stress_atol = num_grid%get_asFloat('eps_stress_atol',defaultVal=1.0e3_pReal) eps_stress_rtol = num_grid%get_asFloat('eps_stress_rtol',defaultVal=0.01_pReal) - num_generic => numerics_root%get('generic',defaultVal=emptyDict) - itmin = num_generic%get_asInt('itmin',defaultVal=1) - itmax = num_generic%get_asInt('itmax',defaultVal=250) + itmin = num_grid%get_asInt('itmin',defaultVal=1) + itmax = num_grid%get_asInt('itmax',defaultVal=250) if (eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol') if (eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol') @@ -533,8 +528,7 @@ subroutine formResidual(in, FandF_tau, & PetscObject :: dummy PetscErrorCode :: ierr class(tNode), pointer :: & - num_grid, & - num_generic + num_grid real(pReal) :: & polarAlpha, & !< polarization scheme parameter 0.0 < alpha < 2.0. alpha = 1.0 ==> AL scheme, alpha = 2.0 ==> accelerated scheme polarBeta !< polarization scheme parameter 0.0 < beta < 2.0. beta = 1.0 ==> AL scheme, beta = 2.0 ==> accelerated scheme @@ -547,10 +541,8 @@ subroutine formResidual(in, FandF_tau, & num_grid => numerics_root%get('grid',defaultVal = emptyDict) polarAlpha = num_grid%get_asFloat('polaralpha',defaultVal=1.0_pReal) polarBeta = num_grid%get_asFloat('polarbeta', defaultVal=1.0_pReal) - - num_generic => numerics_root%get('generic',defaultVal=emptyDict) - itmin = num_generic%get_asInt('itmin',defaultVal=1) - itmax = num_generic%get_asInt('itmax',defaultVal=250) + itmin = num_grid%get_asInt('itmin',defaultVal=1) + itmax = num_grid%get_asInt('itmax',defaultVal=250) if (itmax <= 1) call IO_error(301,ext_msg='itmax') if (itmin > itmax .or. itmin < 1) call IO_error(301,ext_msg='itmin') diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 418330411..ee64566a6 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -63,7 +63,7 @@ subroutine grid_thermal_spectral_init PetscScalar, dimension(:,:,:), pointer :: x_scal PetscErrorCode :: ierr class(tNode), pointer :: & - num_generic + num_grid character(len=pStringLen) :: & petsc_options @@ -74,8 +74,8 @@ subroutine grid_thermal_spectral_init !------------------------------------------------------------------------------------------------- ! read numerical parameter - num_generic => numerics_root%get('generic',defaultVal=emptyDict) - petsc_options = num_generic%get_asString('petsc_options',defaultVal='') + num_grid => numerics_root%get('grid',defaultVal=emptyDict) + petsc_options = num_grid%get_asString('petsc_options',defaultVal='') !-------------------------------------------------------------------------------------------------- ! set default and user defined options for PETSc @@ -147,7 +147,7 @@ function grid_thermal_spectral_solution(timeinc,timeinc_old) result(solution) itmax !< maximum number of iterations type(tSolutionState) :: solution class(tNode), pointer :: & - num_generic + num_grid PetscInt :: devNull PetscReal :: T_min, T_max, stagNorm, solnNorm @@ -156,8 +156,8 @@ function grid_thermal_spectral_solution(timeinc,timeinc_old) result(solution) !------------------------------------------------------------------- ! reading numerical parameter and do sanity check - num_generic => numerics_root%get('generic',defaultVal=emptyDict) - itmax = num_generic%get_asInt('itmax',defaultVal=250) + num_grid => numerics_root%get('grid',defaultVal=emptyDict) + itmax = num_grid%get_asInt('itmax',defaultVal=250) if (itmax <= 1) call IO_error(301,ext_msg='itmax') solution%converged =.false. diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index d8e922343..f09759009 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -192,8 +192,7 @@ subroutine spectral_utilities_init character(len=pStringLen) :: & petsc_options class (tNode) , pointer :: & - num_grid, & - num_generic + num_grid write(6,'(/,a)') ' <<<+- spectral_utilities init -+>>>' @@ -220,8 +219,8 @@ subroutine spectral_utilities_init trim(PETScDebug), & ' add more using the PETSc_Options keyword in numerics.config '; flush(6) - num_generic => numerics_root%get('generic',defaultVal=emptyDict) - petsc_options = num_generic%get_asString('petsc_options',defaultVal='') + num_grid => numerics_root%get('grid',defaultVal=emptyDict) + petsc_options = num_grid%get_asString('petsc_options',defaultVal='') call PETScOptionsClear(PETSC_NULL_OPTIONS,ierr) CHKERRQ(ierr) @@ -236,8 +235,6 @@ subroutine spectral_utilities_init write(6,'(/,a,3(i12 ))') ' grid a b c: ', grid write(6,'(a,3(es12.5))') ' size x y z: ', geomSize - num_grid => numerics_root%get('grid',defaultVal=emptyDict) - num%memory_efficient = num_grid%get_asInt ('memory_efficient', defaultVal=1) > 0 num%FFTW_timelimit = num_grid%get_asFloat ('fftw_timelimit', defaultVal=-1.0_pReal) num%divergence_correction = num_grid%get_asInt ('divergence_correction', defaultVal=2) diff --git a/src/mesh/DAMASK_mesh.f90 b/src/mesh/DAMASK_mesh.f90 index 0e4a02cf2..72fa53566 100644 --- a/src/mesh/DAMASK_mesh.f90 +++ b/src/mesh/DAMASK_mesh.f90 @@ -60,7 +60,7 @@ program DAMASK_mesh stagItMax, & !< max number of field level staggered iterations component class(tNode), pointer :: & - num_generic + num_mesh character(len=pStringLen), dimension(:), allocatable :: fileContent character(len=pStringLen) :: & incInfo, & @@ -80,9 +80,9 @@ program DAMASK_mesh !--------------------------------------------------------------------- ! reading field information from numerics file and do sanity checks - num_generic => numerics_root%get('generic', defaultVal=emptyDict) - stagItMax = num_generic%get_asInt('maxStaggeredIter',defaultVal=10) - maxCutBack = num_generic%get_asInt('maxCutBack',defaultVal=3) + num_mesh => numerics_root%get('mesh', defaultVal=emptyDict) + stagItMax = num_mesh%get_asInt('maxStaggeredIter',defaultVal=10) + maxCutBack = num_mesh%get_asInt('maxCutBack',defaultVal=3) if (stagItMax < 0) call IO_error(301,ext_msg='maxStaggeredIter') if (maxCutBack < 0) call IO_error(301,ext_msg='maxCutBack') diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index f7e00f42c..a8479913e 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -103,8 +103,7 @@ subroutine FEM_utilities_init character(len=pStringLen) :: petsc_optionsOrder class(tNode), pointer :: & - num_mesh, & - num_generic + num_mesh integer :: structOrder !< order of displacement shape functions character(len=pStringLen) :: & petsc_options @@ -114,9 +113,7 @@ subroutine FEM_utilities_init num_mesh => numerics_root%get('mesh',defaultVal=emptyDict) structOrder = num_mesh%get_asInt('structOrder',defaultVal = 2) - - num_generic => numerics_root%get('generic',defaultVal=emptyDict) - petsc_options = num_generic%get_asString('petsc_options', defaultVal='') + petsc_options = num_mesh%get_asString('petsc_options', defaultVal='') !-------------------------------------------------------------------------------------------------- ! set debugging parameters diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index a094a4d38..715b02c0d 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -96,8 +96,7 @@ subroutine FEM_mech_init(fieldBC) PetscErrorCode :: ierr class(tNode), pointer :: & - num_mesh, & - num_generic + num_mesh integer :: & integrationOrder, & !< order of quadrature rule required itmax @@ -109,8 +108,8 @@ subroutine FEM_mech_init(fieldBC) num_mesh => numerics_root%get('mesh',defaultVal=emptyDict) integrationOrder = num_mesh%get_asInt('integrationorder',defaultVal = 2) - num_generic => numerics_root%get('generic',defaultVal=emptyDict) - itmax = num_generic%get_asInt('itmax',defaultVal=250) + num_mesh => numerics_root%get('mesh',defaultVal=emptyDict) + itmax = num_mesh%get_asInt('itmax',defaultVal=250) if (itmax <= 1) call IO_error(301,ext_msg='itmax') !-------------------------------------------------------------------------------------------------- @@ -286,15 +285,15 @@ type(tSolutionState) function FEM_mech_solution( & integer :: & itmax class(tNode), pointer :: & - num_generic + num_mesh PetscErrorCode :: ierr SNESConvergedReason :: reason !-------------------------------------------------------------------------------------------------- ! read numerical parameter and do sanity check - num_generic => numerics_root%get('generic',defaultVal=emptyDict) - itmax = num_generic%get_asInt('itmax',defaultVal=250) + num_mesh => numerics_root%get('mesh',defaultVal=emptyDict) + itmax = num_mesh%get_asInt('itmax',defaultVal=250) if (itmax <= 1) call IO_error(301,ext_msg='itmax') !------------------------------------------------------------------------------------------------- From 055fa64f5f0270feb4f41e4384586b6a70928f57 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 19 Jun 2020 12:25:46 +0200 Subject: [PATCH 23/43] better readable --- python/damask/_rotation.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index b1b1cd5ad..7e9a03c3a 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -530,7 +530,7 @@ class Rotation: np.zeros(qu.shape[:-1]+(2,))]), np.where(np.abs(q03_s) < 1.0e-8, np.block([np.arctan2( 2.0*qu[...,1:2]*qu[...,2:3],qu[...,1:2]**2-qu[...,2:3]**2), - np.broadcast_to(np.pi,qu.shape[:-1]+(1,)), + np.broadcast_to(np.pi,qu[...,0:1].shape), np.zeros(qu.shape[:-1]+(1,))]), np.block([np.arctan2((-_P*q02+q13)*chi, (-_P*q01-q23)*chi), np.arctan2( 2.0*chi, q03_s-q12_s ), @@ -553,7 +553,7 @@ class Rotation: s = np.sign(qu[...,0:1])/np.sqrt(qu[...,1:2]**2+qu[...,2:3]**2+qu[...,3:4]**2) omega = 2.0 * np.arccos(np.clip(qu[...,0:1],-1.0,1.0)) ax = np.where(np.broadcast_to(qu[...,0:1] < 1.0e-8,qu.shape), - np.block([qu[...,1:4],np.broadcast_to(np.pi,qu.shape[:-1]+(1,))]), + np.block([qu[...,1:4],np.broadcast_to(np.pi,qu[...,0:1].shape)]), np.block([qu[...,1:4]*s,omega])) ax[np.isclose(qu[...,0],1.,rtol=0.0)] = [0.0, 0.0, 1.0, 0.0] return ax @@ -564,7 +564,7 @@ class Rotation: with np.errstate(invalid='ignore',divide='ignore'): s = np.linalg.norm(qu[...,1:4],axis=-1,keepdims=True) ro = np.where(np.broadcast_to(np.abs(qu[...,0:1]) < 1.0e-12,qu.shape), - np.block([qu[...,1:2], qu[...,2:3], qu[...,3:4], np.broadcast_to(np.inf,qu.shape[:-1]+(1,))]), + np.block([qu[...,1:2], qu[...,2:3], qu[...,3:4], np.broadcast_to(np.inf,qu[...,0:1].shape)]), np.block([qu[...,1:2]/s,qu[...,2:3]/s,qu[...,3:4]/s, np.tan(np.arccos(np.clip(qu[...,0:1],-1.0,1.0))) ]) From d93ed2bc5c9dfa4a48c4539c74d4dee48087c6db Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 20 Jun 2020 12:20:43 +0200 Subject: [PATCH 24/43] several improvements: - vectorized from_directions - more tests (96% coverage, only random functionality is untested) - updated documentation, folloing numpy standard - inverse operator '~' introduced --- python/damask/_rotation.py | 293 ++++++++++++++++++++++++++-------- python/tests/test_Rotation.py | 41 ++++- 2 files changed, 263 insertions(+), 71 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 7e9a03c3a..a33b4bc51 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -20,8 +20,8 @@ class Rotation: when viewing from the end point of the rotation axis towards the origin. - rotations will be interpreted in the passive sense. - Euler angle triplets are implemented using the Bunge convention, - with the angular ranges as [0, 2π],[0, π],[0, 2π]. - - the rotation angle ω is limited to the interval [0, π]. + with the angular ranges as [0,2π], [0,π], [0,2π]. + - the rotation angle ω is limited to the interval [0,π]. - the real part of a quaternion is positive, Re(q) > 0 - P = -1 (as default). @@ -49,7 +49,8 @@ class Rotation: Parameters ---------- quaternion : numpy.ndarray, optional - Unit quaternion that follows the conventions. Use .from_quaternion to perform a sanity check. + Unit quaternion in positive real hemisphere. + Use .from_quaternion to perform a sanity check. """ self.quaternion = quaternion.copy() @@ -73,7 +74,7 @@ class Rotation: raise NotImplementedError('Support for multiple rotations missing') return '\n'.join([ 'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)), - 'Matrix:\n{}'.format(self.as_matrix()), + 'Matrix:\n{}'.format(np.round(self.as_matrix(),8)), 'Bunge Eulers / deg: ({:3.2f}, {:3.2f}, {:3.2f})'.format(*self.as_Eulers(degrees=True)), ]) @@ -87,10 +88,6 @@ class Rotation: other : numpy.ndarray or Rotation Vector, second or fourth order tensor, or rotation object that is rotated. - Todo - ---- - Check rotation of 4th order tensor - """ if isinstance(other, Rotation): q_m = self.quaternion[...,0:1] @@ -99,7 +96,7 @@ class Rotation: p_o = other.quaternion[...,1:] q = (q_m*q_o - np.einsum('...i,...i',p_m,p_o).reshape(self.shape+(1,))) p = q_m*p_o + q_o*p_m + _P * np.cross(p_m,p_o) - return self.__class__(np.block([q,p])).standardize() + return self.__class__(np.block([q,p]))._standardize() elif isinstance(other,np.ndarray): if self.shape + (3,) == other.shape: @@ -124,24 +121,24 @@ class Rotation: else: raise TypeError('Cannot rotate {}'.format(type(other))) - def inverse(self): - """In-place inverse rotation/backward rotation.""" - self.quaternion[...,1:] *= -1 - return self - def inversed(self): - """Inverse rotation/backward rotation.""" - return self.copy().inverse() - - - def standardize(self): - """In-place quaternion representation with positive real part.""" + def _standardize(self): + """Standardize (ensure positive real hemisphere).""" self.quaternion[self.quaternion[...,0] < 0.0] *= -1 return self - def standardized(self): - """Quaternion representation with positive real part.""" - return self.copy().standardize() + def inverse(self): + """In-place inverse rotation (backward rotation).""" + self.quaternion[...,1:] *= -1 + return self + + def __invert__(self): + """Inverse rotation (backward rotation).""" + return self.copy().inverse() + + def inversed(self): + """Inverse rotation (backward rotation).""" + return ~ self def misorientation(self,other): @@ -154,7 +151,7 @@ class Rotation: Rotation to which the misorientation is computed. """ - return other*self.inversed() + return other@~self def broadcast_to(self,shape): @@ -169,7 +166,7 @@ class Rotation: return self.__class__(q) - def average(self,other): + def average(self,other): #ToDo: discuss calling for vectors """ Calculate the average rotation. @@ -189,25 +186,31 @@ class Rotation: def as_quaternion(self): """ - Unit quaternion [q, p_1, p_2, p_3]. + Represent as unit quaternion. - Parameters - ---------- - quaternion : bool, optional - return quaternion as DAMASK object. + Returns + ------- + q : numpy.ndarray of shape (...,4) + Unit quaternion in positive real hemisphere: (q_0, q_1, q_2, q_3), |q|=1, q_0 ≥ 0. """ - return self.quaternion + return self.quaternion.copy() def as_Eulers(self, degrees = False): """ - Bunge-Euler angles: (φ_1, ϕ, φ_2). + Represent as Bunge-Euler angles. Parameters ---------- degrees : bool, optional - return angles in degrees. + Return angles in degrees. + + Returns + ------- + phi : numpy.ndarray of shape (...,3) + Bunge-Euler angles: (φ_1, ϕ, φ_2), φ_1 ∈ [0,2π], ϕ ∈ [0,π], φ_2 ∈ [0,2π] + unless degrees == True: φ_1 ∈ [0,360], ϕ ∈ [0,180], φ_2 ∈ [0,360] """ eu = Rotation._qu2eu(self.quaternion) @@ -218,14 +221,21 @@ class Rotation: degrees = False, pair = False): """ - Axis angle representation [n_1, n_2, n_3, ω] unless pair == True: ([n_1, n_2, n_3], ω). + Represent as axis angle pair. Parameters ---------- degrees : bool, optional - return rotation angle in degrees. + Return rotation angle in degrees. Defaults to False. pair : bool, optional - return tuple of axis and angle. + Return tuple of axis and angle. Defaults to False. + + Returns + ------- + axis_angle : numpy.ndarray of shape (...,4) unless pair == True: + tuple containing numpy.ndarray of shapes (...,3) and (...) + Axis angle pair: (n_1, n_2, n_3, ω), |n| = 1 and ω ∈ [0,π] + unless degrees = True: ω ∈ [0,180]. """ ax = Rotation._qu2ax(self.quaternion) @@ -233,29 +243,60 @@ class Rotation: return (ax[...,:3],ax[...,3]) if pair else ax def as_matrix(self): - """Rotation matrix.""" + """ + Represent as rotation matrix. + + Returns + ------- + R : numpy.ndarray of shape (...,3,3) + Rotation matrix R, det(R) = 1, R.T∙R=I. + + """ return Rotation._qu2om(self.quaternion) def as_Rodrigues(self, vector = False): """ - Rodrigues-Frank vector representation [n_1, n_2, n_3, tan(ω/2)] unless vector == True: [n_1, n_2, n_3] * tan(ω/2). + Represent as Rodrigues-Frank vector with separated axis and angle argument. Parameters ---------- vector : bool, optional - return as actual Rodrigues--Frank vector, i.e. rotation axis scaled by tan(ω/2). + 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,π]. """ ro = Rotation._qu2ro(self.quaternion) return ro[...,:3]*ro[...,3] if vector else ro def as_homochoric(self): - """Homochoric vector: (h_1, h_2, h_3).""" + """ + Represent as homochoric vector. + + Returns + ------- + h : numpy.ndarray of shape (...,3) + Homochoric vector: (h_1, h_2, h_3), |h| < 1/2*π^(2/3). + + """ return Rotation._qu2ho(self.quaternion) def as_cubochoric(self): - """Cubochoric vector: (c_1, c_2, c_3).""" + """ + Represent as cubocoric vector. + + Returns + ------- + c : numpy.ndarray of shape (...,3) + Cubochoric vector: (c_1, c_2, c_3), max(c_i) < 1/2*π^(2/3). + + """ return Rotation._qu2cu(self.quaternion) def M(self): # ToDo not sure about the name: as_M or M? we do not have a from_M @@ -275,18 +316,34 @@ class Rotation: # Static constructors. The input data needs to follow the conventions, options allow to # relax the conventions. @staticmethod - def from_quaternion(quaternion, + def from_quaternion(q, accept_homomorph = False, P = -1, - acceptHomomorph = None): + acceptHomomorph = None): # old name (for compatibility) + """ + Initialize from quaternion. + Parameters + ---------- + q: numpy.ndarray of shape (...,4) + Unit quaternion in positive real hemisphere: (q_0, q_1, q_2, q_3), + |q|=1, q_0 ≥ 0. + accept_homomorph: boolean, optional + Allow homomorphic variants, i.e. q_0 < 0 (negative real hemisphere). + Defaults to False. + P: integer ∈ {-1,1}, optional + Convention used. Defaults to -1. + + """ if acceptHomomorph is not None: - accept_homomorph = acceptHomomorph - qu = np.array(quaternion,dtype=float) + accept_homomorph = acceptHomomorph # for compatibility + qu = np.array(q,dtype=float) if qu.shape[:-2:-1] != (4,): raise ValueError('Invalid shape.') + if abs(P) != 1: + raise ValueError('P ∉ {-1,1}') - if P > 0: qu[...,1:4] *= -1 # convert from P=1 to P=-1 + if P == 1: qu[...,1:4] *= -1 if accept_homomorph: qu[qu[...,0] < 0.0] *= -1 else: @@ -298,10 +355,21 @@ class Rotation: return Rotation(qu) @staticmethod - def from_Eulers(eulers, + def from_Eulers(phi, degrees = False): + """ + Initialize from Bunge-Euler angles. - eu = np.array(eulers,dtype=float) + Parameters + ---------- + phi: numpy.ndarray of shape (...,3) + Bunge-Euler angles: (φ_1, ϕ, φ_2), φ_1 ∈ [0,2π], ϕ ∈ [0,π], φ_2 ∈ [0,2π] + unless degrees == True: φ_1 ∈ [0,360], ϕ ∈ [0,180], φ_2 ∈ [0,360]. + degrees: boolean, optional + Bunge-Euler angles are given in degrees. Defaults to False. + + """ + eu = np.array(phi,dtype=float) if eu.shape[:-2:-1] != (3,): raise ValueError('Invalid shape.') @@ -316,12 +384,29 @@ class Rotation: degrees = False, normalise = False, P = -1): + """ + Initialize from Axis angle pair. + Parameters + ---------- + axis_angle: numpy.ndarray of shape (...,4) + Axis angle pair: [n_1, n_2, n_3, ω], |n| = 1 and ω ∈ [0,π] + unless degrees = True: ω ∈ [0,180]. + degrees: boolean, optional + Angle ω is given in degrees. Defaults to False. + normalize: boolean, optional + Allow |n| ≠ 1. Defaults to False. + P: integer ∈ {-1,1}, optional + Convention used. Defaults to -1. + + """ ax = np.array(axis_angle,dtype=float) if ax.shape[:-2:-1] != (4,): raise ValueError('Invalid shape.') + if abs(P) != 1: + raise ValueError('P ∉ {-1,1}') - if P > 0: ax[...,0:3] *= -1 # convert from P=1 to P=-1 + if P == 1: ax[...,0:3] *= -1 if degrees: ax[..., 3] = np.radians(ax[...,3]) if normalise: ax[...,0:3] /= np.linalg.norm(ax[...,0:3],axis=-1) if np.any(ax[...,3] < 0.0) or np.any(ax[...,3] > np.pi): @@ -335,7 +420,19 @@ class Rotation: def from_basis(basis, orthonormal = True, reciprocal = False): + """ + Initialize from tbd. + Parameters + ---------- + basis: numpy.ndarray of shape (...,3,3) + tbd + orthonormal: boolean, optional + tbd. Defaults to True. + reciprocal: boolean, optional + tbd. Defaults to False. + + """ om = np.array(basis,dtype=float) if om.shape[:-3:-1] != (3,3): raise ValueError('Invalid shape.') @@ -356,20 +453,64 @@ class Rotation: return Rotation(Rotation._om2qu(om)) @staticmethod - def from_matrix(om): + def from_directions(hkl,uvw): + """ + Initialize from pair of directions/planes. - return Rotation.from_basis(om) + Parameters + ---------- + hkl: numpy.ndarray of shape (...,3) + Direction parallel to z direction, i.e. (h k l) || (0,0,1). + uvw: numpy.ndarray of shape (...,3) + Direction parallel to x direction, i.e. || (1,0,0). + + """ + hkl_ = hkl/np.linalg.norm(hkl,axis=-1,keepdims=True) + uvw_ = uvw/np.linalg.norm(uvw,axis=-1,keepdims=True) + v_1 = np.block([uvw_,np.cross(hkl_,uvw_),hkl_]).reshape(hkl_.shape+(3,)) + v_2 = np.block([uvw_,np.cross(uvw_,hkl_),hkl_]).reshape(hkl_.shape+(3,)) + R = np.where(np.broadcast_to(np.expand_dims(np.expand_dims(np.linalg.det(v_1)>0,-1),-1),v_1.shape), + v_1,v_2) + return Rotation.from_basis(np.swapaxes(R,axis2=-2,axis1=-1)) @staticmethod - def from_Rodrigues(rodrigues, + def from_matrix(R): + """ + Initialize from rotation matrix. + + Parameters + ---------- + R: numpy.ndarray of shape (...,3,3) + Rotation matrix: det(R) = 1, R.T∙R=I. + + """ + return Rotation.from_basis(R) + + @staticmethod + def from_Rodrigues(rho, normalise = False, P = -1): + """ + Initialize from Rodrigues-Frank vector. - ro = np.array(rodrigues,dtype=float) + Parameters + ---------- + rho: numpy.ndarray of shape (...,4) + Rodrigues-Frank vector (angle separated from axis). + (n_1, n_2, n_3, tan(ω/2)), |n| = 1 and ω ∈ [0,π]. + normalize: boolean, optional + Allow |n| ≠ 1. Defaults to False. + P: integer ∈ {-1,1}, optional + Convention used. Defaults to -1. + + """ + ro = np.array(rho,dtype=float) if ro.shape[:-2:-1] != (4,): raise ValueError('Invalid shape.') + if abs(P) != 1: + raise ValueError('P ∉ {-1,1}') - if P > 0: ro[...,0:3] *= -1 # convert from P=1 to P=-1 + if P == 1: ro[...,0:3] *= -1 if normalise: ro[...,0:3] /= np.linalg.norm(ro[...,0:3],axis=-1) if np.any(ro[...,3] < 0.0): raise ValueError('Rodrigues vector rotation angle not positive.') @@ -379,14 +520,26 @@ class Rotation: return Rotation(Rotation._ro2qu(ro)) @staticmethod - def from_homochoric(homochoric, + def from_homochoric(h, P = -1): + """ + Initialize from homochoric vector. - ho = np.array(homochoric,dtype=float) + Parameters + ---------- + h: numpy.ndarray of shape (...,3) + Homochoric vector: (h_1, h_2, h_3), |h| < (3/4*π)^(1/3). + P: integer ∈ {-1,1}, optional + Convention used. Defaults to -1. + + """ + ho = np.array(h,dtype=float) if ho.shape[:-2:-1] != (3,): raise ValueError('Invalid shape.') + if abs(P) != 1: + raise ValueError('P ∉ {-1,1}') - if P > 0: ho *= -1 # convert from P=1 to P=-1 + if P == 1: ho *= -1 if np.any(np.linalg.norm(ho,axis=-1) >_R1+1e-9): raise ValueError('Homochoric coordinate outside of the sphere.') @@ -394,18 +547,30 @@ class Rotation: return Rotation(Rotation._ho2qu(ho)) @staticmethod - def from_cubochoric(cubochoric, - P = -1): + def from_cubochoric(c, + P = -1): + """ + Initialize from cubochoric vector. - cu = np.array(cubochoric,dtype=float) + Parameters + ---------- + c: numpy.ndarray of shape (...,3) + Cubochoric vector: (c_1, c_2, c_3), max(c_i) < 1/2*π^(2/3). + P: integer ∈ {-1,1}, optional + Convention used. Defaults to -1. + + """ + cu = np.array(c,dtype=float) if cu.shape[:-2:-1] != (3,): raise ValueError('Invalid shape.') + if abs(P) != 1: + raise ValueError('P ∉ {-1,1}') - if np.abs(np.max(cu))>np.pi**(2./3.) * 0.5+1e-9: - raise ValueError('Cubochoric coordinate outside of the cube: {} {} {}.'.format(*cu)) + if np.abs(np.max(cu)) > np.pi**(2./3.) * 0.5+1e-9: + raise ValueError('Cubochoric coordinate outside of the cube.') ho = Rotation._cu2ho(cu) - if P > 0: ho *= -1 # convert from P=1 to P=-1 + if P == 1: ho *= -1 return Rotation(Rotation._ho2qu(ho)) @@ -458,7 +623,7 @@ class Rotation: np.cos(2.0*np.pi*r[...,1])*B, np.sin(2.0*np.pi*r[...,0])*A],axis=-1) - return Rotation(q.reshape(r.shape[:-1]+(4,)) if shape is not None else q).standardize() + return Rotation(q.reshape(r.shape[:-1]+(4,)) if shape is not None else q)._standardize() # for compatibility (old names do not follow convention) diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 5a1cd145d..a6b4fea69 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -14,7 +14,7 @@ scatter=1.e-2 @pytest.fixture def default(): - """A set of n random rotations.""" + """A set of n rotations (corner cases and random).""" specials = np.array([ [1.0, 0.0, 0.0, 0.0], #---------------------- @@ -170,7 +170,7 @@ def qu2ax(qu): Modified version of the original formulation, should be numerically more stable """ - if np.isclose(qu[0],1.,rtol=0.0): # set axis to [001] if the angle is 0/360 + if np.isclose(qu[0],1.,rtol=0.0): # set axis to [001] if the angle is 0/360 ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) elif qu[0] > 1.e-8: s = np.sign(qu[0])/np.sqrt(qu[1]**2+qu[2]**2+qu[3]**2) @@ -534,7 +534,7 @@ def mul(me, other): other_p = other.quaternion[1:] R = me.__class__(np.append(me_q*other_q - np.dot(me_p,other_p), me_q*other_p + other_q*me_p + _P * np.cross(me_p,other_p))) - return R.standardize() + return R._standardize() elif isinstance(other, np.ndarray): if other.shape == (3,): A = me.quaternion[0]**2.0 - np.dot(me.quaternion[1:],me.quaternion[1:]) @@ -554,7 +554,7 @@ def mul(me, other): axes = ((0, 2, 4, 6), (0, 1, 2, 3)) return np.tensordot(RRRR, other, axes) else: - raise ValueError('Can only rotate vectors, 2nd order ternsors, and 4th order tensors') + raise ValueError('Can only rotate vectors, 2nd order tensors, and 4th order tensors') else: raise TypeError('Cannot rotate {}'.format(type(other))) @@ -854,6 +854,15 @@ class TestRotation: rot = Rotation.from_basis(om,False,reciprocal=reciprocal) assert np.isclose(np.linalg.det(rot.as_matrix()),1.0) + def test_directions(self): + hkl = np.array([0.,0.,1.]) + uvw = np.array([1.,0.,0.]) + assert np.allclose(Rotation.from_directions(hkl,uvw).as_matrix(),np.eye(3)) + + @pytest.mark.parametrize('shape',[None,1,(4,4)]) + def test_random(self,shape): + Rotation.from_random(shape) + @pytest.mark.parametrize('function',[Rotation.from_quaternion, Rotation.from_Eulers, Rotation.from_axis_angle, @@ -866,6 +875,16 @@ class TestRotation: with pytest.raises(ValueError): function(invalid_shape) + @pytest.mark.parametrize('fr,to',[(Rotation.from_quaternion,'as_quaternion'), + (Rotation.from_axis_angle,'as_axis_angle'), + (Rotation.from_Rodrigues, 'as_Rodrigues'), + (Rotation.from_homochoric,'as_homochoric'), + (Rotation.from_cubochoric,'as_cubochoric')]) + def test_invalid_P(self,fr,to): + R = Rotation.from_random(np.random.randint(8,32,(3))) # noqa + with pytest.raises(ValueError): + fr(eval('R.{}()'.format(to)),P=-30) + @pytest.mark.parametrize('shape',[None,(3,),(4,2)]) def test_broadcast(self,shape): rot = Rotation.from_random(shape) @@ -932,14 +951,18 @@ class TestRotation: phi_2 = 2*np.pi - phi_1 R_1 = Rotation.from_Eulers(np.array([phi_1,0.,0.])) R_2 = Rotation.from_Eulers(np.array([0.,0.,phi_2])) - assert np.allclose(data,R_2*(R_1*data)) + assert np.allclose(data,R_2@(R_1@data)) + + def test_rotate_inverse(self): + R = Rotation.from_random() + assert np.allclose(np.eye(3),(R.inversed()@R).as_matrix()) @pytest.mark.parametrize('data',[np.random.rand(3), np.random.rand(3,3), np.random.rand(3,3,3,3)]) - def test_rotate_inverse(self,data): + def test_rotate_inverse_array(self,data): R = Rotation.from_random() - assert np.allclose(data,R.inversed()*(R*data)) + assert np.allclose(data,R.inversed()@(R@data)) @pytest.mark.parametrize('data',[np.random.rand(4), np.random.rand(3,2), @@ -956,3 +979,7 @@ class TestRotation: R = Rotation.from_random() with pytest.raises(TypeError): R*data + + def test_misorientation(self): + R = Rotation.from_random() + assert np.allclose(R.misorientation(R).as_matrix(),np.eye(3)) From 3f63a4fdbc56f96dae3f87952b91c9c7d8a667db Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 20 Jun 2020 18:13:34 +0200 Subject: [PATCH 25/43] [skip ci] typo --- python/damask/_rotation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index a33b4bc51..09aeab256 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -289,7 +289,7 @@ class Rotation: def as_cubochoric(self): """ - Represent as cubocoric vector. + Represent as cubochoric vector. Returns ------- From 15b43bcebf9ae932d54e5a623b82177351214569 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 20 Jun 2020 19:57:49 +0200 Subject: [PATCH 26/43] from_directions is not general, removed polishing --- python/damask/_rotation.py | 67 ++++++++++++----------------------- python/tests/test_Rotation.py | 5 --- 2 files changed, 23 insertions(+), 49 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 09aeab256..7c469a2b5 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -325,13 +325,13 @@ class Rotation: Parameters ---------- - q: numpy.ndarray of shape (...,4) + q : numpy.ndarray of shape (...,4) Unit quaternion in positive real hemisphere: (q_0, q_1, q_2, q_3), |q|=1, q_0 ≥ 0. - accept_homomorph: boolean, optional + accept_homomorph : boolean, optional Allow homomorphic variants, i.e. q_0 < 0 (negative real hemisphere). Defaults to False. - P: integer ∈ {-1,1}, optional + P : integer ∈ {-1,1}, optional Convention used. Defaults to -1. """ @@ -362,10 +362,10 @@ class Rotation: Parameters ---------- - phi: numpy.ndarray of shape (...,3) + phi : numpy.ndarray of shape (...,3) Bunge-Euler angles: (φ_1, ϕ, φ_2), φ_1 ∈ [0,2π], ϕ ∈ [0,π], φ_2 ∈ [0,2π] unless degrees == True: φ_1 ∈ [0,360], ϕ ∈ [0,180], φ_2 ∈ [0,360]. - degrees: boolean, optional + degrees : boolean, optional Bunge-Euler angles are given in degrees. Defaults to False. """ @@ -389,14 +389,14 @@ class Rotation: Parameters ---------- - axis_angle: numpy.ndarray of shape (...,4) + axis_angle : numpy.ndarray of shape (...,4) Axis angle pair: [n_1, n_2, n_3, ω], |n| = 1 and ω ∈ [0,π] unless degrees = True: ω ∈ [0,180]. - degrees: boolean, optional + degrees : boolean, optional Angle ω is given in degrees. Defaults to False. normalize: boolean, optional Allow |n| ≠ 1. Defaults to False. - P: integer ∈ {-1,1}, optional + P : integer ∈ {-1,1}, optional Convention used. Defaults to -1. """ @@ -421,16 +421,16 @@ class Rotation: orthonormal = True, reciprocal = False): """ - Initialize from tbd. + Initialize from lattice basis vectors. Parameters ---------- - basis: numpy.ndarray of shape (...,3,3) - tbd - orthonormal: boolean, optional - tbd. Defaults to True. - reciprocal: boolean, optional - tbd. Defaults to False. + basis : numpy.ndarray of shape (...,3,3) + Three lattice basis vectors in three dimensions. + orthonormal : boolean, optional + Basis is strictly orthonormal, i.e. is free of stretch components. Defaults to True. + reciprocal : boolean, optional + Basis vectors are given in reciprocal (instead of real) space. Defaults to False. """ om = np.array(basis,dtype=float) @@ -452,27 +452,6 @@ class Rotation: return Rotation(Rotation._om2qu(om)) - @staticmethod - def from_directions(hkl,uvw): - """ - Initialize from pair of directions/planes. - - Parameters - ---------- - hkl: numpy.ndarray of shape (...,3) - Direction parallel to z direction, i.e. (h k l) || (0,0,1). - uvw: numpy.ndarray of shape (...,3) - Direction parallel to x direction, i.e. || (1,0,0). - - """ - hkl_ = hkl/np.linalg.norm(hkl,axis=-1,keepdims=True) - uvw_ = uvw/np.linalg.norm(uvw,axis=-1,keepdims=True) - v_1 = np.block([uvw_,np.cross(hkl_,uvw_),hkl_]).reshape(hkl_.shape+(3,)) - v_2 = np.block([uvw_,np.cross(uvw_,hkl_),hkl_]).reshape(hkl_.shape+(3,)) - R = np.where(np.broadcast_to(np.expand_dims(np.expand_dims(np.linalg.det(v_1)>0,-1),-1),v_1.shape), - v_1,v_2) - return Rotation.from_basis(np.swapaxes(R,axis2=-2,axis1=-1)) - @staticmethod def from_matrix(R): """ @@ -480,7 +459,7 @@ class Rotation: Parameters ---------- - R: numpy.ndarray of shape (...,3,3) + R : numpy.ndarray of shape (...,3,3) Rotation matrix: det(R) = 1, R.T∙R=I. """ @@ -495,12 +474,12 @@ class Rotation: Parameters ---------- - rho: numpy.ndarray of shape (...,4) + rho : numpy.ndarray of shape (...,4) Rodrigues-Frank vector (angle separated from axis). (n_1, n_2, n_3, tan(ω/2)), |n| = 1 and ω ∈ [0,π]. - normalize: boolean, optional + normalize : boolean, optional Allow |n| ≠ 1. Defaults to False. - P: integer ∈ {-1,1}, optional + P : integer ∈ {-1,1}, optional Convention used. Defaults to -1. """ @@ -527,9 +506,9 @@ class Rotation: Parameters ---------- - h: numpy.ndarray of shape (...,3) + h : numpy.ndarray of shape (...,3) Homochoric vector: (h_1, h_2, h_3), |h| < (3/4*π)^(1/3). - P: integer ∈ {-1,1}, optional + P : integer ∈ {-1,1}, optional Convention used. Defaults to -1. """ @@ -554,9 +533,9 @@ class Rotation: Parameters ---------- - c: numpy.ndarray of shape (...,3) + c : numpy.ndarray of shape (...,3) Cubochoric vector: (c_1, c_2, c_3), max(c_i) < 1/2*π^(2/3). - P: integer ∈ {-1,1}, optional + P : integer ∈ {-1,1}, optional Convention used. Defaults to -1. """ diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index a6b4fea69..933719d6f 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -854,11 +854,6 @@ class TestRotation: rot = Rotation.from_basis(om,False,reciprocal=reciprocal) assert np.isclose(np.linalg.det(rot.as_matrix()),1.0) - def test_directions(self): - hkl = np.array([0.,0.,1.]) - uvw = np.array([1.,0.,0.]) - assert np.allclose(Rotation.from_directions(hkl,uvw).as_matrix(),np.eye(3)) - @pytest.mark.parametrize('shape',[None,1,(4,4)]) def test_random(self,shape): Rotation.from_random(shape) From 4c5939ef233ff4021d0f836f836df3d21ae25c91 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 20 Jun 2020 22:51:00 +0200 Subject: [PATCH 27/43] small polishing --- python/damask/util.py | 14 +++++++------- python/tests/test_Rotation.py | 4 ++-- python/tests/test_util.py | 11 +++++++++++ src/discretization.f90 | 2 +- 4 files changed, 21 insertions(+), 10 deletions(-) create mode 100644 python/tests/test_util.py diff --git a/python/damask/util.py b/python/damask/util.py index cb1d8757d..a3fc71d3a 100644 --- a/python/damask/util.py +++ b/python/damask/util.py @@ -93,7 +93,7 @@ def strikeout(what): def execute(cmd, - streamIn = None, + stream_in = None, wd = './', env = None): """ @@ -103,7 +103,7 @@ def execute(cmd, ---------- cmd : str Command to be executed. - streanIn :, optional + stream_in : file object, optional Input (via pipe) for executed process. wd : str, optional Working directory of process. Defaults to ./ . @@ -119,14 +119,14 @@ def execute(cmd, stderr = subprocess.PIPE, stdin = subprocess.PIPE, env = myEnv) - out,error = [i for i in (process.communicate() if streamIn is None - else process.communicate(streamIn.read().encode('utf-8')))] + stdout, stderr = [i for i in (process.communicate() if stream_in is None + else process.communicate(stream_in.read().encode('utf-8')))] os.chdir(initialPath) - out = out.decode('utf-8').replace('\x08','') - error = error.decode('utf-8').replace('\x08','') + stdout = stdout.decode('utf-8').replace('\x08','') + stderr = stderr.decode('utf-8').replace('\x08','') if process.returncode != 0: raise RuntimeError('{} failed with returncode {}'.format(cmd,process.returncode)) - return out,error + return stdout, stderr def show_progress(iterable,N_iter=None,prefix='',bar_length=50): diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 933719d6f..86e5fa2a7 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -950,14 +950,14 @@ class TestRotation: def test_rotate_inverse(self): R = Rotation.from_random() - assert np.allclose(np.eye(3),(R.inversed()@R).as_matrix()) + assert np.allclose(np.eye(3),(~R@R).as_matrix()) @pytest.mark.parametrize('data',[np.random.rand(3), np.random.rand(3,3), np.random.rand(3,3,3,3)]) def test_rotate_inverse_array(self,data): R = Rotation.from_random() - assert np.allclose(data,R.inversed()@(R@data)) + assert np.allclose(data,~R@(R@data)) @pytest.mark.parametrize('data',[np.random.rand(4), np.random.rand(3,2), diff --git a/python/tests/test_util.py b/python/tests/test_util.py new file mode 100644 index 000000000..e67478d82 --- /dev/null +++ b/python/tests/test_util.py @@ -0,0 +1,11 @@ +from damask import util + +class TestUtil: + + def test_execute_direct(self): + out,err = util.execute('echo test') + assert out=='test\n' and err=='' + + def test_execute_env(self): + out,err = util.execute('sh -c "echo $test_for_execute"',env={'test_for_execute':'test'}) + assert out=='test\n' and err=='' diff --git a/src/discretization.f90 b/src/discretization.f90 index 16f76b158..e52784c6d 100644 --- a/src/discretization.f90 +++ b/src/discretization.f90 @@ -49,7 +49,7 @@ subroutine discretization_init(homogenizationAt,microstructureAt,& IPcoords0, & NodeCoords0 integer, optional, intent(in) :: & - sharedNodesBegin ! index of first node shared among different processes (MPI) + sharedNodesBegin !< index of first node shared among different processes (MPI) write(6,'(/,a)') ' <<<+- discretization init -+>>>'; flush(6) From d4efadb333f864b4ebef645373cdbeccb2356e83 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 21 Jun 2020 10:03:52 +0200 Subject: [PATCH 28/43] should be availabe outside of this module --- src/quaternions.f90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/quaternions.f90 b/src/quaternions.f90 index f5f5276e8..c337bd681 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -101,6 +101,8 @@ module quaternions assignment(=), & conjg, aimag, & log, exp, & + abs, dot_product, & + inverse, & real contains From c6a5bb8a3b873a66ba968339d76c754bc63c0d4d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 21 Jun 2020 10:04:45 +0200 Subject: [PATCH 29/43] is 2020 --- python/damask/_rotation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 7c469a2b5..f433b2b37 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -615,7 +615,7 @@ class Rotation: #################################################################################################### # Code below available according to the following conditions on https://github.com/MarDiehl/3Drotations #################################################################################################### -# Copyright (c) 2017-2019, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH +# Copyright (c) 2017-2020, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH # Copyright (c) 2013-2014, Marc De Graef/Carnegie Mellon University # All rights reserved. # From deaf055f8876a55c6bdfee0e0089cb8bfefbadb0 Mon Sep 17 00:00:00 2001 From: Test User Date: Mon, 22 Jun 2020 15:59:31 +0200 Subject: [PATCH 30/43] [skip ci] updated version information after successful test of v2.0.3-2680-g17212b27 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index b25957f75..d041f794f 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-2664-ge959aaab +v2.0.3-2680-g17212b27 From e2a0e98267b90db4e2d71eb04860f1be82646c8d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 23 Jun 2020 11:06:41 +0200 Subject: [PATCH 31/43] small polishing --- src/CPFEM.f90 | 1 - src/YAML_parse.f90 | 2 +- src/YAML_types.f90 | 4 ++-- src/crystallite.f90 | 2 +- src/grid/spectral_utilities.f90 | 2 +- src/numerics.f90 | 8 +++----- 6 files changed, 8 insertions(+), 11 deletions(-) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index fce375259..58beb72a0 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -140,7 +140,6 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS num_commercialFEM => numerics_root%get('commercialFEM',defaultVal=emptyDict) iJacoStiffness = num_commercialFEM%get_asInt('ijacostiffness',defaultVal=1) if (iJacoStiffness < 1) call IO_error(301,ext_msg='iJacoStiffness') -!------------------------------------------------------------------------------ elCP = mesh_FEM2DAMASK_elem(elFE) diff --git a/src/YAML_parse.f90 b/src/YAML_parse.f90 index 2b56ba946..8850a62fc 100644 --- a/src/YAML_parse.f90 +++ b/src/YAML_parse.f90 @@ -51,7 +51,7 @@ recursive function parse_flow(flow_string,defaultVal) result(node) if (len_trim(flow_string) == 0 .and. present(defaultVal)) then node => defaultVal return - elseif (flow_string(1:1) == '{') then ! start of a dictionary + elseif (flow_string(1:1) == '{') then ! start of a dictionary e = 1 allocate(tDict::node) do while (e < len_trim(flow_string)) diff --git a/src/YAML_types.f90 b/src/YAML_types.f90 index 6a2706269..ad3db9d47 100644 --- a/src/YAML_types.f90 +++ b/src/YAML_types.f90 @@ -155,9 +155,9 @@ module YAML_types final :: tItem_finalize end type tItem - type(tDict), target,public :: & + type(tDict), target, public :: & emptyDict - type(tList), target,public :: & + type(tList), target, public :: & emptyList abstract interface diff --git a/src/crystallite.f90 b/src/crystallite.f90 index eb10e8550..9113fa12f 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -129,7 +129,7 @@ subroutine crystallite_init eMax, & !< maximum number of elements myNcomponents !< number of components at current IP - class(tNode) , pointer :: & + class(tNode), pointer :: & num_crystallite write(6,'(/,a)') ' <<<+- crystallite init -+>>>' diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index f09759009..cc9d3dd65 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -191,7 +191,7 @@ subroutine spectral_utilities_init tensorSize = 9_C_INTPTR_T character(len=pStringLen) :: & petsc_options - class (tNode) , pointer :: & + class(tNode), pointer :: & num_grid write(6,'(/,a)') ' <<<+- spectral_utilities init -+>>>' diff --git a/src/numerics.f90 b/src/numerics.f90 index b4d029d00..bb861a6d0 100644 --- a/src/numerics.f90 +++ b/src/numerics.f90 @@ -1,6 +1,7 @@ !-------------------------------------------------------------------------------------------------- !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH +!> @author Sharan Roongta, Max-Planck-Institut für Eisenforschung GmbH !> @brief Managing of parameters related to numerics !-------------------------------------------------------------------------------------------------- module numerics @@ -64,16 +65,13 @@ subroutine numerics_init numerics_root => emptyDict inquire(file='numerics.yaml', exist=fexist) - fileExists: if (fexist) then + if (fexist) then write(6,'(a,/)') ' using values from config file' flush(6) numerics_input = IO_read('numerics.yaml') numerics_inFlow = to_flow(numerics_input) numerics_root => parse_flow(numerics_inFlow,defaultVal=emptyDict) - else fileExists - write(6,'(a,/)') ' using standard values' - flush(6) - endif fileExists + endif !-------------------------------------------------------------------------------------------------- ! openMP parameter From fdf7887b47a6e8f6b30942856a529368d9ba0bd8 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Wed, 24 Jun 2020 11:56:21 +0200 Subject: [PATCH 32/43] handle default internally --- src/YAML_parse.f90 | 7 +++---- src/numerics.f90 | 2 +- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/src/YAML_parse.f90 b/src/YAML_parse.f90 index 8850a62fc..cfbcf11f3 100644 --- a/src/YAML_parse.f90 +++ b/src/YAML_parse.f90 @@ -34,10 +34,9 @@ end subroutine YAML_init !> @brief reads the flow style string and stores it in the form of dictionaries, lists and scalars. !> @details A node type pointer can either point to a dictionary, list or scalar type entities. !-------------------------------------------------------------------------------------------------- -recursive function parse_flow(flow_string,defaultVal) result(node) +recursive function parse_flow(flow_string) result(node) character(len=*), intent(inout) :: flow_string - class(tDict), intent(in), optional, target :: defaultVal class (tNode), pointer :: node class (tNode), pointer :: myVal @@ -48,8 +47,8 @@ recursive function parse_flow(flow_string,defaultVal) result(node) d !> position of key: value separator (':') flow_string = trim(adjustl(flow_string(:))) - if (len_trim(flow_string) == 0 .and. present(defaultVal)) then - node => defaultVal + if (len_trim(flow_string) == 0) then + node => emptyDict return elseif (flow_string(1:1) == '{') then ! start of a dictionary e = 1 diff --git a/src/numerics.f90 b/src/numerics.f90 index bb861a6d0..132c01f87 100644 --- a/src/numerics.f90 +++ b/src/numerics.f90 @@ -70,7 +70,7 @@ subroutine numerics_init flush(6) numerics_input = IO_read('numerics.yaml') numerics_inFlow = to_flow(numerics_input) - numerics_root => parse_flow(numerics_inFlow,defaultVal=emptyDict) + numerics_root => parse_flow(numerics_inFlow) endif !-------------------------------------------------------------------------------------------------- From 692fc98fd5605c265017e5d8b0bb99b27fbc84bf Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Wed, 24 Jun 2020 16:31:45 +0200 Subject: [PATCH 33/43] 'num' structure for data to avoid multiple reading of parameters --- src/CPFEM.f90 | 31 +++++++++++++++++++------------ src/damage_local.f90 | 27 ++++++++++++++++++--------- src/damage_nonlocal.f90 | 26 +++++++++++++++++--------- src/grid/DAMASK_grid.f90 | 21 +++++++++++++-------- src/math.f90 | 17 ++++++++++++----- 5 files changed, 79 insertions(+), 43 deletions(-) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 58beb72a0..6647581e3 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -5,7 +5,6 @@ !-------------------------------------------------------------------------------------------------- module CPFEM use prec - use numerics use debug use FEsolving use math @@ -45,6 +44,13 @@ module CPFEM CPFEM_BACKUPJACOBIAN = 2_pInt**2_pInt, & CPFEM_RESTOREJACOBIAN = 2_pInt**3_pInt + type, private :: tNumerics + integer :: & + iJacoStiffness !< frequency of stiffness update + end type tNumerics + + type(tNumerics), private :: num + public :: & CPFEM_general, & CPFEM_initAll, & @@ -86,6 +92,9 @@ end subroutine CPFEM_initAll !-------------------------------------------------------------------------------------------------- subroutine CPFEM_init + class(tNode), pointer :: & + num_commercialFEM + write(6,'(/,a)') ' <<<+- CPFEM init -+>>>' flush(6) @@ -93,6 +102,13 @@ subroutine CPFEM_init allocate(CPFEM_dcsdE( 6,6,discretization_nIP,discretization_nElem), source= 0.0_pReal) allocate(CPFEM_dcsdE_knownGood(6,6,discretization_nIP,discretization_nElem), source= 0.0_pReal) +!------------------------------------------------------------------------------ +! read numerical parameters and do sanity check + num_commercialFEM => numerics_root%get('commercialFEM',defaultVal=emptyDict) + num%iJacoStiffness = num_commercialFEM%get_asInt('ijacostiffness',defaultVal=1) + if (num%iJacoStiffness < 1) call IO_error(301,ext_msg='iJacoStiffness') +!------------------------------------------------------------------------------ + if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0) then write(6,'(a32,1x,6(i8,1x))') 'CPFEM_cs: ', shape(CPFEM_cs) write(6,'(a32,1x,6(i8,1x))') 'CPFEM_dcsdE: ', shape(CPFEM_dcsdE) @@ -125,21 +141,12 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS H integer(pInt) elCP, & ! crystal plasticity element number - i, j, k, l, m, n, ph, homog, mySource, & - iJacoStiffness !< frequency of stiffness update + i, j, k, l, m, n, ph, homog, mySource logical updateJaco ! flag indicating if Jacobian has to be updated real(pReal), parameter :: ODD_STRESS = 1e15_pReal, & !< return value for stress if terminallyIll ODD_JACOBIAN = 1e50_pReal !< return value for jacobian if terminallyIll - class(tNode), pointer :: & - num_commercialFEM - -!------------------------------------------------------------------------------ -! read numerical parameters and do sanity check - num_commercialFEM => numerics_root%get('commercialFEM',defaultVal=emptyDict) - iJacoStiffness = num_commercialFEM%get_asInt('ijacostiffness',defaultVal=1) - if (iJacoStiffness < 1) call IO_error(301,ext_msg='iJacoStiffness') elCP = mesh_FEM2DAMASK_elem(elFE) @@ -179,7 +186,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_identity2nd(6) else validCalculation - updateJaco = mod(cycleCounter,iJacoStiffness) == 0 + updateJaco = mod(cycleCounter,num%iJacoStiffness) == 0 FEsolving_execElem = elCP FEsolving_execIP = ip if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) & diff --git a/src/damage_local.f90 b/src/damage_local.f90 index 892768fdc..e0dfa34eb 100644 --- a/src/damage_local.f90 +++ b/src/damage_local.f90 @@ -23,9 +23,16 @@ module damage_local output end type tParameters + type, private :: tNumerics + real(pReal) :: & + residualStiffness !< non-zero residual damage + end type tNumerics + type(tparameters), dimension(:), allocatable :: & param + type(tNumerics), private :: num + public :: & damage_local_init, & damage_local_updateState, & @@ -40,9 +47,17 @@ contains subroutine damage_local_init integer :: Ninstance,NofMyHomog,h + class(tNode), pointer :: num_generic write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_local_label//' init -+>>>'; flush(6) +!---------------------------------------------------------------------------------------------- +! read numerics parameter and do sanity check + num_generic => numerics_root%get('generic',defaultVal=emptyDict) + num%residualStiffness = num_generic%get_asFloat('residualStiffness', defaultVal=1.0e-6_pReal) + if (num%residualStiffness < 0.0_pReal) call IO_error(301,ext_msg='residualStiffness') +!---------------------------------------------------------------------------------------------- + Ninstance = count(damage_type == DAMAGE_local_ID) allocate(param(Ninstance)) @@ -85,20 +100,14 @@ function damage_local_updateState(subdt, ip, el) homog, & offset real(pReal) :: & - phi, phiDot, dPhiDot_dPhi, & - residualStiffness !< non-zero residual damage - class(tNode), pointer :: & - num_generic + phi, phiDot, dPhiDot_dPhi - num_generic => numerics_root%get('generic',defaultVal=emptyDict) - residualStiffness = num_generic%get_asFloat('residualStiffness', defaultVal=1.0e-6_pReal) - if (residualStiffness < 0.0_pReal) call IO_error(301,ext_msg='residualStiffness') - + homog = material_homogenizationAt(el) offset = material_homogenizationMemberAt(ip,el) phi = damageState(homog)%subState0(1,offset) call damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el) - phi = max(residualStiffness,min(1.0_pReal,phi + subdt*phiDot)) + phi = max(num%residualStiffness,min(1.0_pReal,phi + subdt*phiDot)) damage_local_updateState = [ abs(phi - damageState(homog)%state(1,offset)) & <= 1.0e-2_pReal & diff --git a/src/damage_nonlocal.f90 b/src/damage_nonlocal.f90 index f2f26889d..466261225 100644 --- a/src/damage_nonlocal.f90 +++ b/src/damage_nonlocal.f90 @@ -24,8 +24,15 @@ module damage_nonlocal output end type tParameters + type, private :: tNumerics + real(pReal) :: & + charLength !< characteristic length scale for gradient problems + end type tNumerics + type(tparameters), dimension(:), allocatable :: & param + type(tNumerics), private :: & + num public :: & damage_nonlocal_init, & @@ -44,9 +51,17 @@ contains subroutine damage_nonlocal_init integer :: Ninstance,NofMyHomog,h - + class(tNode), pointer :: & + num_generic + write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_nonlocal_label//' init -+>>>'; flush(6) +!------------------------------------------------------------------------------------ +! read numerics parameter + num_generic => numerics_root%get('generic',defaultVal= emptyDict) + num%charLength = num_generic%get_asFloat('charLength',defaultVal=1.0_pReal) +!------------------------------------------------------------------------------------ + Ninstance = count(damage_type == DAMAGE_nonlocal_ID) allocate(param(Ninstance)) @@ -139,13 +154,6 @@ function damage_nonlocal_getDiffusion(ip,el) integer :: & homog, & grain - real(pReal) :: & - charLength !< characteristic length scale for gradient problems - class(tNode), pointer :: & - num_generic - - num_generic => numerics_root%get('generic',defaultVal= emptyDict) - charLength = num_generic%get_asFloat('charLength',defaultVal=1.0_pReal) homog = material_homogenizationAt(el) damage_nonlocal_getDiffusion = 0.0_pReal @@ -155,7 +163,7 @@ function damage_nonlocal_getDiffusion(ip,el) enddo damage_nonlocal_getDiffusion = & - charLength**2*damage_nonlocal_getDiffusion/real(homogenization_Ngrains(homog),pReal) + num%charLength**2*damage_nonlocal_getDiffusion/real(homogenization_Ngrains(homog),pReal) end function damage_nonlocal_getDiffusion diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index e151222c2..bc153693d 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -61,7 +61,6 @@ program DAMASK_grid i, j, k, l, field, & errorID = 0, & cutBackLevel = 0, & !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$ - maxCutBack, & !< max number of cut backs stepFraction = 0 !< fraction of current time interval integer :: & currentLoadcase = 0, & !< current load case @@ -69,12 +68,18 @@ program DAMASK_grid totalIncsCounter = 0, & !< total # of increments statUnit = 0, & !< file unit for statistics output stagIter, & - stagItMax, & !< max number of field level staggered iterations nActiveFields = 0 character(len=pStringLen), dimension(:), allocatable :: fileContent character(len=pStringLen) :: & incInfo, & loadcase_string + type :: tNumerics + integer :: & + maxCutBack, & !< max number of cut backs + stagItMax !< max number of field level staggered iterations + end type tNumerics + + type(tNumerics) :: num type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases type(tLoadCase) :: newLoadCase type(tSolutionState), allocatable, dimension(:) :: solres @@ -114,11 +119,11 @@ program DAMASK_grid !------------------------------------------------------------------------------------------------- ! reading field paramters from numerics file and do sanity checks num_grid => numerics_root%get('grid', defaultVal=emptyDict) - stagItMax = num_grid%get_asInt('maxStaggeredIter',defaultVal=10) - maxCutBack = num_grid%get_asInt('maxCutBack',defaultVal=3) + num%stagItMax = num_grid%get_asInt('maxStaggeredIter',defaultVal=10) + num%maxCutBack = num_grid%get_asInt('maxCutBack',defaultVal=3) - if (stagItMax < 0) call IO_error(301,ext_msg='maxStaggeredIter') - if (maxCutBack < 0) call IO_error(301,ext_msg='maxCutBack') + if (num%stagItMax < 0) call IO_error(301,ext_msg='maxStaggeredIter') + if (num%maxCutBack < 0) call IO_error(301,ext_msg='maxCutBack') !-------------------------------------------------------------------------------------------------- ! assign mechanics solver depending on selected type @@ -449,7 +454,7 @@ program DAMASK_grid enddo stagIter = stagIter + 1 - stagIterate = stagIter < stagItMax & + stagIterate = stagIter < num%stagItMax & .and. all(solres(:)%converged) & .and. .not. all(solres(:)%stagConverged) ! stationary with respect to staggered iteration enddo @@ -468,7 +473,7 @@ program DAMASK_grid solres%converged, solres%iterationsNeeded flush(statUnit) endif - elseif (cutBackLevel < maxCutBack) then ! further cutbacking tolerated? + elseif (cutBackLevel < num%maxCutBack) then ! further cutbacking tolerated? cutBack = .true. stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator cutBackLevel = cutBackLevel + 1 diff --git a/src/math.f90 b/src/math.f90 index 66076455f..3d686fe7a 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -72,6 +72,14 @@ module math 3,2, & 3,3 & ],shape(MAPPLAIN)) !< arrangement in Plain notation + + type, private :: tNumerics + integer :: & + randomSeed !< fixed seeding for pseudo-random number generator, Default 0: use random seed + + end type + + type(tNumerics), private :: num interface math_eye module procedure math_identity2nd @@ -91,8 +99,7 @@ subroutine math_init real(pReal), dimension(4) :: randTest integer :: & - randSize, & - randomSeed !< fixed seeding for pseudo-random number generator, Default 0: use random seed + randSize integer, dimension(:), allocatable :: randInit class(tNode), pointer :: & num_generic @@ -100,12 +107,12 @@ subroutine math_init write(6,'(/,a)') ' <<<+- math init -+>>>'; flush(6) num_generic => numerics_root%get('generic',defaultVal=emptyDict) - randomSeed = num_generic%get_asInt('random_seed', defaultVal = 0) + num%randomSeed = num_generic%get_asInt('random_seed', defaultVal = 0) call random_seed(size=randSize) allocate(randInit(randSize)) - if (randomSeed > 0) then - randInit = randomSeed + if (num%randomSeed > 0) then + randInit = num%randomSeed else call random_seed() call random_seed(get = randInit) From 6062cc43c4ef77554547cdaf4391a2f83468d676 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Wed, 24 Jun 2020 16:37:30 +0200 Subject: [PATCH 34/43] extending num structure to other modules; hard coding of tol variables was incorrect --- src/grid/grid_damage_spectral.f90 | 61 +++++++++++++++++-------------- 1 file changed, 33 insertions(+), 28 deletions(-) diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 659b111a3..5c0f98f9f 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -21,6 +21,18 @@ module grid_damage_spectral implicit none private + type, private :: tNumerics + integer :: & + itmax !< max number of iterations + real(pReal) :: & + residualStiffness, & !< non-zero residual damage + eps_damage_atol, & !< absolute tolerance for damage evolution + eps_damage_rtol !< relative tolerance for damage evolution + character(len=:), allocatable :: & + petsc_options + end type tNumerics + + type(tNumerics), private :: num !-------------------------------------------------------------------------------------------------- ! derived types type(tSolutionParams), private :: params @@ -61,10 +73,10 @@ subroutine grid_damage_spectral_init Vec :: uBound, lBound PetscErrorCode :: ierr class(tNode), pointer :: & - num_grid + num_grid, & + num_generic character(len=pStringLen) :: & - snes_type, & - petsc_options + snes_type write(6,'(/,a)') ' <<<+- grid_spectral_damage init -+>>>' @@ -72,16 +84,25 @@ subroutine grid_damage_spectral_init write(6,'(a)') ' https://doi.org/10.1007/978-981-10-6855-3_80' !------------------------------------------------------------------------------------------------- -! read numerical parameter +! read numerical parameters and do sanity checks num_grid => numerics_root%get('grid',defaultVal=emptyDict) - petsc_options = num_grid%get_asString('petsc_options',defaultVal='') + num%petsc_options = num_grid%get_asString('petsc_options',defaultVal='') + num%itmax = num_grid%get_asInt ('itmax',defaultVal=250) + num%eps_damage_atol = num_grid%get_asFloat ('eps_damage_atol',defaultVal=1.0e-2_pReal) + num%eps_damage_rtol = num_grid%get_asFloat ('eps_damage_rtol',defaultVal=1.0e-6_pReal) + + num_generic => numerics_root%get('generic',defaultVal=emptyDict) + num%residualStiffness = num_generic%get_asFloat('residualStiffness', defaultVal=1.0e-6_pReal) + + if (num%residualStiffness < 0.0_pReal) call IO_error(301,ext_msg='residualStiffness') + if (num%itmax <= 1) call IO_error(301,ext_msg='itmax') !-------------------------------------------------------------------------------------------------- ! set default and user defined options for PETSc call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-damage_snes_type newtonls -damage_snes_mf & &-damage_snes_ksp_ew -damage_ksp_type fgmres',ierr) CHKERRQ(ierr) - call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,num%petsc_options,ierr) CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- @@ -146,23 +167,14 @@ function grid_damage_spectral_solution(timeinc,timeinc_old) result(solution) real(pReal), intent(in) :: & timeinc, & !< increment in time for current solution timeinc_old !< increment in time of last increment - integer :: i, j, k, cell, & - itmax !< maximum number of iterations + integer :: i, j, k, cell type(tSolutionState) :: solution - class(tNode), pointer :: & - num_grid PetscInt :: devNull PetscReal :: phi_min, phi_max, stagNorm, solnNorm PetscErrorCode :: ierr SNESConvergedReason :: reason -!------------------------------------------------------------------- -! reading numerical parameter and do sanity check - num_grid => numerics_root%get('grid',defaultVal=emptyDict) - itmax = num_grid%get_asInt('itmax',defaultVal=250) - if (itmax <= 1) call IO_error(301,ext_msg='itmax') - solution%converged =.false. !-------------------------------------------------------------------------------------------------- @@ -175,7 +187,7 @@ function grid_damage_spectral_solution(timeinc,timeinc_old) result(solution) if (reason < 1) then solution%converged = .false. - solution%iterationsNeeded = itmax + solution%iterationsNeeded = num%itmax else solution%converged = .true. solution%iterationsNeeded = totalIter @@ -185,7 +197,7 @@ function grid_damage_spectral_solution(timeinc,timeinc_old) result(solution) call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,solnNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) phi_stagInc = phi_current - solution%stagConverged = stagNorm < max(1.0e-2_pReal, 1.0e-6_pReal*solnNorm) + solution%stagConverged = stagNorm < max(num%eps_damage_atol, num%eps_damage_rtol*solnNorm) !-------------------------------------------------------------------------------------------------- ! updating damage state @@ -256,14 +268,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) PetscObject :: dummy PetscErrorCode :: ierr integer :: i, j, k, cell - real(pReal) :: phiDot, dPhiDot_dPhi, mobility, & - residualStiffness !< non-zero residual damage - class(tNode), pointer :: & - num_generic - - num_generic => numerics_root%get('generic',defaultVal=emptyDict) - residualStiffness = num_generic%get_asFloat('residualStiffness', defaultVal=1.0e-6_pReal) - if (residualStiffness < 0.0_pReal) call IO_error(301,ext_msg='residualStiffness') + real(pReal) :: phiDot, dPhiDot_dPhi, mobility phi_current = x_scal !-------------------------------------------------------------------------------------------------- @@ -299,8 +304,8 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) call utilities_FFTscalarBackward where(scalarField_real(1:grid(1),1:grid(2),1:grid3) > phi_lastInc) & scalarField_real(1:grid(1),1:grid(2),1:grid3) = phi_lastInc - where(scalarField_real(1:grid(1),1:grid(2),1:grid3) < residualStiffness) & - scalarField_real(1:grid(1),1:grid(2),1:grid3) = residualStiffness + where(scalarField_real(1:grid(1),1:grid(2),1:grid3) < num%residualStiffness) & + scalarField_real(1:grid(1),1:grid(2),1:grid3) = num%residualStiffness !-------------------------------------------------------------------------------------------------- ! constructing residual From 434bfffc4649583f060989f2534d70d71a2636c0 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Wed, 24 Jun 2020 16:39:09 +0200 Subject: [PATCH 35/43] hard coding of tolerance variables in solvers not correct --- src/grid/grid_thermal_spectral.f90 | 22 ++++++++++++++++------ 1 file changed, 16 insertions(+), 6 deletions(-) diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index ee64566a6..acf3b1f65 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -26,6 +26,16 @@ module grid_thermal_spectral ! derived types type(tSolutionParams), private :: params + type, private :: tNumerics + real(pReal) :: & + eps_thermal_atol, & !< absolute tolerance for thermal equilibrium + eps_thermal_rtol !< relative tolerance for thermal equilibrium + character(len=:), allocatable :: & + petsc_options + end type tNumerics + + type(tNumerics), private :: num + !-------------------------------------------------------------------------------------------------- ! PETSc data SNES, private :: thermal_snes @@ -64,8 +74,6 @@ subroutine grid_thermal_spectral_init PetscErrorCode :: ierr class(tNode), pointer :: & num_grid - character(len=pStringLen) :: & - petsc_options write(6,'(/,a)') ' <<<+- grid_thermal_spectral init -+>>>' @@ -75,13 +83,15 @@ subroutine grid_thermal_spectral_init !------------------------------------------------------------------------------------------------- ! read numerical parameter num_grid => numerics_root%get('grid',defaultVal=emptyDict) - petsc_options = num_grid%get_asString('petsc_options',defaultVal='') - + num%petsc_options = num_grid%get_asString('petsc_options',defaultVal='') + num%eps_thermal_atol = num_grid%get_asFloat ('eps_thermal_atol',defaultVal=1.0e-2_pReal) + num%eps_thermal_rtol = num_grid%get_asFloat ('eps_thermal_rtol',defaultVal=1.0e-6_pReal) + !-------------------------------------------------------------------------------------------------- ! set default and user defined options for PETSc call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-thermal_snes_type ngmres',ierr) CHKERRQ(ierr) - call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,num%petsc_options,ierr) CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- @@ -182,7 +192,7 @@ function grid_thermal_spectral_solution(timeinc,timeinc_old) result(solution) call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,solnNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) T_stagInc = T_current - solution%stagConverged = stagNorm < max(1.0e-2_pReal, 1.0e-6_pReal*solnNorm) + solution%stagConverged = stagNorm < max(num%eps_thermal_atol, num%eps_thermal_rtol*solnNorm) !-------------------------------------------------------------------------------------------------- ! updating thermal state From be84561e2e41825e2b90b8e3de22149f1b2155fa Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Wed, 24 Jun 2020 16:45:13 +0200 Subject: [PATCH 36/43] num structure replicated; polishing --- src/crystallite.f90 | 4 +- src/grid/grid_mech_FEM.f90 | 91 ++++++------- src/grid/grid_mech_spectral_basic.f90 | 93 ++++++-------- src/grid/grid_mech_spectral_polarisation.f90 | 128 ++++++++----------- src/grid/spectral_utilities.f90 | 10 +- src/marc/discretization_marc.f90 | 4 +- src/mesh/DAMASK_mesh.f90 | 22 ++-- src/mesh/FEM_utilities.f90 | 6 +- src/mesh/mesh_mech_FEM.f90 | 71 +++++----- 9 files changed, 190 insertions(+), 239 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 9113fa12f..2454abaab 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -83,7 +83,7 @@ module crystallite iJacoLpresiduum, & !< frequency of Jacobian update of residuum in Lp nState, & !< state loop limit nStress !< stress loop limit - character(len=pStringLen) :: & + character(len=:), allocatable :: & integrator !< integrator scheme real(pReal) :: & subStepMinCryst, & !< minimum (relative) size of sub-step allowed during cutback @@ -196,7 +196,7 @@ subroutine crystallite_init if(num%nStress< 1) call IO_error(301,ext_msg='nStress') - select case(trim(num%integrator)) + select case(num%integrator) case('FPI') integrateState => integrateStateFPI case('Euler') diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index bdd76b423..775bc913d 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -30,6 +30,23 @@ module grid_mech_FEM ! derived types type(tSolutionParams), private :: params + type, private :: tNumerics + integer :: & + itmin, & !< minimum number of iterations + itmax !< maximum number of iterations + real(pReal) :: & + err_div, & + divTol, & + BCTol, & + eps_div_atol, & !< absolute tolerance for equilibrium + eps_div_rtol, & !< relative tolerance for equilibrium + eps_stress_atol, & !< absolute tolerance for fullfillment of stress BC + eps_stress_rtol !< relative tolerance for fullfillment of stress BC + character(len=:), allocatable :: & + petsc_options + end type tNumerics + + type(tNumerics), private :: num !-------------------------------------------------------------------------------------------------- ! PETSc data DM, private :: mech_grid @@ -97,8 +114,7 @@ subroutine grid_mech_FEM_init PetscErrorCode :: ierr integer(HID_T) :: fileHandle, groupHandle character(len=pStringLen) :: & - fileName, & - petsc_options + fileName class(tNode), pointer :: & num_grid real(pReal), dimension(3,3,3,3) :: devNull @@ -108,16 +124,30 @@ subroutine grid_mech_FEM_init write(6,'(/,a)') ' <<<+- grid_mech_FEM init -+>>>'; flush(6) !------------------------------------------------------------------------------------------------- -! read numerical parameter +! read numerical parameter and do sanity checks num_grid => numerics_root%get('grid',defaultVal=emptyDict) - petsc_options = num_grid%get_asString('petsc_options',defaultVal='') + num%petsc_options = num_grid%get_asString('petsc_options', defaultVal='') + num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal) + num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal) + num%eps_stress_atol = num_grid%get_asFloat ('eps_stress_atol', defaultVal=1.0e3_pReal) + num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol', defaultVal=0.01_pReal) + + num%itmin = num_grid%get_asInt ('itmin',defaultVal=1) + num%itmax = num_grid%get_asInt ('itmax',defaultVal=250) + + if (num%eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol') + if (num%eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol') + if (num%eps_stress_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_stress_atol') + if (num%eps_stress_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_stress_rtol') + if (num%itmax <= 1) call IO_error(301,ext_msg='itmax') + if (num%itmin > num%itmax .or. num%itmin < 1) call IO_error(301,ext_msg='itmin') !-------------------------------------------------------------------------------------------------- ! set default and user defined options for PETSc call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type newtonls -mech_ksp_type fgmres & &-mech_ksp_max_it 25 -mech_pc_type ml -mech_mg_levels_ksp_type chebyshev',ierr) CHKERRQ(ierr) - call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,num%petsc_options,ierr) CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- @@ -419,49 +449,23 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,i SNESConvergedReason :: reason PetscObject :: dummy PetscErrorCode :: ierr - integer :: & - itmin, & !< minimum number of iterations - itmax !< maximum number of iterations real(pReal) :: & err_div, & divTol, & - BCTol, & - eps_div_atol, & !< absolute tolerance for equilibrium - eps_div_rtol, & !< relative tolerance for equilibrium - eps_stress_atol, & !< absolute tolerance for fullfillment of stress BC - eps_stress_rtol !< relative tolerance for fullfillment of stress BC - class(tNode), pointer :: & - num_grid + BCTol -!----------------------------------------------------------------------------------- -! reading numerical parameters and do sanity check - num_grid => numerics_root%get('grid',defaultVal=emptyDict) - eps_div_atol = num_grid%get_asFloat('eps_div_atol',defaultVal=1.0e-4_pReal) - eps_div_rtol = num_grid%get_asFloat('eps_div_rtol',defaultVal=5.0e-4_pReal) - eps_stress_atol = num_grid%get_asFloat('eps_stress_atol',defaultVal=1.0e3_pReal) - eps_stress_rtol = num_grid%get_asFloat('eps_stress_rtol',defaultVal=0.01_pReal) - - itmin = num_grid%get_asInt('itmin',defaultVal=1) - itmax = num_grid%get_asInt('itmax',defaultVal=250) - - if (eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol') - if (eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol') - if (eps_stress_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_stress_atol') - if (eps_stress_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_stress_rtol') - if (itmax <= 1) call IO_error(301,ext_msg='itmax') - if (itmin > itmax .or. itmin < 1) call IO_error(301,ext_msg='itmin') !------------------------------------------------------------------------------------ err_div = fnorm*sqrt(wgt)*geomSize(1)/scaledGeomSize(1)/detJ - divTol = max(maxval(abs(P_av))*eps_div_rtol ,eps_div_atol) - BCTol = max(maxval(abs(P_av))*eps_stress_rtol,eps_stress_atol) + divTol = max(maxval(abs(P_av))*num%eps_div_rtol ,num%eps_div_atol) + BCTol = max(maxval(abs(P_av))*num%eps_stress_rtol,num%eps_stress_atol) - if ((totalIter >= itmin .and. & + if ((totalIter >= num%itmin .and. & all([ err_div/divTol, & err_BC /BCTol ] < 1.0_pReal)) & .or. terminallyIll) then reason = 1 - elseif (totalIter >= itmax) then + elseif (totalIter >= num%itmax) then reason = -1 else reason = 0 @@ -499,19 +503,6 @@ subroutine formResidual(da_local,x_local, & PetscObject :: dummy PetscErrorCode :: ierr real(pReal), dimension(3,3,3,3) :: devNull - integer :: & - itmin, & - itmax - class(tNode), pointer :: & - num_grid - -!---------------------------------------------------------------------- -! read numerical paramteters and do sanity checks - num_grid => numerics_root%get('grid',defaultVal=emptyDict) - itmin = num_grid%get_asInt('itmin',defaultVal=1) - itmax = num_grid%get_asInt('itmax',defaultVal=250) - if (itmax <= 1) call IO_error(301,ext_msg='itmax') - if (itmin > itmax .or. itmin < 1) call IO_error(301,ext_msg='itmin') call SNESGetNumberFunctionEvals(mech_snes,nfuncs,ierr); CHKERRQ(ierr) call SNESGetIterationNumber(mech_snes,PETScIter,ierr); CHKERRQ(ierr) @@ -522,7 +513,7 @@ subroutine formResidual(da_local,x_local, & ! begin of new iteration newIteration: if (totalIter <= PETScIter) then totalIter = totalIter + 1 - write(6,'(1x,a,3(a,i0))') trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter+1, '≤', itmax + write(6,'(1x,a,3(a,i0))') trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter+1, '≤', num%itmax if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 61ac18779..b1c4d085a 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -32,9 +32,19 @@ module grid_mech_spectral_basic type, private :: tNumerics logical :: update_gamma !< update gamma operator with current stiffness + integer :: & + itmin, & !< minimum number of iterations + itmax !< maximum number of iterations + real(pReal) :: & + eps_div_atol, & !< absolute tolerance for equilibrium + eps_div_rtol, & !< relative tolerance for equilibrium + eps_stress_atol, & !< absolute tolerance for fullfillment of stress BC + eps_stress_rtol !< relative tolerance for fullfillment of stress BC + character(len=:), allocatable :: & + petsc_options end type tNumerics - type(tNumerics) :: num ! numerics parameters. Better name? + type(tNumerics), private :: num ! numerics parameters. Better name? !-------------------------------------------------------------------------------------------------- ! PETSc data @@ -98,8 +108,7 @@ subroutine grid_mech_spectral_basic_init integer(HID_T) :: fileHandle, groupHandle integer :: fileUnit character(len=pStringLen) :: & - fileName, & - petsc_options + fileName write(6,'(/,a)') ' <<<+- grid_mech_spectral_basic init -+>>>'; flush(6) @@ -110,16 +119,31 @@ subroutine grid_mech_spectral_basic_init write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006' !------------------------------------------------------------------------------------------------- -! read numerical parameters +! read numerical parameters and do sanity checks num_grid => numerics_root%get('grid',defaultVal=emptyDict) - petsc_options = num_grid%get_asString('petsc_options',defaultVal='') - num%update_gamma = num_grid%get_asBool('update_gamma',defaultVal=.false.) - + + num%petsc_options = num_grid%get_asString ('petsc_options', defaultVal='') + num%update_gamma = num_grid%get_asBool ('update_gamma', defaultVal=.false.) + num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal) + num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal) + num%eps_stress_atol = num_grid%get_asFloat ('eps_stress_atol',defaultVal=1.0e3_pReal) + num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol',defaultVal=0.01_pReal) + + num%itmin = num_grid%get_asInt ('itmin',defaultVal=1) + num%itmax = num_grid%get_asInt ('itmax',defaultVal=250) + + if (num%eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol') + if (num%eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol') + if (num%eps_stress_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_stress_atol') + if (num%eps_stress_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_stress_rtol') + if (num%itmax <= 1) call IO_error(301,ext_msg='itmax') + if (num%itmin > num%itmax .or. num%itmin < 1) call IO_error(301,ext_msg='itmin') + !-------------------------------------------------------------------------------------------------- ! set default and user defined options for PETSc call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr) CHKERRQ(ierr) - call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,num%petsc_options,ierr) CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- @@ -389,47 +413,19 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm SNESConvergedReason :: reason PetscObject :: dummy PetscErrorCode :: ierr - integer :: & - itmin, & !< minimum number of iterations - itmax !< maximum number of iterations real(pReal) :: & divTol, & - BCTol, & - eps_div_atol, & !< absolute tolerance for equilibrium - eps_div_rtol, & !< relative tolerance for equilibrium - eps_stress_atol, & !< absolute tolerance for fullfillment of stress BC - eps_stress_rtol !< relative tolerance for fullfillment of stress BC - class(tNode), pointer :: & - num_grid + BCTol -!----------------------------------------------------------------------------------- -! reading numerical parameters and do sanity check - num_grid => numerics_root%get('grid',defaultVal=emptyDict) - eps_div_atol = num_grid%get_asFloat('eps_div_atol',defaultVal=1.0e-4_pReal) - eps_div_rtol = num_grid%get_asFloat('eps_div_rtol',defaultVal=5.0e-4_pReal) - eps_stress_atol = num_grid%get_asFloat('eps_stress_atol',defaultVal=1.0e3_pReal) - eps_stress_rtol = num_grid%get_asFloat('eps_stress_rtol',defaultVal=0.01_pReal) + divTol = max(maxval(abs(P_av))*num%eps_div_rtol ,num%eps_div_atol) + BCTol = max(maxval(abs(P_av))*num%eps_stress_rtol,num%eps_stress_atol) - itmin = num_grid%get_asInt('itmin',defaultVal=1) - itmax = num_grid%get_asInt('itmax',defaultVal=250) - - if (eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol') - if (eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol') - if (eps_stress_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_stress_atol') - if (eps_stress_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_stress_rtol') - if (itmax <= 1) call IO_error(301,ext_msg='itmax') - if (itmin > itmax .or. itmin < 1) call IO_error(301,ext_msg='itmin') -!------------------------------------------------------------------------------------ - - divTol = max(maxval(abs(P_av))*eps_div_rtol ,eps_div_atol) - BCTol = max(maxval(abs(P_av))*eps_stress_rtol,eps_stress_atol) - - if ((totalIter >= itmin .and. & + if ((totalIter >= num%itmin .and. & all([ err_div/divTol, & err_BC /BCTol ] < 1.0_pReal)) & .or. terminallyIll) then reason = 1 - elseif (totalIter >= itmax) then + elseif (totalIter >= num%itmax) then reason = -1 else reason = 0 @@ -466,19 +462,6 @@ subroutine formResidual(in, F, & nfuncs PetscObject :: dummy PetscErrorCode :: ierr - integer :: & - itmin, & - itmax - class(tNode), pointer :: & - num_grid - -!---------------------------------------------------------------------- -! read numerical paramteter and do sanity checks - num_grid => numerics_root%get('grid',defaultVal=emptyDict) - itmin = num_grid%get_asInt('itmin',defaultVal=1) - itmax = num_grid%get_asInt('itmax',defaultVal=250) - if (itmax <= 1) call IO_error(301,ext_msg='itmax') - if (itmin > itmax .or. itmin < 1) call IO_error(301,ext_msg='itmin') call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) @@ -488,7 +471,7 @@ subroutine formResidual(in, F, & ! begin of new iteration newIteration: if (totalIter <= PETScIter) then totalIter = totalIter + 1 - write(6,'(1x,a,3(a,i0))') trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax + write(6,'(1x,a,3(a,i0))') trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index a32216615..87cb69579 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -33,9 +33,24 @@ module grid_mech_spectral_polarisation type, private :: tNumerics logical :: update_gamma !< update gamma operator with current stiffness + character(len=:), allocatable :: & + petsc_options + integer :: & + itmin, & !< minimum number of iterations + itmax !< maximum number of iterations + real(pReal) :: & + eps_div_atol, & !< absolute tolerance for equilibrium + eps_div_rtol, & !< relative tolerance for equilibrium + eps_curl_atol, & !< absolute tolerance for compatibility + eps_curl_rtol, & !< relative tolerance for compatibility + eps_stress_atol, & !< absolute tolerance for fullfillment of stress BC + eps_stress_rtol !< relative tolerance for fullfillment of stress BC + real(pReal) :: & + polarAlpha, & !< polarization scheme parameter 0.0 < alpha < 2.0. alpha = 1.0 ==> AL scheme, alpha = 2.0 ==> accelerated scheme + polarBeta !< polarization scheme parameter 0.0 < beta < 2.0. beta = 1.0 ==> AL scheme, beta = 2.0 ==> accelerated scheme end type tNumerics - type(tNumerics) :: num ! numerics parameters. Better name? + type(tNumerics), private :: num ! numerics parameters. Better name? !-------------------------------------------------------------------------------------------------- ! PETSc data @@ -107,8 +122,7 @@ subroutine grid_mech_spectral_polarisation_init integer(HID_T) :: fileHandle, groupHandle integer :: fileUnit character(len=pStringLen) :: & - fileName, & - petsc_options + fileName write(6,'(/,a)') ' <<<+- grid_mech_spectral_polarisation init -+>>>'; flush(6) @@ -118,14 +132,36 @@ subroutine grid_mech_spectral_polarisation_init !------------------------------------------------------------------------------------------------- ! read numerical parameters num_grid => numerics_root%get('grid',defaultVal=emptyDict) - petsc_options = num_grid%get_asString('petsc_options',defaultVal='') - num%update_gamma = num_grid%get_asBool('update_gamma',defaultVal=.false.) + + num%petsc_options = num_grid%get_asString('petsc_options', defaultVal='') + num%update_gamma = num_grid%get_asBool ('update_gamma', defaultVal=.false.) + num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal) + num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal) + num%eps_curl_atol = num_grid%get_asFloat ('eps_curl_atol', defaultVal=1.0e-10_pReal) + num%eps_curl_rtol = num_grid%get_asFloat ('eps_curl_rtol', defaultVal=5.0e-4_pReal) + num%eps_stress_atol = num_grid%get_asFloat ('eps_stress_atol', defaultVal=1.0e3_pReal) + num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol', defaultVal=0.01_pReal) + num%itmin = num_grid%get_asInt ('itmin', defaultVal=1) + num%itmax = num_grid%get_asInt ('itmax', defaultVal=250) + num%polarAlpha = num_grid%get_asFloat ('polaralpha', defaultVal=1.0_pReal) + num%polarBeta = num_grid%get_asFloat ('polarbeta', defaultVal=1.0_pReal) + + if (num%eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol') + if (num%eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol') + if (num%eps_curl_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_curl_atol') + if (num%eps_curl_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_curl_rtol') + if (num%eps_stress_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_stress_atol') + if (num%eps_stress_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_stress_rtol') + if (num%itmax <= 1) call IO_error(301,ext_msg='itmax') + if (num%itmin > num%itmax .or. num%itmin < 1) call IO_error(301,ext_msg='itmin') + if (num%polarAlpha <= 0.0_pReal .or. num%polarAlpha > 2.0_pReal) call IO_error(301,ext_msg='polarAlpha') + if (num%polarBeta < 0.0_pReal .or. num%polarBeta > 2.0_pReal) call IO_error(301,ext_msg='polarBeta') !-------------------------------------------------------------------------------------------------- ! set default and user defined options for PETSc call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr) CHKERRQ(ierr) - call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,num%petsc_options,ierr) CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- @@ -436,56 +472,22 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm SNESConvergedReason :: reason PetscObject :: dummy PetscErrorCode :: ierr - integer :: & - itmin, & !< minimum number of iterations - itmax !< maximum number of iterations real(pReal) :: & curlTol, & divTol, & - BCTol, & - eps_div_atol, & !< absolute tolerance for equilibrium - eps_div_rtol, & !< relative tolerance for equilibrium - eps_curl_atol, & !< absolute tolerance for compatibility - eps_curl_rtol, & !< relative tolerance for compatibility - eps_stress_atol, & !< absolute tolerance for fullfillment of stress BC - eps_stress_rtol !< relative tolerance for fullfillment of stress BC - class(tNode), pointer :: & - num_grid + BCTol -!----------------------------------------------------------------------------------- -! reading numerical parameters and do sanity check - num_grid => numerics_root%get('grid',defaultVal=emptyDict) - eps_div_atol = num_grid%get_asFloat('eps_div_atol',defaultVal=1.0e-4_pReal) - eps_div_rtol = num_grid%get_asFloat('eps_div_rtol',defaultVal=5.0e-4_pReal) - eps_curl_atol = num_grid%get_asFloat('eps_curl_atol',defaultVal=1.0e-10_pReal) - eps_curl_rtol = num_grid%get_asFloat('eps_curl_rtol',defaultVal=5.0e-4_pReal) - eps_stress_atol = num_grid%get_asFloat('eps_stress_atol',defaultVal=1.0e3_pReal) - eps_stress_rtol = num_grid%get_asFloat('eps_stress_rtol',defaultVal=0.01_pReal) + curlTol = max(maxval(abs(F_aim-math_I3))*num%eps_curl_rtol ,num%eps_curl_atol) + divTol = max(maxval(abs(P_av)) *num%eps_div_rtol ,num%eps_div_atol) + BCTol = max(maxval(abs(P_av)) *num%eps_stress_rtol,num%eps_stress_atol) - itmin = num_grid%get_asInt('itmin',defaultVal=1) - itmax = num_grid%get_asInt('itmax',defaultVal=250) - - if (eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol') - if (eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol') - if (eps_curl_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_curl_atol') - if (eps_curl_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_curl_rtol') - if (eps_stress_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_stress_atol') - if (eps_stress_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_stress_rtol') - if (itmax <= 1) call IO_error(301,ext_msg='itmax') - if (itmin > itmax .or. itmin < 1) call IO_error(301,ext_msg='itmin') -!------------------------------------------------------------------------------------ - - curlTol = max(maxval(abs(F_aim-math_I3))*eps_curl_rtol ,eps_curl_atol) - divTol = max(maxval(abs(P_av)) *eps_div_rtol ,eps_div_atol) - BCTol = max(maxval(abs(P_av)) *eps_stress_rtol,eps_stress_atol) - - if ((totalIter >= itmin .and. & + if ((totalIter >= num%itmin .and. & all([ err_div /divTol, & err_curl/curlTol, & err_BC /BCTol ] < 1.0_pReal)) & .or. terminallyIll) then reason = 1 - elseif (totalIter >= itmax) then + elseif (totalIter >= num%itmax) then reason = -1 else reason = 0 @@ -527,27 +529,9 @@ subroutine formResidual(in, FandF_tau, & nfuncs PetscObject :: dummy PetscErrorCode :: ierr - class(tNode), pointer :: & - num_grid - real(pReal) :: & - polarAlpha, & !< polarization scheme parameter 0.0 < alpha < 2.0. alpha = 1.0 ==> AL scheme, alpha = 2.0 ==> accelerated scheme - polarBeta !< polarization scheme parameter 0.0 < beta < 2.0. beta = 1.0 ==> AL scheme, beta = 2.0 ==> accelerated scheme - integer :: & - i, j, k, e, & - itmin, itmax + integer :: & + i, j, k, e -!-------------------------------------------------------------------------------------------------- -! read numerical paramteters and do sanity checks - num_grid => numerics_root%get('grid',defaultVal = emptyDict) - polarAlpha = num_grid%get_asFloat('polaralpha',defaultVal=1.0_pReal) - polarBeta = num_grid%get_asFloat('polarbeta', defaultVal=1.0_pReal) - itmin = num_grid%get_asInt('itmin',defaultVal=1) - itmax = num_grid%get_asInt('itmax',defaultVal=250) - - if (itmax <= 1) call IO_error(301,ext_msg='itmax') - if (itmin > itmax .or. itmin < 1) call IO_error(301,ext_msg='itmin') - if (polarAlpha <= 0.0_pReal .or. polarAlpha > 2.0_pReal) call IO_error(301,ext_msg='polarAlpha') - if (polarBeta < 0.0_pReal .or. polarBeta > 2.0_pReal) call IO_error(301,ext_msg='polarBeta') !--------------------------------------------------------------------------------------------------- F => FandF_tau(1:3,1:3,1,& @@ -570,7 +554,7 @@ subroutine formResidual(in, FandF_tau, & ! begin of new iteration newIteration: if (totalIter <= PETScIter) then totalIter = totalIter + 1 - write(6,'(1x,a,3(a,i0))') trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax + write(6,'(1x,a,3(a,i0))') trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) @@ -584,26 +568,26 @@ subroutine formResidual(in, FandF_tau, & tensorField_real = 0.0_pReal do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) tensorField_real(1:3,1:3,i,j,k) = & - polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& - polarAlpha*matmul(F(1:3,1:3,i,j,k), & + num%polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& + num%polarAlpha*matmul(F(1:3,1:3,i,j,k), & math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3)) enddo; enddo; enddo !-------------------------------------------------------------------------------------------------- ! doing convolution in Fourier space call utilities_FFTtensorForward - call utilities_fourierGammaConvolution(params%rotation_BC%rotate(polarBeta*F_aim,active=.true.)) + call utilities_fourierGammaConvolution(params%rotation_BC%rotate(num%polarBeta*F_aim,active=.true.)) call utilities_FFTtensorBackward !-------------------------------------------------------------------------------------------------- ! constructing residual - residual_F_tau = polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) + residual_F_tau = num%polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response call utilities_constitutiveResponse(residual_F, & ! "residuum" gets field of first PK stress (to save memory) P_av,C_volAvg,C_minMaxAvg, & - F - residual_F_tau/polarBeta,params%timeinc,params%rotation_BC) + F - residual_F_tau/num%polarBeta,params%timeinc,params%rotation_BC) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) !-------------------------------------------------------------------------------------------------- diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index cc9d3dd65..08e5d633c 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -121,10 +121,10 @@ module spectral_utilities character(len=:), allocatable :: & spectral_derivative, & !< approximation used for derivatives in Fourier space FFTW_plan_mode, & !< FFTW plan mode, see www.fftw.org - PETSc_options + petsc_options end type tNumerics - type(tNumerics) :: num ! numerics parameters. Better name? + type(tNumerics), private :: num ! numerics parameters. Better name? enum, bind(c); enumerator :: & DERIVATIVE_CONTINUOUS_ID, & @@ -189,8 +189,6 @@ subroutine spectral_utilities_init scalarSize = 1_C_INTPTR_T, & vecSize = 3_C_INTPTR_T, & tensorSize = 9_C_INTPTR_T - character(len=pStringLen) :: & - petsc_options class(tNode), pointer :: & num_grid @@ -220,13 +218,13 @@ subroutine spectral_utilities_init ' add more using the PETSc_Options keyword in numerics.config '; flush(6) num_grid => numerics_root%get('grid',defaultVal=emptyDict) - petsc_options = num_grid%get_asString('petsc_options',defaultVal='') + num%petsc_options = num_grid%get_asString('petsc_options',defaultVal='') call PETScOptionsClear(PETSC_NULL_OPTIONS,ierr) CHKERRQ(ierr) if(debugPETSc) call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),ierr) CHKERRQ(ierr) - call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,num%petsc_options,ierr) CHKERRQ(ierr) grid1Red = grid(1)/2 + 1 diff --git a/src/marc/discretization_marc.f90 b/src/marc/discretization_marc.f90 index 8945c2c18..a4a187c38 100644 --- a/src/marc/discretization_marc.f90 +++ b/src/marc/discretization_marc.f90 @@ -34,8 +34,8 @@ module discretization_marc mesh_unitlength !< physical length of one unit in mesh integer, dimension(:), allocatable, public :: & - mesh_FEM2DAMASK_elem, & !< DAMASK element ID for Marc element ID - mesh_FEM2DAMASK_node !< DAMASK node ID for Marc node ID + mesh_FEM2DAMASK_elem, & !< DAMASK element ID for Marc element ID + mesh_FEM2DAMASK_node !< DAMASK node ID for Marc node ID public :: & discretization_marc_init diff --git a/src/mesh/DAMASK_mesh.f90 b/src/mesh/DAMASK_mesh.f90 index 72fa53566..185f4c375 100644 --- a/src/mesh/DAMASK_mesh.f90 +++ b/src/mesh/DAMASK_mesh.f90 @@ -49,7 +49,6 @@ program DAMASK_mesh i, & errorID, & cutBackLevel = 0, & !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$ - maxCutBack, & !< max number of cutbacks stepFraction = 0, & !< fraction of current time interval currentLoadcase = 0, & !< current load case currentFace = 0, & @@ -57,7 +56,6 @@ program DAMASK_mesh totalIncsCounter = 0, & !< total # of increments statUnit = 0, & !< file unit for statistics output stagIter, & - stagItMax, & !< max number of field level staggered iterations component class(tNode), pointer :: & num_mesh @@ -65,6 +63,14 @@ program DAMASK_mesh character(len=pStringLen) :: & incInfo, & loadcase_string + type :: tNumerics + integer :: & + stagItMax, & !< max number of field level staggered iterations + maxCutBack !< max number of cutbacks + end type tNumerics + + type(tNumerics) :: num + type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases type(tSolutionState), allocatable, dimension(:) :: solres PetscInt :: faceSet, currentFaceSet, field, dimPlex @@ -81,11 +87,11 @@ program DAMASK_mesh !--------------------------------------------------------------------- ! reading field information from numerics file and do sanity checks num_mesh => numerics_root%get('mesh', defaultVal=emptyDict) - stagItMax = num_mesh%get_asInt('maxStaggeredIter',defaultVal=10) - maxCutBack = num_mesh%get_asInt('maxCutBack',defaultVal=3) + num%stagItMax = num_mesh%get_asInt('maxStaggeredIter',defaultVal=10) + num%maxCutBack = num_mesh%get_asInt('maxCutBack',defaultVal=3) - if (stagItMax < 0) call IO_error(301,ext_msg='maxStaggeredIter') - if (maxCutBack < 0) call IO_error(301,ext_msg='maxCutBack') + if (num%stagItMax < 0) call IO_error(301,ext_msg='maxStaggeredIter') + if (num%maxCutBack < 0) call IO_error(301,ext_msg='maxCutBack') ! reading basic information from load case file and allocate data structure containing load cases call DMGetDimension(geomMesh,dimPlex,ierr); CHKERRA(ierr) !< dimension of mesh (2D or 3D) @@ -327,7 +333,7 @@ program DAMASK_mesh enddo stagIter = stagIter + 1 - stagIterate = stagIter < stagItMax & + stagIterate = stagIter < num%stagItMax & .and. all(solres(:)%converged) & .and. .not. all(solres(:)%stagConverged) ! stationary with respect to staggered iteration enddo @@ -335,7 +341,7 @@ program DAMASK_mesh ! check solution cutBack = .False. if(.not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found - if (cutBackLevel < maxCutBack) then ! do cut back + if (cutBackLevel < num%maxCutBack) then ! do cut back write(6,'(/,a)') ' cut back detected' cutBack = .True. stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index a8479913e..a6236175c 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -105,14 +105,14 @@ subroutine FEM_utilities_init class(tNode), pointer :: & num_mesh integer :: structOrder !< order of displacement shape functions - character(len=pStringLen) :: & + character(len=:), allocatable :: & petsc_options PetscErrorCode :: ierr write(6,'(/,a)') ' <<<+- DAMASK_FEM_utilities init -+>>>' num_mesh => numerics_root%get('mesh',defaultVal=emptyDict) - structOrder = num_mesh%get_asInt('structOrder',defaultVal = 2) + structOrder = num_mesh%get_asInt ('structOrder', defaultVal = 2) petsc_options = num_mesh%get_asString('petsc_options', defaultVal='') !-------------------------------------------------------------------------------------------------- @@ -134,7 +134,7 @@ subroutine FEM_utilities_init &-mech_pc_type ml -mech_mg_levels_ksp_type chebyshev & &-mech_mg_levels_pc_type sor -mech_pc_ml_nullspace user',ierr) CHKERRQ(ierr) - call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) + call PetscOptionsInsertString(PETSC_NULL_OPTIONS,petsc_options,ierr) CHKERRQ(ierr) write(petsc_optionsOrder,'(a,i0)') '-mechFE_petscspace_degree ', structOrder call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsOrder),ierr) diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 715b02c0d..918914432 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -37,6 +37,18 @@ module mesh_mech_FEM type(tSolutionParams) :: params + type, private :: tNumerics + integer :: & + integrationOrder, & !< order of quadrature rule required + itmax + logical :: & + BBarStabilisation + real(pReal) :: & + eps_struct_atol, & !< absolute tolerance for mechanical equilibrium + eps_struct_rtol !< relative tolerance for mechanical equilibrium + end type tNumerics + + type(tNumerics), private :: num !-------------------------------------------------------------------------------------------------- ! PETSc data SNES :: mech_snes @@ -97,20 +109,21 @@ subroutine FEM_mech_init(fieldBC) class(tNode), pointer :: & num_mesh - integer :: & - integrationOrder, & !< order of quadrature rule required - itmax write(6,'(/,a)') ' <<<+- FEM_mech init -+>>>'; flush(6) !----------------------------------------------------------------------------- ! read numerical parametes and do sanity checks num_mesh => numerics_root%get('mesh',defaultVal=emptyDict) - integrationOrder = num_mesh%get_asInt('integrationorder',defaultVal = 2) + num%integrationOrder = num_mesh%get_asInt('integrationorder',defaultVal = 2) + num%itmax = num_mesh%get_asInt('itmax',defaultVal=250) + num%BBarStabilisation = num_mesh%get_asBool('bbarstabilisation',defaultVal = .false.) + num%eps_struct_atol = num_mesh%get_asFloat('eps_struct_atol', defaultVal = 1.0e-10_pReal) + num%eps_struct_rtol = num_mesh%get_asFloat('eps_struct_rtol', defaultVal = 1.0e-4_pReal) - num_mesh => numerics_root%get('mesh',defaultVal=emptyDict) - itmax = num_mesh%get_asInt('itmax',defaultVal=250) - if (itmax <= 1) call IO_error(301,ext_msg='itmax') + if (num%itmax <= 1) call IO_error(301,ext_msg='itmax') + if (num%eps_struct_rtol <= 0.0_pReal) call IO_error(301,ext_msg='eps_struct_rtol') + if (num%eps_struct_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_struct_atol') !-------------------------------------------------------------------------------------------------- ! Setup FEM mech mesh @@ -119,9 +132,9 @@ subroutine FEM_mech_init(fieldBC) !-------------------------------------------------------------------------------------------------- ! Setup FEM mech discretization - qPoints = FEM_quadrature_points( dimPlex,integrationOrder)%p - qWeights = FEM_quadrature_weights(dimPlex,integrationOrder)%p - nQuadrature = FEM_nQuadrature( dimPlex,integrationOrder) + qPoints = FEM_quadrature_points( dimPlex,num%integrationOrder)%p + qWeights = FEM_quadrature_weights(dimPlex,num%integrationOrder)%p + nQuadrature = FEM_nQuadrature( dimPlex,num%integrationOrder) qPointsP => qPoints qWeightsP => qWeights call PetscQuadratureCreate(PETSC_COMM_SELF,mechQuad,ierr); CHKERRQ(ierr) @@ -130,7 +143,7 @@ subroutine FEM_mech_init(fieldBC) call PetscQuadratureSetData(mechQuad,dimPlex,nc,nQuadrature,qPointsP,qWeightsP,ierr) CHKERRQ(ierr) call PetscFECreateDefault(PETSC_COMM_SELF,dimPlex,nc,PETSC_TRUE,prefix, & - integrationOrder,mechFE,ierr); CHKERRQ(ierr) + num%integrationOrder,mechFE,ierr); CHKERRQ(ierr) call PetscFESetQuadrature(mechFE,mechQuad,ierr); CHKERRQ(ierr) call PetscFEGetDimension(mechFE,nBasis,ierr); CHKERRQ(ierr) nBasis = nBasis/nc @@ -221,7 +234,7 @@ subroutine FEM_mech_init(fieldBC) call SNESSetMaxLinearSolveFailures(mech_snes, huge(1), ierr); CHKERRQ(ierr) !< ignore linear solve failures call SNESSetConvergenceTest(mech_snes,FEM_mech_converged,PETSC_NULL_VEC,PETSC_NULL_FUNCTION,ierr) CHKERRQ(ierr) - call SNESSetTolerances(mech_snes,1.0,0.0,0.0,itmax,itmax,ierr) + call SNESSetTolerances(mech_snes,1.0,0.0,0.0,num%itmax,num%itmax,ierr) CHKERRQ(ierr) call SNESSetFromOptions(mech_snes,ierr); CHKERRQ(ierr) @@ -281,22 +294,10 @@ type(tSolutionState) function FEM_mech_solution( & character(len=*), intent(in) :: & incInfoIn -!-------------------------------------------------------------------------------------------------- - integer :: & - itmax - class(tNode), pointer :: & - num_mesh PetscErrorCode :: ierr SNESConvergedReason :: reason -!-------------------------------------------------------------------------------------------------- -! read numerical parameter and do sanity check - num_mesh => numerics_root%get('mesh',defaultVal=emptyDict) - itmax = num_mesh%get_asInt('itmax',defaultVal=250) - if (itmax <= 1) call IO_error(301,ext_msg='itmax') -!------------------------------------------------------------------------------------------------- - incInfo = incInfoIn FEM_mech_solution%converged =.false. !-------------------------------------------------------------------------------------------------- @@ -311,7 +312,7 @@ type(tSolutionState) function FEM_mech_solution( & if (reason < 1) then ! 0: still iterating (will not occur), negative -> convergence error FEM_mech_solution%converged = .false. - FEM_mech_solution%iterationsNeeded = itmax + FEM_mech_solution%iterationsNeeded = num%itmax else ! >= 1 proper convergence (or terminally ill) FEM_mech_solution%converged = .true. call SNESGetIterationNumber(mech_snes,FEM_mech_solution%iterationsNeeded,ierr) @@ -349,12 +350,6 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr) PetscInt :: bcSize IS :: bcPoints - class(tNode), pointer :: & - num_mesh - logical :: BBarStabilisation - - num_mesh => numerics_root%get('mesh',defaultVal=emptyDict) - BBarStabilisation = num_mesh%get_asBool('bbarstabilisation',defaultVal = .false.) allocate(pV0(dimPlex)) allocate(pcellJ(dimPlex**2)) @@ -409,7 +404,7 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr) materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1) = & reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1]) enddo - if (BBarStabilisation) then + if (num%BBarStabilisation) then detFAvg = math_det33(sum(materialpoint_F(1:3,1:3,1:nQuadrature,cell+1),dim=3)/real(nQuadrature)) do qPt = 1, nQuadrature materialpoint_F(1:dimPlex,1:dimPlex,qPt,cell+1) = & @@ -498,12 +493,6 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr) IS :: bcPoints - class(tNode), pointer :: & - num_mesh - logical :: BBarStabilisation - - num_mesh => numerics_root%get('mesh',defaultVal=emptyDict) - BBarStabilisation = num_mesh%get_asBool('bbarstabilisation',defaultVal = .false.) allocate(pV0(dimPlex)) allocate(pcellJ(dimPlex**2)) @@ -560,7 +549,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr) MatA = matmul(reshape(reshape(materialpoint_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,qPt+1,cell+1), & shape=[dimPlex,dimPlex,dimPlex,dimPlex], order=[2,1,4,3]), & shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1) - if (BBarStabilisation) then + if (num%BBarStabilisation) then F(1:dimPlex,1:dimPlex) = reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex]) FInv = math_inv33(F) K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0/real(dimPlex)) @@ -577,7 +566,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr) K_eA = K_eA + matmul(transpose(BMat),MatA) endif enddo - if (BBarStabilisation) then + if (num%BBarStabilisation) then FInv = math_inv33(FAvg) K_e = K_eA*math_det33(FAvg/real(nQuadrature))**(1.0/real(dimPlex)) + & (matmul(matmul(transpose(BMatAvg), & @@ -687,7 +676,7 @@ subroutine FEM_mech_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dumm !-------------------------------------------------------------------------------------------------- ! report - divTol = max(maxval(abs(P_av(1:dimPlex,1:dimPlex)))*1.0e-4_pReal,1.0e-10_pReal) + divTol = max(maxval(abs(P_av(1:dimPlex,1:dimPlex)))*num%eps_struct_rtol,num%eps_struct_atol) call SNESConvergedDefault(snes_local,PETScIter,xnorm,snorm,fnorm/divTol,reason,dummy,ierr) CHKERRQ(ierr) if (terminallyIll) reason = SNES_DIVERGED_FUNCTION_DOMAIN From 445d8b4f745527b9433b254deb9b154afc51a40f Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Wed, 24 Jun 2020 17:09:15 +0200 Subject: [PATCH 37/43] sanity checks --- src/grid/grid_damage_spectral.f90 | 2 ++ src/grid/grid_thermal_spectral.f90 | 7 +++++-- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 5c0f98f9f..eec65a89f 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -96,6 +96,8 @@ subroutine grid_damage_spectral_init if (num%residualStiffness < 0.0_pReal) call IO_error(301,ext_msg='residualStiffness') if (num%itmax <= 1) call IO_error(301,ext_msg='itmax') + if (num%eps_damage_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_damage_atol') + if (num%eps_damage_rtol <= 0.0_pReal) call IO_error(301,ext_msg='eps_damage_rtol') !-------------------------------------------------------------------------------------------------- ! set default and user defined options for PETSc diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index acf3b1f65..b15167cee 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -81,12 +81,15 @@ subroutine grid_thermal_spectral_init write(6,'(a)') ' https://doi.org/10.1007/978-981-10-6855-3_80' !------------------------------------------------------------------------------------------------- -! read numerical parameter +! read numerical parameter and do sanity checks num_grid => numerics_root%get('grid',defaultVal=emptyDict) num%petsc_options = num_grid%get_asString('petsc_options',defaultVal='') num%eps_thermal_atol = num_grid%get_asFloat ('eps_thermal_atol',defaultVal=1.0e-2_pReal) num%eps_thermal_rtol = num_grid%get_asFloat ('eps_thermal_rtol',defaultVal=1.0e-6_pReal) - + + if (num%eps_thermal_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_thermal_atol') + if (num%eps_thermal_rtol <= 0.0_pReal) call IO_error(301,ext_msg='eps_thermal_rtol') + !-------------------------------------------------------------------------------------------------- ! set default and user defined options for PETSc call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-thermal_snes_type ngmres',ierr) From a71b8aea5caffa0fed3b01dd950c25dc72cc058d Mon Sep 17 00:00:00 2001 From: Test User Date: Wed, 24 Jun 2020 17:29:13 +0200 Subject: [PATCH 38/43] [skip ci] updated version information after successful test of v2.0.3-2692-g3d93a5ff --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index d041f794f..d94a70ab1 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-2680-g17212b27 +v2.0.3-2692-g3d93a5ff From e155bef9a547738794480b4acb0e9c640f912083 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Wed, 24 Jun 2020 17:18:16 +0200 Subject: [PATCH 39/43] better names; polishing --- src/YAML_parse.f90 | 10 ++++----- src/grid/grid_mech_spectral_polarisation.f90 | 22 ++++++++++---------- src/numerics.f90 | 2 +- 3 files changed, 17 insertions(+), 17 deletions(-) diff --git a/src/YAML_parse.f90 b/src/YAML_parse.f90 index cfbcf11f3..626a387fa 100644 --- a/src/YAML_parse.f90 +++ b/src/YAML_parse.f90 @@ -36,15 +36,15 @@ end subroutine YAML_init !-------------------------------------------------------------------------------------------------- recursive function parse_flow(flow_string) result(node) - character(len=*), intent(inout) :: flow_string + character(len=*), intent(inout) :: flow_string !< YAML file in flow style class (tNode), pointer :: node class (tNode), pointer :: myVal character(len=pStringLen) :: key - integer :: e, & !> end position of dictionary or list - s, & !> start position of dictionary or list - d !> position of key: value separator (':') + integer :: e, & ! end position of dictionary or list + s, & ! start position of dictionary or list + d ! position of key: value separator (':') flow_string = trim(adjustl(flow_string(:))) if (len_trim(flow_string) == 0) then @@ -96,7 +96,7 @@ end function parse_flow !-------------------------------------------------------------------------------------------------- integer function find_end(str,e_char) - character(len=*), intent(in) :: str + character(len=*), intent(in) :: str !< chunk of YAML flow string character, intent(in) :: e_char !< end of list/dict ( '}' or ']') integer :: N_sq, & !< number of open square brackets diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 87cb69579..f6e554ed4 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -46,8 +46,8 @@ module grid_mech_spectral_polarisation eps_stress_atol, & !< absolute tolerance for fullfillment of stress BC eps_stress_rtol !< relative tolerance for fullfillment of stress BC real(pReal) :: & - polarAlpha, & !< polarization scheme parameter 0.0 < alpha < 2.0. alpha = 1.0 ==> AL scheme, alpha = 2.0 ==> accelerated scheme - polarBeta !< polarization scheme parameter 0.0 < beta < 2.0. beta = 1.0 ==> AL scheme, beta = 2.0 ==> accelerated scheme + alpha, & !< polarization scheme parameter 0.0 < alpha < 2.0. alpha = 1.0 ==> AL scheme, alpha = 2.0 ==> accelerated scheme + beta !< polarization scheme parameter 0.0 < beta < 2.0. beta = 1.0 ==> AL scheme, beta = 2.0 ==> accelerated scheme end type tNumerics type(tNumerics), private :: num ! numerics parameters. Better name? @@ -143,8 +143,8 @@ subroutine grid_mech_spectral_polarisation_init num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol', defaultVal=0.01_pReal) num%itmin = num_grid%get_asInt ('itmin', defaultVal=1) num%itmax = num_grid%get_asInt ('itmax', defaultVal=250) - num%polarAlpha = num_grid%get_asFloat ('polaralpha', defaultVal=1.0_pReal) - num%polarBeta = num_grid%get_asFloat ('polarbeta', defaultVal=1.0_pReal) + num%alpha = num_grid%get_asFloat ('alpha', defaultVal=1.0_pReal) + num%beta = num_grid%get_asFloat ('beta', defaultVal=1.0_pReal) if (num%eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol') if (num%eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol') @@ -154,8 +154,8 @@ subroutine grid_mech_spectral_polarisation_init if (num%eps_stress_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_stress_rtol') if (num%itmax <= 1) call IO_error(301,ext_msg='itmax') if (num%itmin > num%itmax .or. num%itmin < 1) call IO_error(301,ext_msg='itmin') - if (num%polarAlpha <= 0.0_pReal .or. num%polarAlpha > 2.0_pReal) call IO_error(301,ext_msg='polarAlpha') - if (num%polarBeta < 0.0_pReal .or. num%polarBeta > 2.0_pReal) call IO_error(301,ext_msg='polarBeta') + if (num%alpha <= 0.0_pReal .or. num%alpha > 2.0_pReal) call IO_error(301,ext_msg='alpha') + if (num%beta < 0.0_pReal .or. num%beta > 2.0_pReal) call IO_error(301,ext_msg='beta') !-------------------------------------------------------------------------------------------------- ! set default and user defined options for PETSc @@ -568,26 +568,26 @@ subroutine formResidual(in, FandF_tau, & tensorField_real = 0.0_pReal do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) tensorField_real(1:3,1:3,i,j,k) = & - num%polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& - num%polarAlpha*matmul(F(1:3,1:3,i,j,k), & + num%beta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& + num%alpha*matmul(F(1:3,1:3,i,j,k), & math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3)) enddo; enddo; enddo !-------------------------------------------------------------------------------------------------- ! doing convolution in Fourier space call utilities_FFTtensorForward - call utilities_fourierGammaConvolution(params%rotation_BC%rotate(num%polarBeta*F_aim,active=.true.)) + call utilities_fourierGammaConvolution(params%rotation_BC%rotate(num%beta*F_aim,active=.true.)) call utilities_FFTtensorBackward !-------------------------------------------------------------------------------------------------- ! constructing residual - residual_F_tau = num%polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) + residual_F_tau = num%beta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response call utilities_constitutiveResponse(residual_F, & ! "residuum" gets field of first PK stress (to save memory) P_av,C_volAvg,C_minMaxAvg, & - F - residual_F_tau/num%polarBeta,params%timeinc,params%rotation_BC) + F - residual_F_tau/num%beta,params%timeinc,params%rotation_BC) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) !-------------------------------------------------------------------------------------------------- diff --git a/src/numerics.f90 b/src/numerics.f90 index 132c01f87..56bfc2492 100644 --- a/src/numerics.f90 +++ b/src/numerics.f90 @@ -19,7 +19,7 @@ module numerics implicit none private - class(tNode), pointer, public :: & + class(tNode), pointer, protected, public :: & numerics_root integer, protected, public :: & worldrank = 0, & !< MPI worldrank (/=0 for MPI simulations only) From ff858fd4c82188ca8b8e070030ac4885e9f62151 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Wed, 24 Jun 2020 20:13:09 +0200 Subject: [PATCH 40/43] [skip ci] corrected help string for "srepr" --- python/damask/util.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/util.py b/python/damask/util.py index a3fc71d3a..6c8fabf90 100644 --- a/python/damask/util.py +++ b/python/damask/util.py @@ -27,7 +27,7 @@ __all__=[ #################################################################################################### def srepr(arg,glue = '\n'): r""" - Join arguments as individual lines. + Join arguments with glue string. Parameters ---------- From 78bf8b0ab7f941d774f963d210db64e729ea7304 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Thu, 25 Jun 2020 11:25:39 +0200 Subject: [PATCH 41/43] bugfix: close file before returning --- src/IO.f90 | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/IO.f90 b/src/IO.f90 index 08f483a42..aae5f69cc 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -130,7 +130,10 @@ function IO_read(fileName) result(fileContent) status='old', position='rewind', action='read',iostat=myStat) if(myStat /= 0) call IO_error(100,ext_msg=trim(fileName)) allocate(character(len=fileLength)::fileContent) - if(fileLength==0) return + if(fileLength==0) then + close(fileUnit) + return + endif read(fileUnit,iostat=myStat) fileContent if(myStat /= 0) call IO_error(102,ext_msg=trim(fileName)) From 6e6e4dcdfd40c4d9b06c2611cfd5e518145bf651 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Thu, 25 Jun 2020 15:15:48 +0200 Subject: [PATCH 42/43] typo during resolve merge conflict --- src/mesh/discretization_mesh.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index dd7b76ff8..7dbd9350e 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -19,6 +19,7 @@ module discretization_mesh use numerics use FEsolving use FEM_quadrature + use YAML_types use prec implicit none From 6f9f494654ec6409c37dc2001708a6235007a1d8 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Fri, 26 Jun 2020 12:22:33 +0200 Subject: [PATCH 43/43] minor polishing --- src/CPFEM.f90 | 2 +- src/homogenization.f90 | 8 ++--- src/homogenization_mech_RGC.f90 | 56 ++++++++++++++++----------------- src/numerics.f90 | 8 ++--- 4 files changed, 37 insertions(+), 37 deletions(-) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 6647581e3..00afbbc01 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -46,7 +46,7 @@ module CPFEM type, private :: tNumerics integer :: & - iJacoStiffness !< frequency of stiffness update + iJacoStiffness !< frequency of stiffness update end type tNumerics type(tNumerics), private :: num diff --git a/src/homogenization.f90 b/src/homogenization.f90 index e055c6f06..a5ba56e9f 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -62,7 +62,7 @@ module homogenization module subroutine mech_RGC_init(num_homogMech) class(tNode), pointer, intent(in) :: & - num_homogMech + num_homogMech !< pointer to mechanical homogenization numerics data end subroutine mech_RGC_init @@ -139,8 +139,8 @@ subroutine homogenization_init num_homogMech, & num_homogGeneric - num_homog => numerics_root%get('homogenization',defaultVal=emptyDict) - num_homogMech => num_homog%get('mech',defaultVal=emptyDict) + num_homog => numerics_root%get('homogenization',defaultVal=emptyDict) + num_homogMech => num_homog%get('mech',defaultVal=emptyDict) num_homogGeneric => num_homog%get('generic',defaultVal=emptyDict) if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call mech_none_init @@ -169,7 +169,7 @@ subroutine homogenization_init if (debug_g < 1 .or. debug_g > homogenization_Ngrains(material_homogenizationAt(debug_e))) & call IO_error(602,ext_msg='constituent', el=debug_e, g=debug_g) - num%nMPstate = num_homogGeneric%get_asInt( 'nMPstate', defaultVal=10) + num%nMPstate = num_homogGeneric%get_asInt ('nMPstate', defaultVal=10) num%subStepMinHomog = num_homogGeneric%get_asFloat('subStepMin', defaultVal=1.0e-3_pReal) num%subStepSizeHomog = num_homogGeneric%get_asFloat('subStepSize', defaultVal=0.25_pReal) num%stepIncreaseHomog = num_homogGeneric%get_asFloat('stepIncrease', defaultVal=1.5_pReal) diff --git a/src/homogenization_mech_RGC.f90 b/src/homogenization_mech_RGC.f90 index 7ef73b130..90e08b71a 100644 --- a/src/homogenization_mech_RGC.f90 +++ b/src/homogenization_mech_RGC.f90 @@ -45,19 +45,19 @@ submodule(homogenization) homogenization_mech_RGC type :: tNumerics_RGC real(pReal) :: & - atol, & !< absolute tolerance of RGC residuum - rtol, & !< relative tolerance of RGC residuum - absMax, & !< absolute maximum of RGC residuum - relMax, & !< relative maximum of RGC residuum - pPert, & !< perturbation for computing RGC penalty tangent - xSmoo, & !< RGC penalty smoothing parameter (hyperbolic tangent) - viscPower, & !< power (sensitivity rate) of numerical viscosity in RGC scheme, Default 1.0e0: Newton viscosity (linear model) - viscModus, & !< stress modulus of RGC numerical viscosity, Default 0.0e0: No viscosity is applied - refRelaxRate, & !< reference relaxation rate in RGC viscosity - maxdRelax, & !< threshold of maximum relaxation vector increment (if exceed this then cutback) - maxVolDiscr, & !< threshold of maximum volume discrepancy allowed - volDiscrMod, & !< stiffness of RGC volume discrepancy (zero = without volume discrepancy constraint) - volDiscrPow !< powerlaw penalty for volume discrepancy + atol, & !< absolute tolerance of RGC residuum + rtol, & !< relative tolerance of RGC residuum + absMax, & !< absolute maximum of RGC residuum + relMax, & !< relative maximum of RGC residuum + pPert, & !< perturbation for computing RGC penalty tangent + xSmoo, & !< RGC penalty smoothing parameter (hyperbolic tangent) + viscPower, & !< power (sensitivity rate) of numerical viscosity in RGC scheme, Default 1.0e0: Newton viscosity (linear model) + viscModus, & !< stress modulus of RGC numerical viscosity, Default 0.0e0: No viscosity is applied + refRelaxRate, & !< reference relaxation rate in RGC viscosity + maxdRelax, & !< threshold of maximum relaxation vector increment (if exceed this then cutback) + maxVolDiscr, & !< threshold of maximum volume discrepancy allowed + volDiscrMod, & !< stiffness of RGC volume discrepancy (zero = without volume discrepancy constraint) + volDiscrPow !< powerlaw penalty for volume discrepancy end type tNumerics_RGC type(tparameters), dimension(:), allocatable :: & @@ -78,7 +78,7 @@ contains module subroutine mech_RGC_init(num_homogMech) class(tNode), pointer, intent(in) :: & - num_homogMech + num_homogMech !< pointer to mechanical homogenization numerics data integer :: & Ninstance, & @@ -87,7 +87,7 @@ module subroutine mech_RGC_init(num_homogMech) sizeState, nIntFaceTot class (tNode), pointer :: & - num_RGC + num_RGC ! pointer to RGC numerics data write(6,'(/,a)') ' <<<+- homogenization_'//HOMOGENIZATION_RGC_label//' init -+>>>'; flush(6) @@ -108,19 +108,19 @@ module subroutine mech_RGC_init(num_homogMech) num_RGC => num_homogMech%get('RGC',defaultVal=emptyDict) - num%atol = num_RGC%get_asFloat('atol', defaultVal=1.0e+4_pReal) - num%rtol = num_RGC%get_asFloat('rtol', defaultVal=1.0e-3_pReal) - num%absMax = num_RGC%get_asFloat('amax', defaultVal=1.0e+10_pReal) - num%relMax = num_RGC%get_asFloat('rmax', defaultVal=1.0e+2_pReal) - num%pPert = num_RGC%get_asFloat('perturbpenalty', defaultVal=1.0e-7_pReal) - num%xSmoo = num_RGC%get_asFloat('relvantmismatch', defaultVal=1.0e-5_pReal) - num%viscPower = num_RGC%get_asFloat('viscositypower', defaultVal=1.0e+0_pReal) - num%viscModus = num_RGC%get_asFloat('viscositymodulus', defaultVal=0.0e+0_pReal) - num%refRelaxRate = num_RGC%get_asFloat('refrelaxationrate',defaultVal=1.0e-3_pReal) - num%maxdRelax = num_RGC%get_asFloat('maxrelaxationrate',defaultVal=1.0e+0_pReal) - num%maxVolDiscr = num_RGC%get_asFloat('maxvoldiscrepancy',defaultVal=1.0e-5_pReal) - num%volDiscrMod = num_RGC%get_asFloat('voldiscrepancymod',defaultVal=1.0e+12_pReal) - num%volDiscrPow = num_RGC%get_asFloat('dicrepancypower', defaultVal=5.0_pReal) + num%atol = num_RGC%get_asFloat('atol', defaultVal=1.0e+4_pReal) + num%rtol = num_RGC%get_asFloat('rtol', defaultVal=1.0e-3_pReal) + num%absMax = num_RGC%get_asFloat('amax', defaultVal=1.0e+10_pReal) + num%relMax = num_RGC%get_asFloat('rmax', defaultVal=1.0e+2_pReal) + num%pPert = num_RGC%get_asFloat('perturbpenalty', defaultVal=1.0e-7_pReal) + num%xSmoo = num_RGC%get_asFloat('relvantmismatch', defaultVal=1.0e-5_pReal) + num%viscPower = num_RGC%get_asFloat('viscositypower', defaultVal=1.0e+0_pReal) + num%viscModus = num_RGC%get_asFloat('viscositymodulus', defaultVal=0.0e+0_pReal) + num%refRelaxRate = num_RGC%get_asFloat('refrelaxationrate', defaultVal=1.0e-3_pReal) + num%maxdRelax = num_RGC%get_asFloat('maxrelaxationrate', defaultVal=1.0e+0_pReal) + num%maxVolDiscr = num_RGC%get_asFloat('maxvoldiscrepancy', defaultVal=1.0e-5_pReal) + num%volDiscrMod = num_RGC%get_asFloat('voldiscrepancymod', defaultVal=1.0e+12_pReal) + num%volDiscrPow = num_RGC%get_asFloat('dicrepancypower', defaultVal=5.0_pReal) if (num%atol <= 0.0_pReal) call IO_error(301,ext_msg='absTol_RGC') diff --git a/src/numerics.f90 b/src/numerics.f90 index 56bfc2492..35436296c 100644 --- a/src/numerics.f90 +++ b/src/numerics.f90 @@ -20,12 +20,12 @@ module numerics private class(tNode), pointer, protected, public :: & - numerics_root + numerics_root !< root pointer storing the numerics YAML structure integer, protected, public :: & - worldrank = 0, & !< MPI worldrank (/=0 for MPI simulations only) - worldsize = 1 !< MPI worldsize (/=1 for MPI simulations only) + worldrank = 0, & !< MPI worldrank (/=0 for MPI simulations only) + worldsize = 1 !< MPI worldsize (/=1 for MPI simulations only) integer(4), protected, public :: & - DAMASK_NumThreadsInt = 0 !< value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive + DAMASK_NumThreadsInt = 0 !< value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive public :: numerics_init