Merge branch 'YAML-compatible-numerics' into YAML-compatible-debug

PRIVATE

@@ -1 +1 @@
-Subproject commit a22709e69a72cd1930385a77048894eea814a7fb
+Subproject commit a33c199e94fc618d5ff6237acef751175d813cac

VERSION

@@ -1 +1 @@
-v2.0.3-2650-g512e54a7
+v2.0.3-2692-g3d93a5ff

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]

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()
+env = damask.Environment()
-baseDir = damaskEnv.relPath('processing/')
+bin_dir = env.root_dir/Path('bin')
-binDir  = damaskEnv.relPath('bin/')
 
-if not os.path.isdir(binDir):
+if not bin_dir.exists():
-  os.mkdir(binDir)
+    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:
+    for the_file in the_dir.glob('*.py'):
-  theDir = os.path.abspath(os.path.join(baseDir,subDir))
+        src = the_dir/the_file
-
+        dst = bin_dir/Path(the_file.with_suffix('').name)
-  sys.stdout.write('\n'+binDir+' ->\n'+theDir+damask.util.deemph(' ...')+'\n')
+        if dst.is_file(): dst.unlink() # dst.unlink(True) for Python >3.8
-
+        dst.symlink_to(src)
-  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)
 
 
 sys.stdout.write('\npruning broken links...\n')
-
+for filename in bin_dir.glob('*'):
-brokenLinks = 0
+    if not filename.is_file():
-
+        filename.unlink
-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')

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])
@@ -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
+  import py_mentat
-  except:
-    damask.util.croak('no valid Mentat release found.')
-    sys.exit(-1)
 
   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
+    py_mentat.py_get_int("nnodes()")
-  except:
+  except py_mentat.InputError as err:
-    damask.util.croak('failed. try setting Tools/Python/"Run as Separate Process" & "Initiate".')
+    damask.util.croak('{}. Try Tools/Python/"Run as Separate Process" & "Initiate".'.format(err))
-    sys.exit()
+    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]))
+  py_mentat.py_send('*open_model "{}"'.format(filenames[0]))
-  except: damask.util.croak('lost connection on sending open command for "{}".'.format(filenames[0]))

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):

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

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 = '.'):
+    @property
-        """Return absolute path from path relative to DAMASK root."""
+    def options(self):
-        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):
         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)

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_modification = False
 
@@ -1056,14 +1057,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())
@@ -1190,7 +1194,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())
 
 
@@ -1266,7 +1270,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],
+            v.write('{}_inc{}'.format(self.fname.stem,inc[3:].zfill(N_digits)))
-                                            inc[3:].zfill(N_digits))
-
-            v.write(file_out)

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π].
+      with the angular ranges as [0,2π], [0,π], [0,2π].
-    - the rotation angle ω is limited to the interval [0, π].
+    - 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):
+    def _standardize(self):
-        """Inverse rotation/backward rotation."""
+        """Standardize (ensure positive real hemisphere)."""
-        return self.copy().inverse()
-
-
-    def standardize(self):
-        """In-place quaternion representation with positive real part."""
         self.quaternion[self.quaternion[...,0] < 0.0] *= -1
         return self
 
-    def standardized(self):
+    def inverse(self):
-        """Quaternion representation with positive real part."""
+        """In-place inverse rotation (backward rotation)."""
-        return self.copy().standardize()
+        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
+        Returns
-        ----------
+        -------
-        quaternion : bool, optional
+        q : numpy.ndarray of shape (...,4)
-            return quaternion as DAMASK object.
+            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 cubochoric 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
+            accept_homomorph = acceptHomomorph # for compatibility
-        qu = np.array(quaternion,dtype=float)
+        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 lattice basis vectors.
 
+        Parameters
+        ----------
+        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)
         if om.shape[:-3:-1] != (3,3):
             raise ValueError('Invalid shape.')
@@ -356,20 +453,43 @@ class Rotation:
         return Rotation(Rotation._om2qu(om))
 
     @staticmethod
-    def from_matrix(om):
+    def from_matrix(R):
+        """
+        Initialize from rotation matrix.
 
-        return Rotation.from_basis(om)
+        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(rodrigues,
+    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 +499,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 +526,30 @@ class Rotation:
         return Rotation(Rotation._ho2qu(ho))
 
     @staticmethod
-    def from_cubochoric(cubochoric,
+    def from_cubochoric(c,
-                       P = -1):
+                        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:
+        if np.abs(np.max(cu)) > np.pi**(2./3.) * 0.5+1e-9:
-            raise ValueError('Cubochoric coordinate outside of the cube: {} {} {}.'.format(*cu))
+            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 +602,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)
@@ -471,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.
 #
@@ -530,7 +674,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 +697,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 +708,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)))
                                    ])

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]
+        ext = Path(fname).suffix
-        if ext == '.vtk':
+        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))
+            raise ValueError('Given extension {} does not match default .{}'.format(ext,default_ext))
-        writer.SetFileName('{}.{}'.format(name,default_ext))
+        writer.SetFileName(str(Path(fname).with_suffix('.'+default_ext)))
         writer.SetCompressorTypeToZLib()
         writer.SetDataModeToBinary()
         writer.SetInputData(self.geom)

python/damask/solver/_marc.py

@@ -1,7 +1,7 @@
-import os
 import subprocess
 import shlex
 import string
+from pathlib import Path
 
 from .._environment import Environment
 
@@ -19,28 +19,24 @@ class Marc:
 
         """
         self.solver  = 'Marc'
-        try:
+        self.version = version
-            self.version = int(version)
-        except TypeError:
-            self.version = -1
 
-
+    @property
-#--------------------------
+    def library_path(self):
-    def libraryPath(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_dir() else None
 
 
-#--------------------------
+    @property
-    def toolsPath(self):
+    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 os.path.exists(path_tools) else ''
+        return path_tools if path_tools.is_dir() else None
 
 
 #--------------------------
@@ -53,21 +49,21 @@ class Marc:
                   ):
 
 
-        damaskEnv = Environment()
+        env = Environment()
 
-        user = os.path.join(damaskEnv.relPath('src'),'DAMASK_marc{}.{}'.format(self.version,'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 os.path.isfile(user):
+        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'
+        if compile: cmd += ' -u ' + str(user) + ' -save y'
-        else:       cmd += ' -prog ' + os.path.splitext(user)[0]
+        else:       cmd += ' -prog ' + str(user.with_suffix(''))
 
         print('job submission {} compilation: {}'.format('with' if compile else 'without',user))
         if logfile: log = open(logfile, 'w')

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
     ----------
@@ -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
+    stdout, stderr = [i for i in (process.communicate() if stream_in is None
-                                                   else process.communicate(streamIn.read().encode('utf-8')))]
+                             else process.communicate(stream_in.read().encode('utf-8')))]
     os.chdir(initialPath)
-    out   = out.decode('utf-8').replace('\x08','')
+    stdout = stdout.decode('utf-8').replace('\x08','')
-    error = error.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):

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')))

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,10 @@ class TestRotation:
             rot = Rotation.from_basis(om,False,reciprocal=reciprocal)
             assert np.isclose(np.linalg.det(rot.as_matrix()),1.0)
 
+    @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 +870,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 +946,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@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@(R@data))
 
     @pytest.mark.parametrize('data',[np.random.rand(4),
                                      np.random.rand(3,2),
@@ -956,3 +974,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))

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'))
+        vtr = VTK.from_file(os.path.join(tmp_path,'rectilinearGrid.vtr'))
-        assert(s == v.__repr__())
+        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'))
+        vtp = VTK.from_file(os.path.join(tmp_path,'polyData.vtp'))
-        assert(s == v.__repr__())
+        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'))
+        vtu = VTK.from_file(os.path.join(tmp_path,'unstructuredGrid.vtu'))
-        assert(s == v.__repr__())
+        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)

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==''

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,7 +92,10 @@ end subroutine CPFEM_initAll
 !--------------------------------------------------------------------------------------------------
 subroutine CPFEM_init
 
-  class(tNode), pointer :: debug_CPFEM
+  class(tNode), pointer :: &
+    num_commercialFEM, &
+    debug_CPFEM
+
   write(6,'(/,a)')   ' <<<+-  CPFEM init  -+>>>'
   flush(6)
 
@@ -94,6 +103,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')
+!------------------------------------------------------------------------------
+
   debug_CPFEM => debug_root%get('cpfem',defaultVal=emptyList)
   if(debug_CPFEM%contains('basic')) then
     write(6,'(a32,1x,6(i8,1x))')   'CPFEM_cs:              ', shape(CPFEM_cs)
@@ -127,23 +143,15 @@ 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, &
+                                                      i, j, k, l, m, n, ph, homog, mySource
-                                                      iJacoStiffness                                !< frequency of stiffness update
   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, &
     debug_CPFEM
 
-!------------------------------------------------------------------------------
-! 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)
 
@@ -185,7 +193,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 (debug_CPFEM%contains('extensive')) &

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
@@ -106,7 +105,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
 
@@ -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,17 @@ 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) then
+    close(fileUnit)
+    return
+  endif
+
   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
 
 

src/YAML_parse.f90

@@ -34,24 +34,23 @@ 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
+  character(len=*), intent(inout)             :: flow_string                                        !< YAML file in flow style
-  class(tDict), intent(in), optional, target  :: defaultVal
   class (tNode), pointer                      :: node
  
   class (tNode),    pointer       :: myVal
   character(len=pStringLen)       :: key
 
-  integer                         :: e, &                                                           !> end position of dictionary or list
+  integer                         :: e, &                                                           ! end position of dictionary or list
-                                     s, &                                                           !> start position of dictionary or list
+                                     s, &                                                           ! start position of dictionary or list
-                                     d                                                              !> position of key: value separator (':')
+                                     d                                                              ! position of key: value separator (':')
   
   flow_string = trim(adjustl(flow_string(:)))
-  if (len_trim(flow_string) == 0 .and. present(defaultVal)) then
+  if (len_trim(flow_string) == 0) then
-    node => defaultVal
+    node => emptyDict
     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))
@@ -97,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

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

src/constitutive_plastic_dislotwin.f90

@@ -833,7 +833,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')

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
@@ -198,7 +198,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')
@@ -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 :: &

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, &
+    phi, phiDot, dPhiDot_dPhi
-    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')
-  
   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 &

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
 

src/debug.f90

@@ -42,7 +42,7 @@ subroutine debug_init
   fileExists: if (fexist) then
     debug_input  = IO_read('debug.yaml') 
     debug_inFlow = to_flow(debug_input)
-    debug_root   => parse_flow(debug_inFlow,defaultVal=emptyDict)
+    debug_root   => parse_flow(debug_inFlow)
   endif fileExists
 
 end subroutine debug_init

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)
 

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
@@ -93,7 +98,6 @@ program DAMASK_grid
     quit
   class (tNode), pointer :: &
     num_grid, &
-    num_generic, &
     debug_grid                                                                                      ! pointer to grid debug options
 
 !--------------------------------------------------------------------------------------------------
@@ -115,19 +119,17 @@ program DAMASK_grid
 
 !-------------------------------------------------------------------------------------------------
 ! reading field paramters from numerics file and do sanity checks
-  num_generic => numerics_root%get('generic', defaultVal=emptyDict)
+  num_grid => numerics_root%get('grid', defaultVal=emptyDict)
-  stagItMax  = num_generic%get_asInt('maxStaggeredIter',defaultVal=10)
+  num%stagItMax  = num_grid%get_asInt('maxStaggeredIter',defaultVal=10)
-  maxCutBack = num_generic%get_asInt('maxCutBack',defaultVal=3)
+  num%maxCutBack = num_grid%get_asInt('maxCutBack',defaultVal=3)
 
-  if (stagItMax < 0)    call IO_error(301,ext_msg='maxStaggeredIter')
+  if (num%stagItMax < 0)    call IO_error(301,ext_msg='maxStaggeredIter')
-  if (maxCutBack < 0)   call IO_error(301,ext_msg='maxCutBack')
+  if (num%maxCutBack < 0)   call IO_error(301,ext_msg='maxCutBack')
 
 !--------------------------------------------------------------------------------------------------
 ! assign mechanics solver depending on selected type
  
   debug_grid => debug_root%get('grid',defaultVal=emptyList)
-  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
@@ -321,7 +323,6 @@ program DAMASK_grid
     loadCases = [loadCases,newLoadCase]                                                             ! load case is ok, append it
   enddo
 
-
 !--------------------------------------------------------------------------------------------------
 ! doing initialization depending on active solvers
   call spectral_Utilities_init
@@ -454,7 +455,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
@@ -473,7 +474,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

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_generic
   character(len=pStringLen) :: &
-    snes_type, &
+    snes_type
-    petsc_options
  
   write(6,'(/,a)') ' <<<+-  grid_spectral_damage init  -+>>>'
 
@@ -72,16 +84,27 @@ 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)
+  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)
-  petsc_options = num_generic%get_asString('petsc_options',defaultVal='')
+  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')
+  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
  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 +169,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, &
+  integer :: i, j, k, cell
-    itmax                                                                                           !< maximum number of iterations
   type(tSolutionState) :: solution
-  class(tNode), pointer :: &
-    num_generic
   PetscInt  :: devNull
   PetscReal :: phi_min, phi_max, stagNorm, solnNorm
   
   PetscErrorCode :: ierr
   SNESConvergedReason :: reason
 
-!-------------------------------------------------------------------
-! reading numerical parameter and do sanity check
-  num_generic => numerics_root%get('generic',defaultVal=emptyDict)
-  itmax = num_generic%get_asInt('itmax',defaultVal=250)
-  if (itmax <= 1)   call IO_error(301,ext_msg='itmax')
-
   solution%converged =.false.
  
 !--------------------------------------------------------------------------------------------------
@@ -175,7 +189,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 +199,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 +270,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, &
+  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')
 
   phi_current = x_scal 
 !--------------------------------------------------------------------------------------------------
@@ -299,8 +306,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) &
+  where(scalarField_real(1:grid(1),1:grid(2),1:grid3) < num%residualStiffness) &
-        scalarField_real(1:grid(1),1:grid(2),1:grid3) = residualStiffness
+        scalarField_real(1:grid(1),1:grid(2),1:grid3) = num%residualStiffness
  
 !--------------------------------------------------------------------------------------------------
 ! constructing residual

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,10 +114,9 @@ subroutine grid_mech_FEM_init
   PetscErrorCode :: ierr
   integer(HID_T) :: fileHandle, groupHandle
   character(len=pStringLen) :: &
-    fileName, &
+    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
@@ -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_generic => numerics_root%get('generic',defaultVal=emptyDict)
+  num_grid => numerics_root%get('grid',defaultVal=emptyDict)
-  petsc_options = num_generic%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,51 +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, &
+    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, &
-    num_generic
 
-!-----------------------------------------------------------------------------------
-! 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)
-
-  num_generic => numerics_root%get('generic',defaultVal=emptyDict)
-  itmin = num_generic%get_asInt('itmin',defaultVal=1)
-  itmax = num_generic%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)
+  divTol = max(maxval(abs(P_av))*num%eps_div_rtol   ,num%eps_div_atol)
-  BCTol  = max(maxval(abs(P_av))*eps_stress_rtol,eps_stress_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
@@ -501,24 +503,11 @@ 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_generic, &
     debug_grid                                                                                      ! pointer to grid  debug options
 
-!---------------------------------------------------------------------
-! debug pointer to grid
   debug_grid => debug_root%get('grid',defaultVal=emptyList)
 
-!----------------------------------------------------------------------
-! 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)
-  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)
@@ -529,7 +518,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 (debug_grid%contains('rotation')) &
       write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
               ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))

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
@@ -89,8 +99,7 @@ subroutine grid_mech_spectral_basic_init
   real(pReal), dimension(3,3) :: &
     temp33_Real = 0.0_pReal
   class (tNode), pointer :: &
-    num_grid, &
+    num_grid
-    num_generic
  
   PetscErrorCode :: ierr
   PetscScalar, pointer, dimension(:,:,:,:) :: &
@@ -99,8 +108,7 @@ subroutine grid_mech_spectral_basic_init
   integer(HID_T) :: fileHandle, groupHandle
   integer        :: fileUnit
   character(len=pStringLen) :: &
-    fileName, &
+    fileName
-    petsc_options
 
   write(6,'(/,a)') ' <<<+-  grid_mech_spectral_basic init  -+>>>'; flush(6)
 
@@ -111,18 +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_generic => numerics_root%get('generic',defaultVal=emptyDict)
-  petsc_options = num_generic%get_asString('petsc_options',defaultVal='')
-
   num_grid => numerics_root%get('grid',defaultVal=emptyDict)
-  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)
 
 !--------------------------------------------------------------------------------------------------
@@ -392,48 +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, &
+    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, &
-    num_generic
-!-----------------------------------------------------------------------------------
-! 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)
 
-  num_generic => numerics_root%get('generic',defaultVal=emptyDict)
+  divTol = max(maxval(abs(P_av))*num%eps_div_rtol   ,num%eps_div_atol)
-  itmin = num_generic%get_asInt('itmin',defaultVal=1)
+  BCTol  = max(maxval(abs(P_av))*num%eps_stress_rtol,num%eps_stress_atol)
-  itmax = num_generic%get_asInt('itmax',defaultVal=250)
 
-  if (eps_div_atol <= 0.0_pReal)     call IO_error(301,ext_msg='eps_div_atol')
+  if ((totalIter >= num%itmin .and. &
-  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. &
                             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
@@ -470,24 +462,11 @@ subroutine formResidual(in, F, &
     nfuncs
   PetscObject :: dummy
   PetscErrorCode :: ierr
-  integer :: &
-    itmin, &
-    itmax
   class(tNode), pointer :: &
-    num_generic, &
     debug_grid                                                                                      ! pointer to constitutive debug options
 
-!---------------------------------------------------------------------
-! debug pointer to grid
   debug_grid => debug_root%get('grid', defaultVal=emptyList)
 
-!----------------------------------------------------------------------
-! 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)
-  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)
@@ -497,7 +476,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 (debug_grid%contains('rotation')) &
       write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
               ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))

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) :: &
+      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) :: num                                                                            ! numerics parameters. Better name?
+  type(tNumerics), private :: num                                                                   ! numerics parameters. Better name?
 
 !--------------------------------------------------------------------------------------------------
 ! PETSc data
@@ -96,8 +111,7 @@ subroutine grid_mech_spectral_polarisation_init
   real(pReal), dimension(3,3) :: &
     temp33_Real = 0.0_pReal
   class (tNode), pointer :: &
-    num_grid, &
+    num_grid
-    num_generic
 
   PetscErrorCode :: ierr
   PetscScalar, pointer, dimension(:,:,:,:) :: &
@@ -108,8 +122,7 @@ subroutine grid_mech_spectral_polarisation_init
   integer(HID_T) :: fileHandle, groupHandle
   integer        :: fileUnit
   character(len=pStringLen) :: &
-    fileName, &
+    fileName
-    petsc_options
 
   write(6,'(/,a)') ' <<<+-  grid_mech_spectral_polarisation init  -+>>>'; flush(6)
 
@@ -118,17 +131,37 @@ 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)
-  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%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')
+  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%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
   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)
 
 !--------------------------------------------------------------------------------------------------
@@ -439,58 +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, &
+    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, &
-    num_generic
 
-!-----------------------------------------------------------------------------------
+  curlTol    = max(maxval(abs(F_aim-math_I3))*num%eps_curl_rtol  ,num%eps_curl_atol)
-! reading numerical parameters and do sanity check 
+  divTol     = max(maxval(abs(P_av))         *num%eps_div_rtol   ,num%eps_div_atol)
-  num_grid => numerics_root%get('grid',defaultVal=emptyDict)
+  BCTol      = max(maxval(abs(P_av))         *num%eps_stress_rtol,num%eps_stress_atol)
-  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)
 
-  num_generic => numerics_root%get('generic',defaultVal=emptyDict)
+  if ((totalIter >= num%itmin .and. &
-  itmin = num_generic%get_asInt('itmin',defaultVal=1)
-  itmax = num_generic%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. &
                             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
@@ -533,34 +530,12 @@ subroutine formResidual(in, FandF_tau, &
   PetscObject :: dummy
   PetscErrorCode :: ierr
   class(tNode), pointer :: &
-    num_grid, &
-    num_generic, &
     debug_grid                                                                                      !  pointer to grid  debug options
-  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, &
+    i, j, k, e
-    itmin, itmax
 
-!--------------------------------------------------------------------------------------------------
-! debug pointer for grid
   debug_grid => debug_root%get('grid',defaultVal=emptyList)
 
-!--------------------------------------------------------------------------------------------------
-! 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)
-
-  num_generic => numerics_root%get('generic',defaultVal=emptyDict)
-  itmin = num_generic%get_asInt('itmin',defaultVal=1)
-  itmax = num_generic%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,&
@@ -579,11 +554,12 @@ subroutine formResidual(in, FandF_tau, &
   call SNESGetIterationNumber(snes,PETScIter,ierr);  CHKERRQ(ierr)
 
   if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1                                              ! new increment
+
 !--------------------------------------------------------------------------------------------------
 ! 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(debug_grid%contains('rotation')) &
       write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
               ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
@@ -597,26 +573,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) -&
+      num%beta*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%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(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 = 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/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)
 
 !--------------------------------------------------------------------------------------------------

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
@@ -63,9 +73,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
 
   write(6,'(/,a)') ' <<<+-  grid_thermal_spectral init  -+>>>'
 
@@ -73,15 +81,20 @@ 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_generic => numerics_root%get('generic',defaultVal=emptyDict)
+  num_grid => numerics_root%get('grid',defaultVal=emptyDict)
-  petsc_options = num_generic%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)
+
+  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)
  CHKERRQ(ierr)
- call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr)
+ call PETScOptionsInsertString(PETSC_NULL_OPTIONS,num%petsc_options,ierr)
  CHKERRQ(ierr)
  
 !--------------------------------------------------------------------------------------------------
@@ -147,7 +160,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 +169,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)
+  num_grid => numerics_root%get('grid',defaultVal=emptyDict)
-  itmax = num_generic%get_asInt('itmax',defaultVal=250)
+  itmax = num_grid%get_asInt('itmax',defaultVal=250)
   if (itmax <= 1)   call IO_error(301,ext_msg='itmax')
 
   solution%converged =.false.
@@ -182,7 +195,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 

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, &
@@ -190,11 +190,9 @@ subroutine spectral_utilities_init
     vecSize    = 3_C_INTPTR_T, &
     tensorSize = 9_C_INTPTR_T
   character(len=pStringLen) :: &
-    petsc_options, &
     PETSCDEBUG = ' -snes_view -snes_monitor '
   class (tNode) , pointer :: &
     num_grid, &
-    num_generic, &
     debug_grid                                                                                      ! pointer to grid  debug options
 
   write(6,'(/,a)') ' <<<+-  spectral_utilities init  -+>>>'
@@ -224,14 +222,14 @@ 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)
+  num_grid => numerics_root%get('grid',defaultVal=emptyDict)
-  petsc_options = num_generic%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
@@ -240,8 +238,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)

src/homogenization.f90

@@ -66,7 +66,7 @@ module homogenization
 
     module subroutine mech_RGC_init(num_homogMech, debug_homogenization)
       class(tNode), pointer, intent(in) :: &
-        num_homogMech, &
+        num_homogMech, &                                                                            !< pointer to mechanical homogenization numerics data
         debug_homogenization                                                                        !< pointer to debug options for homogenization
     end subroutine mech_RGC_init
 
@@ -186,7 +186,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)

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
+      atol, &                                                                                       !< absolute tolerance of RGC residuum
-     rtol, &                                                                                        !< relative tolerance of RGC residuum
+      rtol, &                                                                                       !< relative tolerance of RGC residuum
-     absMax, &                                                                                      !< absolute maximum of RGC residuum
+      absMax, &                                                                                     !< absolute maximum of RGC residuum
-     relMax, &                                                                                      !< relative maximum of RGC residuum
+      relMax, &                                                                                     !< relative maximum of RGC residuum
-     pPert, &                                                                                       !< perturbation for computing RGC penalty tangent
+      pPert, &                                                                                      !< perturbation for computing RGC penalty tangent
-     xSmoo, &                                                                                       !< RGC penalty smoothing parameter (hyperbolic 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)
+      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
+      viscModus, &                                                                                  !< stress modulus of RGC numerical viscosity, Default 0.0e0: No viscosity is applied
-     refRelaxRate, &                                                                                !< reference relaxation rate in RGC viscosity
+      refRelaxRate, &                                                                               !< reference relaxation rate in RGC viscosity
-     maxdRelax, &                                                                                   !< threshold of maximum relaxation vector increment (if exceed this then cutback)
+      maxdRelax, &                                                                                  !< threshold of maximum relaxation vector increment (if exceed this then cutback)
-     maxVolDiscr, &                                                                                 !< threshold of maximum volume discrepancy allowed
+      maxVolDiscr, &                                                                                !< threshold of maximum volume discrepancy allowed
-     volDiscrMod, &                                                                                 !< stiffness of RGC volume discrepancy (zero = without volume discrepancy constraint)
+      volDiscrMod, &                                                                                !< stiffness of RGC volume discrepancy (zero = without volume discrepancy constraint)
-     volDiscrPow                                                                                    !< powerlaw penalty for volume discrepancy
+      volDiscrPow                                                                                   !< powerlaw penalty for volume discrepancy
   end type tNumerics_RGC
 
   type(tparameters),          dimension(:), allocatable :: &
@@ -78,8 +78,9 @@ contains
 module subroutine mech_RGC_init(num_homogMech,debug_homogenization)
 
   class(tNode), pointer, intent(in) :: &
-    num_homogMech, &
+    num_homogMech, &                                                                                !< pointer to mechanical homogenization numerics data
     debug_homogenization                                                                            !< pointer to debug options for homogenization
+
   integer :: &
     Ninstance, &
     h, &
@@ -89,7 +90,7 @@ module subroutine mech_RGC_init(num_homogMech,debug_homogenization)
     debug_i 
   
   class (tNode), pointer :: &
-    num_RGC
+    num_RGC                                                                                         ! pointer to RGC numerics data
 
   write(6,'(/,a)') ' <<<+-  homogenization_'//HOMOGENIZATION_RGC_label//' init  -+>>>'; flush(6)
 
@@ -110,19 +111,19 @@ module subroutine mech_RGC_init(num_homogMech,debug_homogenization)
  
   num_RGC => num_homogMech%get('RGC',defaultVal=emptyDict)
   
-  num%atol         =  num_RGC%get_asFloat('atol',             defaultVal=1.0e+4_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%rtol         =  num_RGC%get_asFloat('rtol',              defaultVal=1.0e-3_pReal)
-  num%absMax       =  num_RGC%get_asFloat('amax',             defaultVal=1.0e+10_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%relMax       =  num_RGC%get_asFloat('rmax',              defaultVal=1.0e+2_pReal)
-  num%pPert        =  num_RGC%get_asFloat('perturbpenalty',   defaultVal=1.0e-7_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%xSmoo        =  num_RGC%get_asFloat('relvantmismatch',   defaultVal=1.0e-5_pReal)
-  num%viscPower    =  num_RGC%get_asFloat('viscositypower',   defaultVal=1.0e+0_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%viscModus    =  num_RGC%get_asFloat('viscositymodulus',  defaultVal=0.0e+0_pReal)
-  num%refRelaxRate =  num_RGC%get_asFloat('refrelaxationrate',defaultVal=1.0e-3_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%maxdRelax    =  num_RGC%get_asFloat('maxrelaxationrate', defaultVal=1.0e+0_pReal)
-  num%maxVolDiscr  =  num_RGC%get_asFloat('maxvoldiscrepancy',defaultVal=1.0e-5_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%volDiscrMod  =  num_RGC%get_asFloat('voldiscrepancymod', defaultVal=1.0e+12_pReal)
-  num%volDiscrPow  =  num_RGC%get_asFloat('dicrepancypower',  defaultVal=5.0_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')

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_elem, &                                                                         !< DAMASK element ID for Marc element ID
-    mesh_FEM2DAMASK_node                                                                             !< DAMASK node ID for Marc node ID
+    mesh_FEM2DAMASK_node                                                                            !< DAMASK node ID for Marc node ID
 
   public :: &
     discretization_marc_init

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, &
+    randSize
-    randomSeed                                                                                      !< fixed seeding for pseudo-random number generator, Default 0: use random seed
   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
+  if (num%randomSeed > 0) then
-    randInit = randomSeed
+    randInit = num%randomSeed
   else
     call random_seed()
     call random_seed(get = randInit)

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,14 +56,21 @@ 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_generic
+    num_mesh
   character(len=pStringLen), dimension(:), allocatable :: fileContent
   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
@@ -80,12 +86,12 @@ program DAMASK_mesh
 
 !--------------------------------------------------------------------- 
 ! reading field information from numerics file and do sanity checks
-  num_generic => numerics_root%get('generic', defaultVal=emptyDict)
+  num_mesh => numerics_root%get('mesh', defaultVal=emptyDict)
-  stagItMax  = num_generic%get_asInt('maxStaggeredIter',defaultVal=10)
+  num%stagItMax  = num_mesh%get_asInt('maxStaggeredIter',defaultVal=10)
-  maxCutBack = num_generic%get_asInt('maxCutBack',defaultVal=3)
+  num%maxCutBack = num_mesh%get_asInt('maxCutBack',defaultVal=3)
 
-  if (stagItMax < 0)    call IO_error(301,ext_msg='maxStaggeredIter')
+  if (num%stagItMax < 0)    call IO_error(301,ext_msg='maxStaggeredIter')
-  if (maxCutBack < 0)   call IO_error(301,ext_msg='maxCutBack')
+  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

src/mesh/FEM_utilities.f90

@@ -104,21 +104,20 @@ subroutine FEM_utilities_init
   character(len=pStringLen) :: petsc_optionsOrder
   class(tNode), pointer :: &
     num_mesh, &
-    num_generic, &
     debug_mesh                                                                                      ! pointer to mesh debug options
   integer :: structOrder                                                                            !< order of displacement shape functions
+  character(len=:), allocatable :: &
+    petsc_options
   character(len=pStringLen) :: &
-    petsc_options, &
     PETSCDEBUG = ' -snes_view -snes_monitor '
+
   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='')
-  num_generic => numerics_root%get('generic',defaultVal=emptyDict)
-  petsc_options = num_generic%get_asString('petsc_options', defaultVal='')
 
 !--------------------------------------------------------------------------------------------------
 ! set debugging parameters
@@ -141,7 +140,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)

src/mesh/discretization_mesh.f90

@@ -8,54 +8,51 @@ module discretization_mesh
 #include <petsc/finclude/petscdmplex.h>
 #include <petsc/finclude/petscis.h>
 #include <petsc/finclude/petscdmda.h>
- use PETScdmplex
+  use PETScdmplex
- use PETScdmda
+  use PETScdmda
- use PETScis
+  use PETScis
 
- use DAMASK_interface
+  use DAMASK_interface
- use IO
+  use IO
- use debug
+  use debug
- use discretization
+  use discretization
- use numerics
+  use numerics
- use FEsolving
+  use FEsolving
- use FEM_quadrature
+  use FEM_quadrature
- use prec
+  use YAML_types
- use YAML_types
+  use prec
- 
- implicit none
- private
- 
- integer, public, protected :: &
-   mesh_Nboundaries, &
-   mesh_NcpElemsGlobal
 
- integer :: &
+  implicit none
-   mesh_NcpElems                                                                                     !< total number of CP elements in mesh
+  private
+
+  integer, public, protected :: &
+    mesh_Nboundaries, &
+    mesh_NcpElemsGlobal
+
+  integer :: &
+    mesh_NcpElems                                                                                   !< total number of CP elements in mesh
 
 !!!! BEGIN DEPRECATED !!!!!
- integer, public, protected :: &
+  integer, public, protected :: &
-   mesh_maxNips                                                                                     !< max number of IPs in any CP element
+    mesh_maxNips                                                                                    !< max number of IPs in any CP element
 !!!! BEGIN DEPRECATED !!!!!
 
- integer, dimension(:,:), allocatable :: &
+  real(pReal), dimension(:,:), allocatable :: &
-   mesh_element !DEPRECATED
+    mesh_ipVolume, &                                                                                !< volume associated with IP (initially!)
+    mesh_node0                                                                                      !< node x,y,z coordinates (initially!)
 
- real(pReal), dimension(:,:), allocatable :: &
+  real(pReal), dimension(:,:,:), allocatable :: &
-   mesh_ipVolume, &                                                                                 !< volume associated with IP (initially!)
+    mesh_ipCoordinates                                                                              !< IP x,y,z coordinates (after deformation!)
-   mesh_node0                                                                                       !< node x,y,z coordinates (initially!)
- 
- real(pReal), dimension(:,:,:), allocatable :: &
-   mesh_ipCoordinates                                                                               !< IP x,y,z coordinates (after deformation!)
 
- DM, public :: geomMesh
+  DM, public :: geomMesh
- 
- PetscInt, dimension(:), allocatable, public, protected :: &
-   mesh_boundaries
 
- public :: &
+  PetscInt, dimension(:), allocatable, public, protected :: &
-   discretization_mesh_init, &
+    mesh_boundaries
-   mesh_FEM_build_ipVolumes, &
+
-   mesh_FEM_build_ipCoordinates
+  public :: &
+    discretization_mesh_init, &
+    mesh_FEM_build_ipVolumes, &
+    mesh_FEM_build_ipCoordinates
 
 contains
 
@@ -68,24 +65,21 @@ 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, &
     debug_e, debug_i
-  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 :: &
+    homogenizationAt, &
+    microstructureAt
   class(tNode), pointer :: &
     num_mesh
   integer :: integrationOrder                                                                       !< order of quadrature rule required
@@ -104,19 +98,15 @@ subroutine discretization_mesh_init(restart)
   debug_i = debug_root%get_asInt('integrationpoint',defaultVal=1)
 
  
-  ! 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)
@@ -136,7 +126,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
@@ -159,48 +149,42 @@ 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 = FEM_nQuadrature(dimPlex,integrationOrder)
-  mesh_maxNips = FE_Nips(1)
+
-  
-  write(6,*) 'mesh_maxNips',mesh_maxNips
   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( 3,j) = 1                                                                      ! homogenization
-    call DMGetLabelValue(geomMesh,'material',j-1,mesh_element(4,j),ierr)
-    CHKERRQ(ierr)
-  end do 
 
-  if (debug_e < 1 .or. debug_e > mesh_NcpElems) &
+  allocate(microstructureAt(mesh_NcpElems))
-    call IO_error(602,ext_msg='element')                                                        ! selected element does not exist
+  allocate(homogenizationAt(mesh_NcpElems),source=1)
-  if (debug_i < 1 .or. debug_i > FE_Nips(FE_geomtype(mesh_element(2,debug_e)))) &
+  do j = 1, mesh_NcpElems
-    call IO_error(602,ext_msg='IP')                                                             ! selected element does not have requested IP
+    call DMGetLabelValue(geomMesh,'material',j-1,microstructureAt(j),ierr)
-  
+    CHKERRQ(ierr)
-  FEsolving_execElem = [ 1,mesh_NcpElems ]                                                      ! parallel loop bounds set to comprise all DAMASK elements
+  end do
-  FEsolving_execIP   = [1,FE_Nips(FE_geomtype(mesh_element(2,1)))]
+
-  
+  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]
+
   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)
- 
+
 end subroutine discretization_mesh_init
 
 
@@ -208,23 +192,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
 
@@ -236,20 +220,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
@@ -263,7 +247,7 @@ subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints)
       enddo
       qOffset = qOffset + dimPlex
     enddo
-  enddo  
+  enddo
 
 end subroutine mesh_FEM_build_ipCoordinates
 

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
@@ -96,22 +108,22 @@ subroutine FEM_mech_init(fieldBC)
   PetscErrorCode                         :: ierr
 
   class(tNode), pointer :: &
-    num_mesh, &
+    num_mesh
-    num_generic
-  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_generic => numerics_root%get('generic',defaultVal=emptyDict)
+  if (num%itmax <= 1)                       call IO_error(301,ext_msg='itmax')
-  itmax = num_generic%get_asInt('itmax',defaultVal=250)
+  if (num%eps_struct_rtol <= 0.0_pReal)     call IO_error(301,ext_msg='eps_struct_rtol')
-  if (itmax <= 1) call IO_error(301,ext_msg='itmax')
+  if (num%eps_struct_atol <= 0.0_pReal)     call IO_error(301,ext_msg='eps_struct_atol')
 
 !--------------------------------------------------------------------------------------------------
 ! Setup FEM mech mesh
@@ -120,9 +132,9 @@ subroutine FEM_mech_init(fieldBC)
 
 !--------------------------------------------------------------------------------------------------
 ! Setup FEM mech discretization
-  qPoints  = FEM_quadrature_points( dimPlex,integrationOrder)%p
+  qPoints  = FEM_quadrature_points( dimPlex,num%integrationOrder)%p
-  qWeights = FEM_quadrature_weights(dimPlex,integrationOrder)%p
+  qWeights = FEM_quadrature_weights(dimPlex,num%integrationOrder)%p
-  nQuadrature = FEM_nQuadrature(    dimPlex,integrationOrder)
+  nQuadrature = FEM_nQuadrature(    dimPlex,num%integrationOrder)
   qPointsP  => qPoints
   qWeightsP => qWeights
   call PetscQuadratureCreate(PETSC_COMM_SELF,mechQuad,ierr); CHKERRQ(ierr)
@@ -131,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
@@ -222,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)
 
@@ -282,22 +294,10 @@ type(tSolutionState) function FEM_mech_solution( &
   character(len=*), intent(in) :: &
     incInfoIn
 
-!--------------------------------------------------------------------------------------------------
-  integer :: &
-    itmax
-  class(tNode), pointer :: &
-    num_generic
 
   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)
-  if (itmax <= 1) call IO_error(301,ext_msg='itmax')
-!-------------------------------------------------------------------------------------------------
-
   incInfo = incInfoIn
   FEM_mech_solution%converged =.false.
 !--------------------------------------------------------------------------------------------------
@@ -312,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)
@@ -350,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))
@@ -410,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) = &
@@ -499,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))
@@ -561,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))
@@ -578,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), &
@@ -688,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

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
@@ -18,13 +19,13 @@ module numerics
   implicit none
   private
   
-  class(tNode), pointer, public :: &
+  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)
+    worldrank                  =  0, &                                                              !< MPI worldrank (/=0 for MPI simulations only)
-    worldsize                  =  1                                                                  !< MPI worldsize (/=1 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
 
@@ -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)
+    numerics_root =>  parse_flow(numerics_inFlow)
-  else fileExists
+  endif
-    write(6,'(a,/)') ' using standard values'
-    flush(6)
-  endif fileExists
 
 !--------------------------------------------------------------------------------------------------
 ! openMP parameter

src/quaternions.f90

@@ -101,6 +101,8 @@ module quaternions
     assignment(=), &
     conjg, aimag, &
     log, exp, &
+    abs, dot_product, &
+    inverse, &
     real
 
 contains