diff --git a/PRIVATE b/PRIVATE index 25ce39dd0..3b498f5cb 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 25ce39dd0f5bf49dc5c2bec20767d93d2b76d353 +Subproject commit 3b498f5cb3c50e669588106de1b4cdc4c03ffff1 diff --git a/VERSION b/VERSION index f11394f7c..8ec654de6 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v3.0.0-alpha-381-g63f2419e9 +v3.0.0-alpha-406-gca6f7f7a3 diff --git a/examples/SpectralMethod/Polycrystal/material.yaml b/examples/SpectralMethod/Polycrystal/material.yaml index a3b75c5a3..394d16ed8 100644 --- a/examples/SpectralMethod/Polycrystal/material.yaml +++ b/examples/SpectralMethod/Polycrystal/material.yaml @@ -2,105 +2,105 @@ homogenization: SX: mech: {type: none} -microstructure: +material: - constituents: - fraction: 1.0 - orientation: [1.0, 0.0, 0.0, 0.0] + O: [1.0, 0.0, 0.0, 0.0] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.7936696712125002, -0.28765777461664166, -0.3436487135089419, 0.4113964260949434] + O: [0.7936696712125002, -0.28765777461664166, -0.3436487135089419, 0.4113964260949434] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.3986143167493579, -0.7014883552495493, 0.2154871765709027, 0.5500781677772945] + O: [0.3986143167493579, -0.7014883552495493, 0.2154871765709027, 0.5500781677772945] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.28645844315788244, -0.022571491243423537, -0.467933059311115, -0.8357456192708106] + O: [0.28645844315788244, -0.022571491243423537, -0.467933059311115, -0.8357456192708106] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.33012772942625784, -0.6781865350268957, 0.6494525351030648, 0.09638521992649676] + O: [0.33012772942625784, -0.6781865350268957, 0.6494525351030648, 0.09638521992649676] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.43596817439583935, -0.5982537129781701, 0.046599032277502436, 0.6707106499919265] + O: [0.43596817439583935, -0.5982537129781701, 0.046599032277502436, 0.6707106499919265] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.169734823419553, -0.699615227367322, -0.6059581215838098, -0.33844257746495854] + O: [0.169734823419553, -0.699615227367322, -0.6059581215838098, -0.33844257746495854] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.9698864809294915, 0.1729052643205874, -0.15948307917616958, 0.06315956884687175] + O: [0.9698864809294915, 0.1729052643205874, -0.15948307917616958, 0.06315956884687175] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.46205660912967883, 0.3105054068891252, -0.617849551030653, 0.555294529545738] + O: [0.46205660912967883, 0.3105054068891252, -0.617849551030653, 0.555294529545738] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.4512443497461787, -0.7636045534540555, -0.04739348426715133, -0.45939142396805815] + O: [0.4512443497461787, -0.7636045534540555, -0.04739348426715133, -0.45939142396805815] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.2161856212656443, -0.6581450184826598, -0.5498086209601588, 0.4667112513346289] + O: [0.2161856212656443, -0.6581450184826598, -0.5498086209601588, 0.4667112513346289] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.8753220715350803, -0.4561599367657419, -0.13298279533852678, -0.08969369719975541] + O: [0.8753220715350803, -0.4561599367657419, -0.13298279533852678, -0.08969369719975541] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.11908260752431069, 0.18266024809834172, -0.7144822594012615, -0.664807992845101] + O: [0.11908260752431069, 0.18266024809834172, -0.7144822594012615, -0.664807992845101] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.751104669484278, 0.5585633382623958, -0.34579336397009175, 0.06538900566860861] + O: [0.751104669484278, 0.5585633382623958, -0.34579336397009175, 0.06538900566860861] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.08740438971703973, 0.8991264096610437, -0.4156704205935976, 0.10559485570696363] + O: [0.08740438971703973, 0.8991264096610437, -0.4156704205935976, 0.10559485570696363] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.5584325870096193, 0.6016408353068798, -0.14280340445801173, 0.5529814994483859] + O: [0.5584325870096193, 0.6016408353068798, -0.14280340445801173, 0.5529814994483859] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.4052725440888093, 0.25253073423599154, 0.5693263597910454, -0.669215876471182] + O: [0.4052725440888093, 0.25253073423599154, 0.5693263597910454, -0.669215876471182] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.7570164606888676, 0.15265448024694664, -0.5998021466848317, 0.20942796551297105] + O: [0.7570164606888676, 0.15265448024694664, -0.5998021466848317, 0.20942796551297105] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.6987659297138081, -0.132172211261028, -0.19693254724422338, 0.6748883269678543] + O: [0.6987659297138081, -0.132172211261028, -0.19693254724422338, 0.6748883269678543] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.7729330445886478, 0.21682179052722322, -0.5207379472917645, 0.2905078484066341] + O: [0.7729330445886478, 0.21682179052722322, -0.5207379472917645, 0.2905078484066341] phase: Aluminum homogenization: SX diff --git a/processing/pre/geom_fromMinimalSurface.py b/processing/pre/geom_fromMinimalSurface.py index 83d1a3684..2b40940b7 100755 --- a/processing/pre/geom_fromMinimalSurface.py +++ b/processing/pre/geom_fromMinimalSurface.py @@ -4,8 +4,6 @@ import os import sys from optparse import OptionParser -import numpy as np - import damask @@ -13,14 +11,7 @@ scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) -minimal_surfaces = ['primitive','gyroid','diamond'] - -surface = { - 'primitive': lambda x,y,z: np.cos(x)+np.cos(y)+np.cos(z), - 'gyroid': lambda x,y,z: np.sin(x)*np.cos(y)+np.sin(y)*np.cos(z)+np.cos(x)*np.sin(z), - 'diamond': lambda x,y,z: np.cos(x-y)*np.cos(z)+np.sin(x+y)*np.sin(z), - } - +minimal_surfaces = list(damask.Geom._minimal_surface.keys()) # -------------------------------------------------------------------- # MAIN @@ -71,16 +62,8 @@ parser.set_defaults(type = minimal_surfaces[0], name = None if filename == [] else filename[0] damask.util.report(scriptName,name) -x,y,z = np.meshgrid(options.periods*2.0*np.pi*(np.arange(options.grid[0])+0.5)/options.grid[0], - options.periods*2.0*np.pi*(np.arange(options.grid[1])+0.5)/options.grid[1], - options.periods*2.0*np.pi*(np.arange(options.grid[2])+0.5)/options.grid[2], - indexing='xy',sparse=True) - -microstructure = np.where(options.threshold < surface[options.type](x,y,z), - options.microstructure[1],options.microstructure[0]) - -geom=damask.Geom(microstructure,options.size, - comments=[scriptID + ' ' + ' '.join(sys.argv[1:])]) +geom=damask.Geom.from_minimal_surface(options.grid,options.size,options.type,options.threshold, + options.periods,options.microstructure) damask.util.croak(geom) geom.save_ASCII(sys.stdout if name is None else name,compress=False) diff --git a/python/damask/_geom.py b/python/damask/_geom.py index 6a7bc3c0e..2ccbb1988 100644 --- a/python/damask/_geom.py +++ b/python/damask/_geom.py @@ -120,6 +120,29 @@ class Geom: return np.unique(self.material).size + @staticmethod + def load(fname): + """ + Read a VTK rectilinear grid. + + Parameters + ---------- + fname : str or or pathlib.Path + Geometry file to read. + Valid extension is .vtr, it will be appended if not given. + + """ + v = VTK.load(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr') + comments = v.get_comments() + grid = np.array(v.vtk_data.GetDimensions())-1 + bbox = np.array(v.vtk_data.GetBounds()).reshape(3,2).T + + return Geom(material = v.get('material').reshape(grid,order='F'), + size = bbox[1] - bbox[0], + origin = bbox[0], + comments=comments) + + @staticmethod def load_ASCII(fname): """ @@ -184,29 +207,6 @@ class Geom: return Geom(material.reshape(grid,order='F'),size,origin,comments) - @staticmethod - def load(fname): - """ - Read a VTK rectilinear grid. - - Parameters - ---------- - fname : str or or pathlib.Path - Geometry file to read. - Valid extension is .vtr, it will be appended if not given. - - """ - v = VTK.load(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr') - comments = v.get_comments() - grid = np.array(v.vtk_data.GetDimensions())-1 - bbox = np.array(v.vtk_data.GetBounds()).reshape(3,2).T - - return Geom(material = v.get('material').reshape(grid,order='F'), - size = bbox[1] - bbox[0], - origin = bbox[0], - comments=comments) - - @staticmethod def _find_closest_seed(seeds, weights, point): return np.argmin(np.sum((np.broadcast_to(point,(len(seeds),3))-seeds)**2,axis=1) - weights) @@ -292,6 +292,131 @@ class Geom: ) + _minimal_surface = \ + {'Schwarz P': lambda x,y,z: np.cos(x) + np.cos(y) + np.cos(z), + 'Double Primitive': lambda x,y,z: ( 0.5 * (np.cos(x)*np.cos(y) + np.cos(y)*np.cos(z) + np.cos(z)*np.cos(x)) + + 0.2 * (np.cos(2*x) + np.cos(2*y) + np.cos(2*z)) ), + 'Schwarz D': lambda x,y,z: ( np.sin(x)*np.sin(y)*np.sin(z) + + np.sin(x)*np.cos(y)*np.cos(z) + + np.cos(x)*np.cos(y)*np.sin(z) + + np.cos(x)*np.sin(y)*np.cos(z) ), + 'Complementary D': lambda x,y,z: ( np.cos(3*x+y)*np.cos(z) - np.sin(3*x-y)*np.sin(z) + np.cos(x+3*y)*np.cos(z) + + np.sin(x-3*y)*np.sin(z) + np.cos(x-y)*np.cos(3*z) - np.sin(x+y)*np.sin(3*z) ), + 'Double Diamond': lambda x,y,z: 0.5 * (np.sin(x)*np.sin(y) + + np.sin(y)*np.sin(z) + + np.sin(z)*np.sin(x) + + np.cos(x) * np.cos(y) * np.cos(z) ), + 'Dprime': lambda x,y,z: 0.5 * ( np.cos(x)*np.cos(y)*np.cos(z) + + np.cos(x)*np.sin(y)*np.sin(z) + + np.sin(x)*np.cos(y)*np.sin(z) + + np.sin(x)*np.sin(y)*np.cos(z) + - np.sin(2*x)*np.sin(2*y) + - np.sin(2*y)*np.sin(2*z) + - np.sin(2*z)*np.sin(2*x) ) - 0.2, + 'Gyroid': lambda x,y,z: np.cos(x)*np.sin(y) + np.cos(y)*np.sin(z) + np.cos(z)*np.sin(x), + 'Gprime': lambda x,y,z : ( np.sin(2*x)*np.cos(y)*np.sin(z) + + np.sin(2*y)*np.cos(z)*np.sin(x) + + np.sin(2*z)*np.cos(x)*np.sin(y) ) + 0.32, + 'Karcher K': lambda x,y,z: ( 0.3 * ( np.cos(x) + np.cos(y) + np.cos(z) + + np.cos(x)*np.cos(y) + np.cos(y)*np.cos(z) + np.cos(z)*np.cos(x) ) + - 0.4 * ( np.cos(2*x) + np.cos(2*y) + np.cos(2*z) ) ) + 0.2, + 'Lidinoid': lambda x,y,z: 0.5 * ( np.sin(2*x)*np.cos(y)*np.sin(z) + + np.sin(2*y)*np.cos(z)*np.sin(x) + + np.sin(2*z)*np.cos(x)*np.sin(y) + - np.cos(2*x)*np.cos(2*y) + - np.cos(2*y)*np.cos(2*z) + - np.cos(2*z)*np.cos(2*x) ) + 0.15, + 'Neovius': lambda x,y,z: ( 3 * (np.cos(x)+np.cos(y)+np.cos(z)) + + 4 * np.cos(x)*np.cos(y)*np.cos(z) ), + 'Fisher-Koch S': lambda x,y,z: ( np.cos(2*x)*np.sin( y)*np.cos( z) + + np.cos( x)*np.cos(2*y)*np.sin( z) + + np.sin( x)*np.cos( y)*np.cos(2*z) ), + } + + + @staticmethod + def from_minimal_surface(grid,size,surface,threshold=0.0,periods=1,materials=(1,2)): + """ + Generate geometry from definition of triply periodic minimal surface. + + Parameters + ---------- + grid : int numpy.ndarray of shape (3) + Number of grid points in x,y,z direction. + size : list or numpy.ndarray of shape (3) + Physical size of the geometry in meter. + surface : str + Type of the minimal surface. See notes for details. + threshold : float, optional. + Threshold of the minimal surface. Defaults to 0.0. + periods : integer, optional. + Number of periods per unit cell. Defaults to 1. + materials : (int, int), optional + Material IDs. Defaults to (1,2). + + Notes + ----- + The following triply-periodic minimal surfaces are implemented: + - Schwarz P + - Double Primitive + - Schwarz D + - Complementary D + - Double Diamond + - Dprime + - Gyroid + - Gprime + - Karcher K + - Lidinoid + - Neovius + - Fisher-Koch S + + References + ---------- + Surface curvature in triply-periodic minimal surface architectures as + a distinct design parameter in preparing advanced tissue engineering scaffolds + Sébastien B G Blanquer, Maike Werner, Markus Hannula, Shahriar Sharifi, + Guillaume P R Lajoinie, David Eglin, Jari Hyttinen, André A Poot, and Dirk W Grijpma + 10.1088/1758-5090/aa6553 + + Triply Periodic Bicontinuous Cubic Microdomain Morphologies by Symmetries + Meinhard Wohlgemuth, Nataliya Yufa, James Hoffman, and Edwin L. Thomas + 10.1021/ma0019499 + + Minisurf – A minimal surface generator for finite element modeling and additive manufacturing + Meng-Ting Hsieh, Lorenzo Valdevit + 10.1016/j.simpa.2020.100026 + + """ + x,y,z = np.meshgrid(periods*2.0*np.pi*(np.arange(grid[0])+0.5)/grid[0], + periods*2.0*np.pi*(np.arange(grid[1])+0.5)/grid[1], + periods*2.0*np.pi*(np.arange(grid[2])+0.5)/grid[2], + indexing='ij',sparse=True) + return Geom(material = np.where(threshold < Geom._minimal_surface[surface](x,y,z),materials[1],materials[0]), + size = size, + comments = util.execution_stamp('Geom','from_minimal_surface'), + ) + + + def save(self,fname,compress=True): + """ + Generates vtk rectilinear grid. + + Parameters + ---------- + fname : str, optional + Filename to write. If no file is given, a string is returned. + Valid extension is .vtr, it will be appended if not given. + compress : bool, optional + Compress with zlib algorithm. Defaults to True. + + """ + v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin) + v.add(self.material.flatten(order='F'),'material') + v.add_comments(self.comments) + + v.save(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr',parallel=False,compress=compress) + + def save_ASCII(self,fname,compress=None): """ Writes a geom file. @@ -364,26 +489,6 @@ class Geom: f.write(f'{reps} of {former}\n') - def save(self,fname,compress=True): - """ - Generates vtk rectilinear grid. - - Parameters - ---------- - fname : str, optional - Filename to write. If no file is given, a string is returned. - Valid extension is .vtr, it will be appended if not given. - compress : bool, optional - Compress with zlib algorithm. Defaults to True. - - """ - v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin) - v.add(self.material.flatten(order='F'),'material') - v.add_comments(self.comments) - - v.save(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr',parallel=False,compress=compress) - - def show(self): """Show on screen.""" v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin) diff --git a/python/tests/test_Geom.py b/python/tests/test_Geom.py index 21eb38cdd..c2e476f8f 100644 --- a/python/tests/test_Geom.py +++ b/python/tests/test_Geom.py @@ -360,3 +360,45 @@ class TestGeom: elif approach == 'Voronoi': geom = Geom.from_Voronoi_tessellation(grid,size,seeds, periodic=np.random.random()>0.5) assert np.all(geom.material == material) + + + @pytest.mark.parametrize('surface',['Schwarz P', + 'Double Primitive', + 'Schwarz D', + 'Complementary D', + 'Double Diamond', + 'Dprime', + 'Gyroid', + 'Gprime', + 'Karcher K', + 'Lidinoid', + 'Neovius', + 'Fisher-Koch S', + ]) + def test_minimal_surface_basic_properties(self,surface): + grid = np.random.randint(60,100,3) + size = np.ones(3)+np.random.rand(3) + threshold = 2*np.random.rand()-1. + periods = np.random.randint(2)+1 + materials = np.random.randint(0,40,2) + geom = Geom.from_minimal_surface(grid,size,surface,threshold,periods,materials) + assert set(geom.material.flatten()) | set(materials) == set(materials) \ + and (geom.size == size).all() and (geom.grid == grid).all() + + @pytest.mark.parametrize('surface,threshold',[('Schwarz P',0), + ('Double Primitive',-1./6.), + ('Schwarz D',0), + ('Complementary D',0), + ('Double Diamond',-0.133), + ('Dprime',-0.0395), + ('Gyroid',0), + ('Gprime',0.22913), + ('Karcher K',0.17045), + ('Lidinoid',0.14455), + ('Neovius',0), + ('Fisher-Koch S',0), + ]) + def test_minimal_surface_volume(self,surface,threshold): + grid = np.ones(3,dtype=int)*64 + geom = Geom.from_minimal_surface(grid,np.ones(3),surface,threshold) + assert np.isclose(np.count_nonzero(geom.material==1)/np.prod(geom.grid),.5,rtol=1e-3) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 05255e2a1..f93472e51 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -78,11 +78,11 @@ subroutine CPFEM_initAll call DAMASK_interface_init call prec_init call IO_init + call YAML_types_init + call YAML_parse_init call config_init call math_init call rotations_init - call YAML_types_init - call YAML_parse_init call HDF5_utilities_init call results_init(.false.) call discretization_marc_init diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index fed43ba78..994859758 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -48,11 +48,11 @@ subroutine CPFEM_initAll #ifdef Mesh call FEM_quadrature_init #endif + call YAML_types_init + call YAML_parse_init call config_init call math_init call rotations_init - call YAML_types_init - call YAML_parse_init call lattice_init call HDF5_utilities_init call results_init(restart=interface_restartInc>0) diff --git a/src/IO.f90 b/src/IO.f90 index f9e708ead..4a442cb06 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -427,20 +427,20 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) case (146) msg = 'number of values does not match' case (148) - msg = 'Nconstituents mismatch between homogenization and microstructure' + msg = 'Nconstituents mismatch between homogenization and material' !-------------------------------------------------------------------------------------------------- ! material error messages and related messages in mesh case (150) msg = 'index out of bounds' case (151) - msg = 'microstructure has no constituents' + msg = 'material has no constituents' case (153) msg = 'sum of phase fractions differs from 1' case (155) - msg = 'microstructure index out of bounds' + msg = 'material index out of bounds' case (180) - msg = 'missing/invalid microstructure definition via State Variable 2' + msg = 'missing/invalid material definition via State Variable 2' case (190) msg = 'unknown element type:' case (191) @@ -494,6 +494,10 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) msg = 'Unsupported feature' case (706) msg = 'Access by incorrect node type' + case (707) + msg = 'Abrupt end of file' + case (708) + msg = '--- expected after YAML file header' !------------------------------------------------------------------------------------------------- ! errors related to the grid solver @@ -522,7 +526,7 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) case (842) msg = 'incomplete information in grid mesh header' case (843) - msg = 'microstructure count mismatch' + msg = 'material count mismatch' case (844) msg = 'invalid VTR file' case (846) @@ -621,6 +625,9 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg) msg = 'polar decomposition failed' case (700) msg = 'unknown crystal symmetry' + case (709) + msg = 'read only the first document' + case (850) msg = 'max number of cut back exceeded, terminating' case default diff --git a/src/YAML_parse.f90 b/src/YAML_parse.f90 index d655bc2dc..cb5b726dc 100644 --- a/src/YAML_parse.f90 +++ b/src/YAML_parse.f90 @@ -227,6 +227,52 @@ logical function isKey(line) end function isKey +!-------------------------------------------------------------------------------------------------- +! @brief skip empty lines +! @details update start position in the block by skipping empty lines if present. +!-------------------------------------------------------------------------------------------------- +subroutine skip_empty_lines(blck,s_blck) + + character(len=*), intent(in) :: blck + integer, intent(inout) :: s_blck + + logical :: empty + + empty = .true. + do while(empty .and. len_trim(blck(s_blck:)) /= 0) + empty = len_trim(IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2))) == 0 + if(empty) s_blck = s_blck + index(blck(s_blck:),IO_EOL) + enddo + + +end subroutine skip_empty_lines + + +!-------------------------------------------------------------------------------------------------- +! @brief skip file header +! @details update start position in the block by skipping file header if present. +!-------------------------------------------------------------------------------------------------- +subroutine skip_file_header(blck,s_blck) + + character(len=*), intent(in) :: blck + integer, intent(inout) :: s_blck + + character(len=:), allocatable :: line + + line = IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2)) + if(index(adjustl(line),'%YAML') == 1) then + s_blck = s_blck + index(blck(s_blck:),IO_EOL) + call skip_empty_lines(blck,s_blck) + if(trim(IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2))) == '---') then + s_blck = s_blck + index(blck(s_blck:),IO_EOL) + else + call IO_error(708,ext_msg = line) + endif + endif + +end subroutine skip_file_header + + !-------------------------------------------------------------------------------------------------- ! @brief reads a line of YAML block which is already in flow style ! @details Dicts should be enlcosed within '{}' for it to be consistent with DAMASK YAML parser @@ -363,7 +409,9 @@ recursive subroutine lst(blck,flow,s_blck,s_flow,offset) do while (s_blck <= len_trim(blck)) e_blck = s_blck + index(blck(s_blck:),IO_EOL) - 2 line = IO_rmComment(blck(s_blck:e_blck)) - if (len_trim(line) == 0) then + if(trim(line) == '---' .or. trim(line) == '...') then + exit + elseif (len_trim(line) == 0) then s_blck = e_blck + 2 ! forward to next line cycle elseif(indentDepth(line,offset) > indent) then @@ -377,8 +425,10 @@ recursive subroutine lst(blck,flow,s_blck,s_flow,offset) else if(trim(adjustl(line)) == '-') then ! list item in next line s_blck = e_blck + 2 - e_blck = e_blck + index(blck(e_blck+2:),IO_EOL) + call skip_empty_lines(blck,s_blck) + e_blck = s_blck + index(blck(s_blck:),IO_EOL) - 2 line = IO_rmComment(blck(s_blck:e_blck)) + if(trim(line) == '---') call IO_error(707,ext_msg=line) if(indentDepth(line) < indent .or. indentDepth(line) == indent) & call IO_error(701,ext_msg=line) @@ -447,7 +497,9 @@ recursive subroutine dct(blck,flow,s_blck,s_flow,offset) do while (s_blck <= len_trim(blck)) e_blck = s_blck + index(blck(s_blck:),IO_EOL) - 2 line = IO_rmComment(blck(s_blck:e_blck)) - if (len_trim(line) == 0) then + if(trim(line) == '---' .or. trim(line) == '...') then + exit + elseif (len_trim(line) == 0) then s_blck = e_blck + 2 ! forward to next line cycle elseif(indentDepth(line,offset) < indent) then @@ -510,10 +562,12 @@ recursive subroutine decide(blck,flow,s_blck,s_flow,offset) character(len=:), allocatable :: line if(s_blck <= len(blck)) then + call skip_empty_lines(blck,s_blck) e_blck = s_blck + index(blck(s_blck:),IO_EOL) - 2 line = IO_rmComment(blck(s_blck:e_blck)) - - if(len_trim(line) == 0) then + if(trim(line) == '---' .or. trim(line) == '...') then + continue ! end parsing at this point but not stop the simulation + elseif(len_trim(line) == 0) then s_blck = e_blck +2 call decide(blck,flow,s_blck,s_flow,offset) elseif (isListItem(line)) then @@ -548,23 +602,30 @@ function to_flow(blck) character(len=:), allocatable :: to_flow character(len=*), intent(in) :: blck !< YAML mixed style + + character(len=:), allocatable :: line integer :: s_blck, & !< start position in blck s_flow, & !< start position in flow offset, & !< counts leading '- ' in nested lists end_line - if(isFlow(blck)) then - to_flow = trim(adjustl(blck)) - else - allocate(character(len=len(blck)*2)::to_flow) - ! move forward here (skip empty lines) and remove '----' if found - s_flow = 1 - s_blck = 1 - offset = 0 + + allocate(character(len=len(blck)*2)::to_flow) + s_flow = 1 + s_blck = 1 + offset = 0 + + if(len_trim(blck) /= 0) then + call skip_empty_lines(blck,s_blck) + call skip_file_header(blck,s_blck) + line = IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2)) + if(trim(line) == '---') s_blck = s_blck + index(blck(s_blck:),IO_EOL) call decide(blck,to_flow,s_blck,s_flow,offset) - to_flow = trim(to_flow(:s_flow-1)) endif - end_line = index(to_flow,IO_EOL) - if(end_line > 0) to_flow = to_flow(:end_line-1) + line = IO_rmComment(blck(s_blck:s_blck+index(blck(s_blck:),IO_EOL)-2)) + if(trim(line)== '---') call IO_warning(709,ext_msg=line) + to_flow = trim(to_flow(:s_flow-1)) + end_line = index(to_flow,IO_EOL) + if(end_line > 0) to_flow = to_flow(:end_line-1) end function to_flow @@ -636,6 +697,20 @@ subroutine selfTest if (.not. to_flow(block_dict_newline) == flow_dict) error stop 'to_flow' end block basic_dict + only_flow: block + character(len=*), parameter :: flow_dict = & + " {a: [b,c: {d: e}, f: g, e]}"//IO_EOL + character(len=*), parameter :: flow_list = & + " [a,b: c, d,e: {f: g}]"//IO_EOL + character(len=*), parameter :: flow_1 = & + "{a: [b, {c: {d: e}}, {f: g}, e]}" + character(len=*), parameter :: flow_2 = & + "[a, {b: c}, d, {e: {f: g}}]" + + if (.not. to_flow(flow_dict) == flow_1) error stop 'to_flow' + if (.not. to_flow(flow_list) == flow_2) error stop 'to_flow' + end block only_flow + basic_flow: block character(len=*), parameter :: flow_braces = & " source: [{param: 1}, {param: 2}, {param: 3}, {param: 4}]"//IO_EOL @@ -650,12 +725,21 @@ subroutine selfTest basic_mixed: block character(len=*), parameter :: block_flow = & + "%YAML 1.1"//IO_EOL//& + " "//IO_EOL//& + " "//IO_EOL//& + "---"//IO_EOL//& " aa:"//IO_EOL//& " - "//IO_EOL//& - " param_1: [a: b, c, {d: {e: [f: g, h]}}]"//IO_EOL//& + " "//IO_EOL//& + " "//IO_EOL//& + " param_1: [a: b, c, {d: {e: [f: g, h]}}]"//IO_EOL//& " - c: d"//IO_EOL//& " bb:"//IO_EOL//& - " - {param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}"//IO_EOL + " "//IO_EOL//& + " - "//IO_EOL//& + " {param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}"//IO_EOL//& + "..."//IO_EOL character(len=*), parameter :: mixed_flow = & "{aa: [{param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}, {c: d}], bb: [{param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}]}" diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 0bd816142..ad1c2d4a7 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -54,12 +54,10 @@ module crystallite ! crystallite_Lp0, & !< plastic velocitiy grad at start of FE inc crystallite_partionedLp0, & !< plastic velocity grad at start of homog inc - crystallite_subLp0,& !< plastic velocity grad at start of crystallite inc ! crystallite_Li, & !< current intermediate velocitiy grad (end of converged time step) crystallite_Li0, & !< intermediate velocitiy grad at start of FE inc crystallite_partionedLi0, & !< intermediate velocity grad at start of homog inc - crystallite_subLi0, & !< intermediate velocity grad at start of crystallite inc ! crystallite_S0, & !< 2nd Piola-Kirchhoff stress vector at start of FE inc crystallite_partionedS0 !< 2nd Piola-Kirchhoff stress vector at start of homog inc @@ -183,7 +181,6 @@ subroutine crystallite_init crystallite_Li,crystallite_Lp, & crystallite_subF,crystallite_subF0, & crystallite_subFp0,crystallite_subFi0, & - crystallite_subLi0,crystallite_subLp0, & source = crystallite_partionedF) allocate(crystallite_dPdF(3,3,3,3,cMax,iMax,eMax),source=0.0_pReal) @@ -326,34 +323,17 @@ function crystallite_stress() startIP, endIP, & s logical, dimension(homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: todo !ToDo: need to set some values to false for different Ngrains - todo = .false. + real(pReal), dimension(:,:,:,:,:), allocatable :: & + subLp0,& !< plastic velocity grad at start of crystallite inc + subLi0 !< intermediate velocity grad at start of crystallite inc + + + todo = .false. + + subLp0 = crystallite_partionedLp0 + subLi0 = crystallite_partionedLi0 + -#ifdef DEBUG - if (debugCrystallite%selective & - .and. FEsolving_execElem(1) <= debugCrystallite%element & - .and. debugCrystallite%element <= FEsolving_execElem(2)) then - print'(/,a,i8,1x,i2,1x,i3)', '<< CRYST stress >> boundary and initial values at el ip ipc ', & - debugCrystallite%element,debugCrystallite%ip, debugCrystallite%grain - print'(a,/,3(12x,3(f14.9,1x)/))', '<< CRYST stress >> F ', & - transpose(crystallite_partionedF(1:3,1:3,debugCrystallite%grain, & - debugCrystallite%ip,debugCrystallite%element)) - print'(a,/,3(12x,3(f14.9,1x)/))', '<< CRYST stress >> F0 ', & - transpose(crystallite_partionedF0(1:3,1:3,debugCrystallite%grain, & - debugCrystallite%ip,debugCrystallite%element)) - print'(a,/,3(12x,3(f14.9,1x)/))', '<< CRYST stress >> Fp0', & - transpose(crystallite_partionedFp0(1:3,1:3,debugCrystallite%grain, & - debugCrystallite%ip,debugCrystallite%element)) - print'(a,/,3(12x,3(f14.9,1x)/))', '<< CRYST stress >> Fi0', & - transpose(crystallite_partionedFi0(1:3,1:3,debugCrystallite%grain, & - debugCrystallite%ip,debugCrystallite%element)) - print'(a,/,3(12x,3(f14.9,1x)/))', '<< CRYST stress >> Lp0', & - transpose(crystallite_partionedLp0(1:3,1:3,debugCrystallite%grain, & - debugCrystallite%ip,debugCrystallite%element)) - print'(a,/,3(12x,3(f14.9,1x)/))', '<< CRYST stress >> Li0', & - transpose(crystallite_partionedLi0(1:3,1:3,debugCrystallite%grain, & - debugCrystallite%ip,debugCrystallite%element)) - endif -#endif !-------------------------------------------------------------------------------------------------- ! initialize to starting condition @@ -370,9 +350,7 @@ function crystallite_stress() sourceState(material_phaseAt(c,e))%p(s)%partionedState0(:,material_phaseMemberAt(c,i,e)) enddo crystallite_subFp0(1:3,1:3,c,i,e) = crystallite_partionedFp0(1:3,1:3,c,i,e) - crystallite_subLp0(1:3,1:3,c,i,e) = crystallite_partionedLp0(1:3,1:3,c,i,e) crystallite_subFi0(1:3,1:3,c,i,e) = crystallite_partionedFi0(1:3,1:3,c,i,e) - crystallite_subLi0(1:3,1:3,c,i,e) = crystallite_partionedLi0(1:3,1:3,c,i,e) crystallite_subF0(1:3,1:3,c,i,e) = crystallite_partionedF0(1:3,1:3,c,i,e) crystallite_subFrac(c,i,e) = 0.0_pReal crystallite_subStep(c,i,e) = 1.0_pReal/num%subStepSizeCryst @@ -415,8 +393,8 @@ function crystallite_stress() todo(c,i,e) = crystallite_subStep(c,i,e) > 0.0_pReal ! still time left to integrate on? if (todo(c,i,e)) then crystallite_subF0 (1:3,1:3,c,i,e) = crystallite_subF(1:3,1:3,c,i,e) - crystallite_subLp0(1:3,1:3,c,i,e) = crystallite_Lp (1:3,1:3,c,i,e) - crystallite_subLi0(1:3,1:3,c,i,e) = crystallite_Li (1:3,1:3,c,i,e) + subLp0(1:3,1:3,c,i,e) = crystallite_Lp (1:3,1:3,c,i,e) + subLi0(1:3,1:3,c,i,e) = crystallite_Li (1:3,1:3,c,i,e) crystallite_subFp0(1:3,1:3,c,i,e) = crystallite_Fp (1:3,1:3,c,i,e) crystallite_subFi0(1:3,1:3,c,i,e) = crystallite_Fi (1:3,1:3,c,i,e) plasticState( material_phaseAt(c,e))%subState0(:,material_phaseMemberAt(c,i,e)) & @@ -435,8 +413,8 @@ function crystallite_stress() crystallite_Fi (1:3,1:3,c,i,e) = crystallite_subFi0(1:3,1:3,c,i,e) crystallite_S (1:3,1:3,c,i,e) = crystallite_S0 (1:3,1:3,c,i,e) if (crystallite_subStep(c,i,e) < 1.0_pReal) then ! actual (not initial) cutback - crystallite_Lp (1:3,1:3,c,i,e) = crystallite_subLp0(1:3,1:3,c,i,e) - crystallite_Li (1:3,1:3,c,i,e) = crystallite_subLi0(1:3,1:3,c,i,e) + crystallite_Lp (1:3,1:3,c,i,e) = subLp0(1:3,1:3,c,i,e) + crystallite_Li (1:3,1:3,c,i,e) = subLi0(1:3,1:3,c,i,e) endif plasticState (material_phaseAt(c,e))%state( :,material_phaseMemberAt(c,i,e)) & = plasticState(material_phaseAt(c,e))%subState0(:,material_phaseMemberAt(c,i,e)) @@ -460,6 +438,7 @@ function crystallite_stress() math_inv33(crystallite_Fi(1:3,1:3,c,i,e))) crystallite_subdt(c,i,e) = crystallite_subStep(c,i,e) * crystallite_dt(c,i,e) crystallite_converged(c,i,e) = .false. + call integrateState(c,i,e) endif enddo @@ -467,9 +446,10 @@ function crystallite_stress() enddo elementLooping3 !$OMP END PARALLEL DO + call nonlocalConvergenceCheck + !-------------------------------------------------------------------------------------------------- ! integrate --- requires fully defined state array (basic + dependent state) - if (any(todo)) call integrateState(todo) ! TODO: unroll into proper elementloop to avoid N^2 for single point evaluation where(.not. crystallite_converged .and. crystallite_subStep > num%subStepMinCryst) & ! do not try non-converged but fully cutbacked any further todo = .true. ! TODO: again unroll this into proper elementloop to avoid N^2 for single point evaluation @@ -1125,14 +1105,14 @@ end function integrateStress !> @brief integrate stress, state with adaptive 1st order explicit Euler method !> using Fixed Point Iteration to adapt the stepsize !-------------------------------------------------------------------------------------------------- -subroutine integrateStateFPI(todo) +subroutine integrateStateFPI(g,i,e) - logical, dimension(:,:,:), intent(in) :: todo - integer :: & - NiterationState, & !< number of iterations in state loop + integer, intent(in) :: & e, & !< element index in element loop i, & !< integration point index in ip loop - g, & !< grain index in grain loop + g !< grain index in grain loop + integer :: & + NiterationState, & !< number of iterations in state loop p, & c, & s, & @@ -1147,101 +1127,88 @@ subroutine integrateStateFPI(todo) plastic_dotState real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState logical :: & - nonlocalBroken, broken + broken - nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE(size_pl,size_so,r,zeta,p,c,plastic_dotState,source_dotState,broken) - do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1),FEsolving_execIP(2) - do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) - p = material_phaseAt(g,e) - if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then + p = material_phaseAt(g,e) + c = material_phaseMemberAt(g,i,e) - c = material_phaseMemberAt(g,i,e) + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partionedF0, & + crystallite_Fi(1:3,1:3,g,i,e), & + crystallite_partionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + if(broken) return - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partionedF0, & - crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. - if(broken) cycle + size_pl = plasticState(p)%sizeDotState + plasticState(p)%state(1:size_pl,c) = plasticState(p)%subState0(1:size_pl,c) & + + plasticState(p)%dotState (1:size_pl,c) & + * crystallite_subdt(g,i,e) + plastic_dotState(1:size_pl,2) = 0.0_pReal + do s = 1, phase_Nsources(p) + size_so(s) = sourceState(p)%p(s)%sizeDotState + sourceState(p)%p(s)%state(1:size_so(s),c) = sourceState(p)%p(s)%subState0(1:size_so(s),c) & + + sourceState(p)%p(s)%dotState (1:size_so(s),c) & + * crystallite_subdt(g,i,e) + source_dotState(1:size_so(s),2,s) = 0.0_pReal + enddo - size_pl = plasticState(p)%sizeDotState - plasticState(p)%state(1:size_pl,c) = plasticState(p)%subState0(1:size_pl,c) & - + plasticState(p)%dotState (1:size_pl,c) & - * crystallite_subdt(g,i,e) - plastic_dotState(1:size_pl,2) = 0.0_pReal - do s = 1, phase_Nsources(p) - size_so(s) = sourceState(p)%p(s)%sizeDotState - sourceState(p)%p(s)%state(1:size_so(s),c) = sourceState(p)%p(s)%subState0(1:size_so(s),c) & - + sourceState(p)%p(s)%dotState (1:size_so(s),c) & - * crystallite_subdt(g,i,e) - source_dotState(1:size_so(s),2,s) = 0.0_pReal - enddo + iteration: do NiterationState = 1, num%nState - iteration: do NiterationState = 1, num%nState + if(nIterationState > 1) plastic_dotState(1:size_pl,2) = plastic_dotState(1:size_pl,1) + plastic_dotState(1:size_pl,1) = plasticState(p)%dotState(:,c) + do s = 1, phase_Nsources(p) + if(nIterationState > 1) source_dotState(1:size_so(s),2,s) = source_dotState(1:size_so(s),1,s) + source_dotState(1:size_so(s),1,s) = sourceState(p)%p(s)%dotState(:,c) + enddo - if(nIterationState > 1) plastic_dotState(1:size_pl,2) = plastic_dotState(1:size_pl,1) - plastic_dotState(1:size_pl,1) = plasticState(p)%dotState(:,c) - do s = 1, phase_Nsources(p) - if(nIterationState > 1) source_dotState(1:size_so(s),2,s) = source_dotState(1:size_so(s),1,s) - source_dotState(1:size_so(s),1,s) = sourceState(p)%p(s)%dotState(:,c) - enddo + broken = integrateStress(g,i,e) + if(broken) exit iteration - broken = integrateStress(g,i,e) - if(broken) exit iteration + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partionedF0, & + crystallite_Fi(1:3,1:3,g,i,e), & + crystallite_partionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + if(broken) exit iteration - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partionedF0, & - crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) exit iteration + zeta = damper(plasticState(p)%dotState(:,c),plastic_dotState(1:size_pl,1),& + plastic_dotState(1:size_pl,2)) + plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta & + + plastic_dotState(1:size_pl,1) * (1.0_pReal - zeta) + r(1:size_pl) = plasticState(p)%state (1:size_pl,c) & + - plasticState(p)%subState0(1:size_pl,c) & + - plasticState(p)%dotState (1:size_pl,c) * crystallite_subdt(g,i,e) + plasticState(p)%state(1:size_pl,c) = plasticState(p)%state(1:size_pl,c) & + - r(1:size_pl) + crystallite_converged(g,i,e) = converged(r(1:size_pl), & + plasticState(p)%state(1:size_pl,c), & + plasticState(p)%atol(1:size_pl)) + do s = 1, phase_Nsources(p) + zeta = damper(sourceState(p)%p(s)%dotState(:,c), & + source_dotState(1:size_so(s),1,s),& + source_dotState(1:size_so(s),2,s)) + sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) * zeta & + + source_dotState(1:size_so(s),1,s)* (1.0_pReal - zeta) + r(1:size_so(s)) = sourceState(p)%p(s)%state (1:size_so(s),c) & + - sourceState(p)%p(s)%subState0(1:size_so(s),c) & + - sourceState(p)%p(s)%dotState (1:size_so(s),c) * crystallite_subdt(g,i,e) + sourceState(p)%p(s)%state(1:size_so(s),c) = sourceState(p)%p(s)%state(1:size_so(s),c) & + - r(1:size_so(s)) + crystallite_converged(g,i,e) = & + crystallite_converged(g,i,e) .and. converged(r(1:size_so(s)), & + sourceState(p)%p(s)%state(1:size_so(s),c), & + sourceState(p)%p(s)%atol(1:size_so(s))) + enddo - zeta = damper(plasticState(p)%dotState(:,c),plastic_dotState(1:size_pl,1),& - plastic_dotState(1:size_pl,2)) - plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta & - + plastic_dotState(1:size_pl,1) * (1.0_pReal - zeta) - r(1:size_pl) = plasticState(p)%state (1:size_pl,c) & - - plasticState(p)%subState0(1:size_pl,c) & - - plasticState(p)%dotState (1:size_pl,c) * crystallite_subdt(g,i,e) - plasticState(p)%state(1:size_pl,c) = plasticState(p)%state(1:size_pl,c) & - - r(1:size_pl) - crystallite_converged(g,i,e) = converged(r(1:size_pl), & - plasticState(p)%state(1:size_pl,c), & - plasticState(p)%atol(1:size_pl)) - do s = 1, phase_Nsources(p) - zeta = damper(sourceState(p)%p(s)%dotState(:,c), & - source_dotState(1:size_so(s),1,s),& - source_dotState(1:size_so(s),2,s)) - sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) * zeta & - + source_dotState(1:size_so(s),1,s)* (1.0_pReal - zeta) - r(1:size_so(s)) = sourceState(p)%p(s)%state (1:size_so(s),c) & - - sourceState(p)%p(s)%subState0(1:size_so(s),c) & - - sourceState(p)%p(s)%dotState (1:size_so(s),c) * crystallite_subdt(g,i,e) - sourceState(p)%p(s)%state(1:size_so(s),c) = sourceState(p)%p(s)%state(1:size_so(s),c) & - - r(1:size_so(s)) - crystallite_converged(g,i,e) = & - crystallite_converged(g,i,e) .and. converged(r(1:size_so(s)), & - sourceState(p)%p(s)%state(1:size_so(s),c), & - sourceState(p)%p(s)%atol(1:size_so(s))) - enddo + if(crystallite_converged(g,i,e)) then + broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_Fe(1:3,1:3,g,i,e), & + crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) + exit iteration + endif - if(crystallite_converged(g,i,e)) then - broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_Fe(1:3,1:3,g,i,e), & - crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) - exit iteration - endif + enddo iteration - enddo iteration - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. - endif - enddo; enddo; enddo - !$OMP END PARALLEL DO - - if (nonlocalBroken) call nonlocalConvergenceCheck contains @@ -1271,64 +1238,48 @@ end subroutine integrateStateFPI !-------------------------------------------------------------------------------------------------- !> @brief integrate state with 1st order explicit Euler method !-------------------------------------------------------------------------------------------------- -subroutine integrateStateEuler(todo) +subroutine integrateStateEuler(g,i,e) - logical, dimension(:,:,:), intent(in) :: todo - - integer :: & + integer, intent(in) :: & e, & !< element index in element loop i, & !< integration point index in ip loop - g, & !< grain index in grain loop + g !< grain index in grain loop + integer :: & p, & c, & s, & sizeDotState logical :: & - nonlocalBroken, broken + broken - nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE (sizeDotState,p,c,broken) - do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1),FEsolving_execIP(2) - do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) - p = material_phaseAt(g,e) - if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then + p = material_phaseAt(g,e) + c = material_phaseMemberAt(g,i,e) - c = material_phaseMemberAt(g,i,e) + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partionedF0, & + crystallite_Fi(1:3,1:3,g,i,e), & + crystallite_partionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + if(broken) return - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partionedF0, & - crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. - if(broken) cycle + sizeDotState = plasticState(p)%sizeDotState + plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & + + plasticState(p)%dotState (1:sizeDotState,c) & + * crystallite_subdt(g,i,e) + do s = 1, phase_Nsources(p) + sizeDotState = sourceState(p)%p(s)%sizeDotState + sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & + + sourceState(p)%p(s)%dotState (1:sizeDotState,c) & + * crystallite_subdt(g,i,e) + enddo - sizeDotState = plasticState(p)%sizeDotState - plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & - + plasticState(p)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState - sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & - + sourceState(p)%p(s)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) - enddo + broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_Fe(1:3,1:3,g,i,e), & + crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) + if(broken) return - broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_Fe(1:3,1:3,g,i,e), & - crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. - if(broken) cycle - - broken = integrateStress(g,i,e) - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. - crystallite_converged(g,i,e) = .not. broken - endif - enddo; enddo; enddo - !$OMP END PARALLEL DO - - if (nonlocalBroken) call nonlocalConvergenceCheck + broken = integrateStress(g,i,e) + crystallite_converged(g,i,e) = .not. broken end subroutine integrateStateEuler @@ -1336,93 +1287,78 @@ end subroutine integrateStateEuler !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, state with 1st order Euler method with adaptive step size !-------------------------------------------------------------------------------------------------- -subroutine integrateStateAdaptiveEuler(todo) - - logical, dimension(:,:,:), intent(in) :: todo +subroutine integrateStateAdaptiveEuler(g,i,e) + integer, intent(in) :: & + e, & !< element index in element loop + i, & !< integration point index in ip loop + g !< grain index in grain loop integer :: & - e, & ! element index in element loop - i, & ! integration point index in ip loop - g, & ! grain index in grain loop p, & c, & s, & sizeDotState logical :: & - nonlocalBroken, broken + broken real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: residuum_plastic real(pReal), dimension(constitutive_source_maxSizeDotState,maxval(phase_Nsources)) :: residuum_source - nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,residuum_plastic,residuum_source,broken) - do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1),FEsolving_execIP(2) - do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) - broken = .false. - p = material_phaseAt(g,e) - if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then - c = material_phaseMemberAt(g,i,e) + p = material_phaseAt(g,e) + c = material_phaseMemberAt(g,i,e) - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partionedF0, & - crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) cycle + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partionedF0, & + crystallite_Fi(1:3,1:3,g,i,e), & + crystallite_partionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + if(broken) return - sizeDotState = plasticState(p)%sizeDotState + sizeDotState = plasticState(p)%sizeDotState - residuum_plastic(1:sizeDotState) = - plasticState(p)%dotstate(1:sizeDotState,c) * 0.5_pReal * crystallite_subdt(g,i,e) - plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & - + plasticState(p)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState + residuum_plastic(1:sizeDotState) = - plasticState(p)%dotstate(1:sizeDotState,c) * 0.5_pReal * crystallite_subdt(g,i,e) + plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & + + plasticState(p)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) + do s = 1, phase_Nsources(p) + sizeDotState = sourceState(p)%p(s)%sizeDotState - residuum_source(1:sizeDotState,s) = - sourceState(p)%p(s)%dotstate(1:sizeDotState,c) & - * 0.5_pReal * crystallite_subdt(g,i,e) - sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & - + sourceState(p)%p(s)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) - enddo + residuum_source(1:sizeDotState,s) = - sourceState(p)%p(s)%dotstate(1:sizeDotState,c) & + * 0.5_pReal * crystallite_subdt(g,i,e) + sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & + + sourceState(p)%p(s)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) + enddo - broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_Fe(1:3,1:3,g,i,e), & - crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) - if(broken) cycle + broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_Fe(1:3,1:3,g,i,e), & + crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) + if(broken) return - broken = integrateStress(g,i,e) - if(broken) cycle + broken = integrateStress(g,i,e) + if(broken) return - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partionedF0, & - crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) cycle + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partionedF0, & + crystallite_Fi(1:3,1:3,g,i,e), & + crystallite_partionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + if(broken) return - sizeDotState = plasticState(p)%sizeDotState - crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState) & - + 0.5_pReal * plasticState(p)%dotState(:,c) * crystallite_subdt(g,i,e), & - plasticState(p)%state(1:sizeDotState,c), & - plasticState(p)%atol(1:sizeDotState)) + sizeDotState = plasticState(p)%sizeDotState + crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState) & + + 0.5_pReal * plasticState(p)%dotState(:,c) * crystallite_subdt(g,i,e), & + plasticState(p)%state(1:sizeDotState,c), & + plasticState(p)%atol(1:sizeDotState)) - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState - crystallite_converged(g,i,e) = & - crystallite_converged(g,i,e) .and. converged(residuum_source(1:sizeDotState,s) & - + 0.5_pReal*sourceState(p)%p(s)%dotState(:,c)*crystallite_subdt(g,i,e), & - sourceState(p)%p(s)%state(1:sizeDotState,c), & - sourceState(p)%p(s)%atol(1:sizeDotState)) - enddo - - endif - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. - enddo; enddo; enddo - !$OMP END PARALLEL DO - - if (nonlocalBroken) call nonlocalConvergenceCheck + do s = 1, phase_Nsources(p) + sizeDotState = sourceState(p)%p(s)%sizeDotState + crystallite_converged(g,i,e) = & + crystallite_converged(g,i,e) .and. converged(residuum_source(1:sizeDotState,s) & + + 0.5_pReal*sourceState(p)%p(s)%dotState(:,c)*crystallite_subdt(g,i,e), & + sourceState(p)%p(s)%state(1:sizeDotState,c), & + sourceState(p)%p(s)%atol(1:sizeDotState)) + enddo end subroutine integrateStateAdaptiveEuler @@ -1430,9 +1366,9 @@ end subroutine integrateStateAdaptiveEuler !--------------------------------------------------------------------------------------------------- !> @brief Integrate state (including stress integration) with the classic Runge Kutta method !--------------------------------------------------------------------------------------------------- -subroutine integrateStateRK4(todo) +subroutine integrateStateRK4(g,i,e) - logical, dimension(:,:,:), intent(in) :: todo + integer, intent(in) :: g,i,e real(pReal), dimension(3,3), parameter :: & A = reshape([& @@ -1445,7 +1381,7 @@ subroutine integrateStateRK4(todo) real(pReal), dimension(4), parameter :: & B = [1.0_pReal/6.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/6.0_pReal] - call integrateStateRK(todo,A,B,C) + call integrateStateRK(g,i,e,A,B,C) end subroutine integrateStateRK4 @@ -1453,9 +1389,9 @@ end subroutine integrateStateRK4 !--------------------------------------------------------------------------------------------------- !> @brief Integrate state (including stress integration) with the Cash-Carp method !--------------------------------------------------------------------------------------------------- -subroutine integrateStateRKCK45(todo) +subroutine integrateStateRKCK45(g,i,e) - logical, dimension(:,:,:), intent(in) :: todo + integer, intent(in) :: g,i,e real(pReal), dimension(5,5), parameter :: & A = reshape([& @@ -1475,7 +1411,7 @@ subroutine integrateStateRKCK45(todo) [2825.0_pReal/27648.0_pReal, .0_pReal, 18575.0_pReal/48384.0_pReal,& 13525.0_pReal/55296.0_pReal, 277.0_pReal/14336.0_pReal, 1._pReal/4._pReal] - call integrateStateRK(todo,A,B,C,DB) + call integrateStateRK(g,i,e,A,B,C,DB) end subroutine integrateStateRKCK45 @@ -1484,18 +1420,18 @@ end subroutine integrateStateRKCK45 !> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an !! embedded explicit Runge-Kutta method !-------------------------------------------------------------------------------------------------- -subroutine integrateStateRK(todo,A,B,CC,DB) +subroutine integrateStateRK(g,i,e,A,B,CC,DB) - logical, dimension(:,:,:), intent(in) :: todo real(pReal), dimension(:,:), intent(in) :: A real(pReal), dimension(:), intent(in) :: B, CC real(pReal), dimension(:), intent(in), optional :: DB + integer, intent(in) :: & + e, & !< element index in element loop + i, & !< integration point index in ip loop + g !< grain index in grain loop integer :: & - e, & ! element index in element loop - i, & ! integration point index in ip loop - g, & ! grain index in grain loop stage, & ! stage index in integration stage loop n, & p, & @@ -1503,116 +1439,102 @@ subroutine integrateStateRK(todo,A,B,CC,DB) s, & sizeDotState logical :: & - nonlocalBroken, broken + broken real(pReal), dimension(constitutive_source_maxSizeDotState,size(B),maxval(phase_Nsources)) :: source_RKdotState real(pReal), dimension(constitutive_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState - nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,plastic_RKdotState,source_RKdotState,broken) - do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1),FEsolving_execIP(2) - do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) - broken = .false. - p = material_phaseAt(g,e) - if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then + p = material_phaseAt(g,e) + c = material_phaseMemberAt(g,i,e) - c = material_phaseMemberAt(g,i,e) + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partionedF0, & + crystallite_Fi(1:3,1:3,g,i,e), & + crystallite_partionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + if(broken) return - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partionedF0, & - crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) cycle + do stage = 1,size(A,1) + sizeDotState = plasticState(p)%sizeDotState + plastic_RKdotState(1:sizeDotState,stage) = plasticState(p)%dotState(:,c) + plasticState(p)%dotState(:,c) = A(1,stage) * plastic_RKdotState(1:sizeDotState,1) + do s = 1, phase_Nsources(p) + sizeDotState = sourceState(p)%p(s)%sizeDotState + source_RKdotState(1:sizeDotState,stage,s) = sourceState(p)%p(s)%dotState(:,c) + sourceState(p)%p(s)%dotState(:,c) = A(1,stage) * source_RKdotState(1:sizeDotState,1,s) + enddo - do stage = 1,size(A,1) - sizeDotState = plasticState(p)%sizeDotState - plastic_RKdotState(1:sizeDotState,stage) = plasticState(p)%dotState(:,c) - plasticState(p)%dotState(:,c) = A(1,stage) * plastic_RKdotState(1:sizeDotState,1) - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState - source_RKdotState(1:sizeDotState,stage,s) = sourceState(p)%p(s)%dotState(:,c) - sourceState(p)%p(s)%dotState(:,c) = A(1,stage) * source_RKdotState(1:sizeDotState,1,s) - enddo + do n = 2, stage + sizeDotState = plasticState(p)%sizeDotState + plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) & + + A(n,stage) * plastic_RKdotState(1:sizeDotState,n) + do s = 1, phase_Nsources(p) + sizeDotState = sourceState(p)%p(s)%sizeDotState + sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) & + + A(n,stage) * source_RKdotState(1:sizeDotState,n,s) + enddo + enddo - do n = 2, stage - sizeDotState = plasticState(p)%sizeDotState - plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) & - + A(n,stage) * plastic_RKdotState(1:sizeDotState,n) - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState - sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) & - + A(n,stage) * source_RKdotState(1:sizeDotState,n,s) - enddo - enddo - - sizeDotState = plasticState(p)%sizeDotState - plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & - + plasticState(p)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState - sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & - + sourceState(p)%p(s)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) - enddo - - broken = integrateStress(g,i,e,CC(stage)) - if(broken) exit - - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partionedF0, & - crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_partionedFp0, & - crystallite_subdt(g,i,e)*CC(stage), g,i,e,p,c) - if(broken) exit - - enddo - if(broken) cycle - - sizeDotState = plasticState(p)%sizeDotState - - plastic_RKdotState(1:sizeDotState,size(B)) = plasticState (p)%dotState(:,c) - plasticState(p)%dotState(:,c) = matmul(plastic_RKdotState(1:sizeDotState,1:size(B)),B) - plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & - + plasticState(p)%dotState (1:sizeDotState,c) & + sizeDotState = plasticState(p)%sizeDotState + plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & + + plasticState(p)%dotState (1:sizeDotState,c) & + * crystallite_subdt(g,i,e) + do s = 1, phase_Nsources(p) + sizeDotState = sourceState(p)%p(s)%sizeDotState + sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & + + sourceState(p)%p(s)%dotState (1:sizeDotState,c) & * crystallite_subdt(g,i,e) - if(present(DB)) & - broken = .not. converged( matmul(plastic_RKdotState(1:sizeDotState,1:size(DB)),DB) & - * crystallite_subdt(g,i,e), & - plasticState(p)%state(1:sizeDotState,c), & - plasticState(p)%atol(1:sizeDotState)) + enddo - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState + broken = integrateStress(g,i,e,CC(stage)) + if(broken) exit - source_RKdotState(1:sizeDotState,size(B),s) = sourceState(p)%p(s)%dotState(:,c) - sourceState(p)%p(s)%dotState(:,c) = matmul(source_RKdotState(1:sizeDotState,1:size(B),s),B) - sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & - + sourceState(p)%p(s)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) - if(present(DB)) & - broken = broken .or. .not. converged(matmul(source_RKdotState(1:sizeDotState,1:size(DB),s),DB) & - * crystallite_subdt(g,i,e), & - sourceState(p)%p(s)%state(1:sizeDotState,c), & - sourceState(p)%p(s)%atol(1:sizeDotState)) - enddo - if(broken) cycle + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partionedF0, & + crystallite_Fi(1:3,1:3,g,i,e), & + crystallite_partionedFp0, & + crystallite_subdt(g,i,e)*CC(stage), g,i,e,p,c) + if(broken) exit - broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_Fe(1:3,1:3,g,i,e), & - crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) - if(broken) cycle + enddo + if(broken) return - broken = integrateStress(g,i,e) - crystallite_converged(g,i,e) = .not. broken + sizeDotState = plasticState(p)%sizeDotState - endif - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. - enddo; enddo; enddo - !$OMP END PARALLEL DO + plastic_RKdotState(1:sizeDotState,size(B)) = plasticState (p)%dotState(:,c) + plasticState(p)%dotState(:,c) = matmul(plastic_RKdotState(1:sizeDotState,1:size(B)),B) + plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & + + plasticState(p)%dotState (1:sizeDotState,c) & + * crystallite_subdt(g,i,e) + if(present(DB)) & + broken = .not. converged( matmul(plastic_RKdotState(1:sizeDotState,1:size(DB)),DB) & + * crystallite_subdt(g,i,e), & + plasticState(p)%state(1:sizeDotState,c), & + plasticState(p)%atol(1:sizeDotState)) + + do s = 1, phase_Nsources(p) + sizeDotState = sourceState(p)%p(s)%sizeDotState + + source_RKdotState(1:sizeDotState,size(B),s) = sourceState(p)%p(s)%dotState(:,c) + sourceState(p)%p(s)%dotState(:,c) = matmul(source_RKdotState(1:sizeDotState,1:size(B),s),B) + sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & + + sourceState(p)%p(s)%dotState (1:sizeDotState,c) & + * crystallite_subdt(g,i,e) + if(present(DB)) & + broken = broken .or. .not. converged(matmul(source_RKdotState(1:sizeDotState,1:size(DB),s),DB) & + * crystallite_subdt(g,i,e), & + sourceState(p)%p(s)%state(1:sizeDotState,c), & + sourceState(p)%p(s)%atol(1:sizeDotState)) + enddo + if(broken) return + + broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_Fe(1:3,1:3,g,i,e), & + crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) + if(broken) return + + broken = integrateStress(g,i,e) + crystallite_converged(g,i,e) = .not. broken - if(nonlocalBroken) call nonlocalConvergenceCheck end subroutine integrateStateRK @@ -1624,7 +1546,19 @@ end subroutine integrateStateRK subroutine nonlocalConvergenceCheck integer :: e,i,p + logical :: nonlocal_broken + nonlocal_broken = .false. + !$OMP PARALLEL DO PRIVATE(p) + do e = FEsolving_execElem(1),FEsolving_execElem(2) + p = material_phaseAt(1,e) + do i = FEsolving_execIP(1),FEsolving_execIP(2) + if(plasticState(p)%nonlocal .and. .not. crystallite_converged(1,i,e)) nonlocal_broken = .true. + enddo + enddo + !$OMP END PARALLEL DO + + if(.not. nonlocal_broken) return !$OMP PARALLEL DO PRIVATE(p) do e = FEsolving_execElem(1),FEsolving_execElem(2) p = material_phaseAt(1,e) diff --git a/src/discretization.f90 b/src/discretization.f90 index 3bb905cf3..e6e53fcf4 100644 --- a/src/discretization.f90 +++ b/src/discretization.f90 @@ -15,7 +15,7 @@ module discretization discretization_nElem integer, public, protected, dimension(:), allocatable :: & - discretization_microstructureAt + discretization_materialAt real(pReal), public, protected, dimension(:,:), allocatable :: & discretization_IPcoords0, & @@ -37,12 +37,12 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief stores the relevant information in globally accesible variables !-------------------------------------------------------------------------------------------------- -subroutine discretization_init(microstructureAt,& +subroutine discretization_init(materialAt,& IPcoords0,NodeCoords0,& sharedNodesBegin) integer, dimension(:), intent(in) :: & - microstructureAt + materialAt real(pReal), dimension(:,:), intent(in) :: & IPcoords0, & NodeCoords0 @@ -51,10 +51,10 @@ subroutine discretization_init(microstructureAt,& print'(/,a)', ' <<<+- discretization init -+>>>'; flush(6) - discretization_nElem = size(microstructureAt,1) + discretization_nElem = size(materialAt,1) discretization_nIP = size(IPcoords0,2)/discretization_nElem - discretization_microstructureAt = microstructureAt + discretization_materialAt = materialAt discretization_IPcoords0 = IPcoords0 discretization_IPcoords = IPcoords0 diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index 8158996d6..a1bb23375 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -181,7 +181,7 @@ program DAMASK_grid if ((N_def /= N_n) .or. (N_n /= N_t) .or. N_n < 1) & ! sanity check call IO_error(error_ID=837,el=currentLoadCase,ext_msg = trim(interface_loadFile)) ! error message for incomplete loadcase - newLoadCase%stress%myType='stress' + newLoadCase%stress%myType='p' field = 1 newLoadCase%ID(field) = FIELD_MECH_ID ! mechanical active by default thermalActive: if (any(thermal_type == THERMAL_conduction_ID)) then @@ -210,18 +210,16 @@ program DAMASK_grid temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not a * if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable enddo - newLoadCase%deformation%maskLogical = transpose(reshape(temp_maskVector,[ 3,3])) ! logical mask in 3x3 notation - newLoadCase%deformation%maskFloat = merge(ones,zeros,newLoadCase%deformation%maskLogical) ! float (1.0/0.0) mask in 3x3 notation - newLoadCase%deformation%values = math_9to33(temp_valueVector) ! values in 3x3 notation + newLoadCase%deformation%mask = transpose(reshape(temp_maskVector,[ 3,3])) ! mask in 3x3 notation + newLoadCase%deformation%values = math_9to33(temp_valueVector) ! values in 3x3 notation case('p','stress', 's') temp_valueVector = 0.0_pReal do j = 1, 9 temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not an asterisk if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable enddo - newLoadCase%stress%maskLogical = transpose(reshape(temp_maskVector,[ 3,3])) - newLoadCase%stress%maskFloat = merge(ones,zeros,newLoadCase%stress%maskLogical) - newLoadCase%stress%values = math_9to33(temp_valueVector) + newLoadCase%stress%mask = transpose(reshape(temp_maskVector,[ 3,3])) + newLoadCase%stress%values = math_9to33(temp_valueVector) case('t','time','delta') ! increment time newLoadCase%time = IO_floatValue(line,chunkPos,i+1) case('n','incs','increments') ! number of increments @@ -268,8 +266,8 @@ program DAMASK_grid print*, ' drop guessing along trajectory' if (newLoadCase%deformation%myType == 'l') then do j = 1, 3 - if (any(newLoadCase%deformation%maskLogical(j,1:3) .eqv. .true.) .and. & - any(newLoadCase%deformation%maskLogical(j,1:3) .eqv. .false.)) errorID = 832 ! each row should be either fully or not at all defined + if (any(newLoadCase%deformation%mask(j,1:3) .eqv. .true.) .and. & + any(newLoadCase%deformation%mask(j,1:3) .eqv. .false.)) errorID = 832 ! each row should be either fully or not at all defined enddo print*, ' velocity gradient:' else if (newLoadCase%deformation%myType == 'f') then @@ -278,20 +276,19 @@ program DAMASK_grid print*, ' deformation gradient rate:' endif do i = 1, 3; do j = 1, 3 - if(newLoadCase%deformation%maskLogical(i,j)) then + if(newLoadCase%deformation%mask(i,j)) then write(IO_STDOUT,'(2x,f12.7)',advance='no') newLoadCase%deformation%values(i,j) else write(IO_STDOUT,'(2x,12a)',advance='no') ' * ' endif enddo; write(IO_STDOUT,'(/)',advance='no') enddo - if (any(newLoadCase%stress%maskLogical .eqv. & - newLoadCase%deformation%maskLogical)) errorID = 831 ! exclusive or masking only - if (any(newLoadCase%stress%maskLogical .and. transpose(newLoadCase%stress%maskLogical) & - .and. (math_I3<1))) errorID = 838 ! no rotation is allowed by stress BC + if (any(newLoadCase%stress%mask .eqv. newLoadCase%deformation%mask)) errorID = 831 ! exclusive or masking only + if (any(newLoadCase%stress%mask .and. transpose(newLoadCase%stress%mask) .and. (math_I3<1))) & + errorID = 838 ! no rotation is allowed by stress BC print*, ' stress / GPa:' do i = 1, 3; do j = 1, 3 - if(newLoadCase%stress%maskLogical(i,j)) then + if(newLoadCase%stress%mask(i,j)) then write(IO_STDOUT,'(2x,f12.7)',advance='no') newLoadCase%stress%values(i,j)*1e-9_pReal else write(IO_STDOUT,'(2x,12a)',advance='no') ' * ' @@ -368,11 +365,7 @@ program DAMASK_grid timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal) else if (currentLoadCase == 1) then ! 1st load case of logarithmic scale - if (inc == 1) then ! 1st inc of 1st load case of logarithmic scale - timeinc = loadCases(1)%time*(2.0_pReal**real( 1-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd - else ! not-1st inc of 1st load case of logarithmic scale - timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1-loadCases(1)%incs ,pReal)) - endif + timeinc = loadCases(1)%time*(2.0_pReal**real(max(inc-1,1)-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd else ! not-1st load case of logarithmic scale timeinc = time0 * & ( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc ,pReal)/& @@ -432,17 +425,11 @@ program DAMASK_grid do field = 1, nActiveFields select case(loadCases(currentLoadCase)%ID(field)) case(FIELD_MECH_ID) - solres(field) = mech_solution (& - incInfo,timeinc,timeIncOld, & - stress_BC = loadCases(currentLoadCase)%stress, & - rotation_BC = loadCases(currentLoadCase)%rot) - + solres(field) = mech_solution(incInfo) case(FIELD_THERMAL_ID) - solres(field) = grid_thermal_spectral_solution(timeinc,timeIncOld) - + solres(field) = grid_thermal_spectral_solution(timeinc) case(FIELD_DAMAGE_ID) - solres(field) = grid_damage_spectral_solution(timeinc,timeIncOld) - + solres(field) = grid_damage_spectral_solution(timeinc) end select if (.not. solres(field)%converged) exit ! no solution found diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 27dc5c1bd..7e529fc68 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -156,11 +156,10 @@ end subroutine grid_damage_spectral_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the spectral damage scheme with internal iterations !-------------------------------------------------------------------------------------------------- -function grid_damage_spectral_solution(timeinc,timeinc_old) result(solution) +function grid_damage_spectral_solution(timeinc) result(solution) real(pReal), intent(in) :: & - timeinc, & !< increment in time for current solution - timeinc_old !< increment in time of last increment + timeinc !< increment in time for current solution integer :: i, j, k, cell type(tSolutionState) :: solution PetscInt :: devNull @@ -174,7 +173,6 @@ function grid_damage_spectral_solution(timeinc,timeinc_old) result(solution) !-------------------------------------------------------------------------------------------------- ! set module wide availabe data params%timeinc = timeinc - params%timeincOld = timeinc_old call SNESSolve(damage_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) call SNESGetConvergedReason(damage_snes,reason,ierr); CHKERRQ(ierr) diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index e1ceb42b0..b680dafa9 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -25,63 +25,56 @@ module grid_mech_FEM implicit none private -!-------------------------------------------------------------------------------------------------- -! derived types - type(tSolutionParams), private :: params + type(tSolutionParams) :: params - type, private :: tNumerics + type :: 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 end type tNumerics - type(tNumerics), private :: num - logical, private:: & - debugRotation + type(tNumerics) :: num ! numerics parameters. Better name? + + logical :: debugRotation !-------------------------------------------------------------------------------------------------- ! PETSc data - DM, private :: mech_grid - SNES, private :: mech_snes - Vec, private :: solution_current, solution_lastInc, solution_rate + DM :: mech_grid + SNES :: mech_snes + Vec :: solution_current, solution_lastInc, solution_rate !-------------------------------------------------------------------------------------------------- ! common pointwise data - real(pReal), private, dimension(:,:,:,:,:), allocatable :: F, P_current, F_lastInc - real(pReal), private :: detJ - real(pReal), private, dimension(3) :: delta - real(pReal), private, dimension(3,8) :: BMat - real(pReal), private, dimension(8,8) :: HGMat - PetscInt, private :: xstart,ystart,zstart,xend,yend,zend + real(pReal), dimension(:,:,:,:,:), allocatable :: F, P_current, F_lastInc + real(pReal) :: detJ + real(pReal), dimension(3) :: delta + real(pReal), dimension(3,8) :: BMat + real(pReal), dimension(8,8) :: HGMat + PetscInt :: xstart,ystart,zstart,xend,yend,zend !-------------------------------------------------------------------------------------------------- ! stress, stiffness and compliance average etc. - real(pReal), private, dimension(3,3) :: & + real(pReal), dimension(3,3) :: & F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient F_aim = math_I3, & !< current prescribed deformation gradient - F_aim_lastIter = math_I3, & F_aim_lastInc = math_I3, & !< previous average deformation gradient - P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress - - character(len=:), allocatable, private :: incInfo !< time and increment information - - real(pReal), private, dimension(3,3,3,3) :: & + P_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress + P_aim = 0.0_pReal + character(len=:), allocatable :: incInfo !< time and increment information + real(pReal), dimension(3,3,3,3) :: & C_volAvg = 0.0_pReal, & !< current volume average stiffness C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness S = 0.0_pReal !< current compliance (filled up with zeros) - real(pReal), private :: & + real(pReal) :: & err_BC !< deviation from stress BC - integer, private :: & + integer :: & totalIter = 0 !< total iteration in current increment public :: & @@ -99,7 +92,6 @@ contains subroutine grid_mech_FEM_init real(pReal) :: HGCoeff = 0.0e-2_pReal - PetscInt, dimension(0:worldsize-1) :: localK real(pReal), dimension(3,3) :: & temp33_Real = 0.0_pReal real(pReal), dimension(4,8) :: & @@ -111,33 +103,34 @@ subroutine grid_mech_FEM_init -1.0_pReal, 1.0_pReal,-1.0_pReal,-1.0_pReal, & 1.0_pReal,-1.0_pReal,-1.0_pReal,-1.0_pReal, & 1.0_pReal, 1.0_pReal, 1.0_pReal, 1.0_pReal], [4,8]) + real(pReal), dimension(3,3,3,3) :: devNull PetscErrorCode :: ierr + PetscScalar, pointer, dimension(:,:,:,:) :: & + u_current,u_lastInc + PetscInt, dimension(0:worldsize-1) :: localK integer(HID_T) :: fileHandle, groupHandle character(len=pStringLen) :: & fileName class(tNode), pointer :: & num_grid, & debug_grid - real(pReal), dimension(3,3,3,3) :: devNull - PetscScalar, pointer, dimension(:,:,:,:) :: & - u_current,u_lastInc print'(/,a)', ' <<<+- grid_mech_FEM init -+>>>'; flush(IO_STDOUT) -!----------------------------------------------------------------------------------------------- +!------------------------------------------------------------------------------------------------- ! debugging options debug_grid => config_debug%get('grid', defaultVal=emptyList) debugRotation = debug_grid%contains('rotation') - + !------------------------------------------------------------------------------------------------- ! read numerical parameters and do sanity checks num_grid => config_numerics%get('grid',defaultVal=emptyDict) - 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) + 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=1.0e-3_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') @@ -242,6 +235,7 @@ subroutine grid_mech_FEM_init F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) endif restartRead + materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent call utilities_updateCoords(F) call utilities_constitutiveResponse(P_current,temp33_Real,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2 @@ -268,19 +262,12 @@ end subroutine grid_mech_FEM_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the FEM scheme with internal iterations !-------------------------------------------------------------------------------------------------- -function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution) +function grid_mech_FEM_solution(incInfoIn) result(solution) !-------------------------------------------------------------------------------------------------- ! input data for solution character(len=*), intent(in) :: & incInfoIn - real(pReal), intent(in) :: & - timeinc, & !< time increment of current solution - timeinc_old !< time increment of last successful increment - type(tBoundaryCondition), intent(in) :: & - stress_BC - type(rotation), intent(in) :: & - rotation_BC type(tSolutionState) :: & solution !-------------------------------------------------------------------------------------------------- @@ -292,14 +279,7 @@ function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) -!-------------------------------------------------------------------------------------------------- -! set module wide available data - params%stress_mask = stress_BC%maskFloat - params%stress_BC = stress_BC%values - params%rotation_BC = rotation_BC - params%timeinc = timeinc - params%timeincOld = timeinc_old + S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg) !-------------------------------------------------------------------------------------------------- ! solve BVP @@ -341,6 +321,14 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime, PetscScalar, pointer, dimension(:,:,:,:) :: & u_current,u_lastInc +!-------------------------------------------------------------------------------------------------- +! set module wide available data + params%stress_mask = stress_BC%mask + params%rotation_BC = rotation_BC + params%timeinc = timeinc + params%timeincOld = timeinc_old + + call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) @@ -349,20 +337,20 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime, else C_volAvgLastInc = C_volAvg - F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess) + F_aimDot = merge(merge((F_aim-F_aim_lastInc)/timeinc_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess) F_aim_lastInc = F_aim !----------------------------------------------------------------------------------------------- ! calculate rate for aim if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc) + F_aimDot = F_aimDot & + + merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask) elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * deformation_BC%values + F_aimDot = F_aimDot & + + merge(deformation_BC%values,.0_pReal,deformation_BC%mask) elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime + F_aimDot = F_aimDot & + + merge((deformation_BC%values - F_aim_lastInc)/loadCaseTime,.0_pReal,deformation_BC%mask) endif if (guess) then @@ -382,6 +370,12 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime, !-------------------------------------------------------------------------------------------------- ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * timeinc + if (stress_BC%myType=='p') then + P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask) + elseif (stress_BC%myType=='pdot') then !UNTESTED + P_aim = P_aim + merge(stress_BC%values*timeinc,.0_pReal,stress_BC%mask) + endif + call VecAXPY(solution_current,timeinc,solution_rate,ierr); CHKERRQ(ierr) call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr);CHKERRQ(ierr) @@ -497,8 +491,6 @@ subroutine formResidual(da_local,x_local, & PetscScalar, pointer,dimension(:,:,:,:) :: x_scal, f_scal PetscScalar, dimension(8,3) :: x_elem, f_elem PetscInt :: i, ii, j, jj, k, kk, ctr, ele - real(pReal), dimension(3,3) :: & - deltaF_aim PetscInt :: & PETScIter, & nfuncs @@ -547,10 +539,8 @@ subroutine formResidual(da_local,x_local, & !-------------------------------------------------------------------------------------------------- ! stress BC handling - F_aim_lastIter = F_aim - deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC) - F_aim = F_aim - deltaF_aim - err_BC = maxval(abs(params%stress_mask * (P_av - params%stress_BC))) ! mask = 0.0 when no stress bc + F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc + err_BC = maxval(abs(merge(P_av - P_aim,.0_pReal,params%stress_mask))) !-------------------------------------------------------------------------------------------------- ! constructing residual diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index e7fb9bd19..e7e544886 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -24,11 +24,9 @@ module grid_mech_spectral_basic implicit none private -!-------------------------------------------------------------------------------------------------- -! derived types type(tSolutionParams) :: params - type, private :: tNumerics + type :: tNumerics logical :: update_gamma !< update gamma operator with current stiffness integer :: & itmin, & !< minimum number of iterations @@ -42,7 +40,7 @@ module grid_mech_spectral_basic type(tNumerics) :: num ! numerics parameters. Better name? - logical, private :: debugRotation + logical :: debugRotation !-------------------------------------------------------------------------------------------------- ! PETSc data @@ -62,10 +60,10 @@ module grid_mech_spectral_basic F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient F_aim = math_I3, & !< current prescribed deformation gradient F_aim_lastInc = math_I3, & !< previous average deformation gradient - P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress - + P_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress + P_aim = 0.0_pReal character(len=:), allocatable :: incInfo !< time and increment information - real(pReal), private, dimension(3,3,3,3) :: & + real(pReal), dimension(3,3,3,3) :: & C_volAvg = 0.0_pReal, & !< current volume average stiffness C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 stiffness @@ -96,18 +94,17 @@ subroutine grid_mech_spectral_basic_init real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P real(pReal), dimension(3,3) :: & temp33_Real = 0.0_pReal - class (tNode), pointer :: & - num_grid, & - debug_grid - PetscErrorCode :: ierr PetscScalar, pointer, dimension(:,:,:,:) :: & F ! pointer to solution data - PetscInt, dimension(worldsize) :: localK + PetscInt, dimension(0:worldsize-1) :: localK integer(HID_T) :: fileHandle, groupHandle integer :: fileUnit character(len=pStringLen) :: & fileName + class (tNode), pointer :: & + num_grid, & + debug_grid print'(/,a)', ' <<<+- grid_mech_spectral_basic init -+>>>'; flush(IO_STDOUT) @@ -126,13 +123,13 @@ subroutine grid_mech_spectral_basic_init ! read numerical parameters and do sanity checks num_grid => config_numerics%get('grid',defaultVal=emptyDict) - 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) + 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=1.0e-3_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') @@ -158,7 +155,7 @@ subroutine grid_mech_spectral_basic_init call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) localK = 0 - localK(worldrank+1) = grid3 + localK(worldrank) = grid3 call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) call DMDACreate3d(PETSC_COMM_WORLD, & DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary @@ -202,8 +199,8 @@ subroutine grid_mech_spectral_basic_init endif restartRead materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent - call Utilities_updateCoords(reshape(F,shape(F_lastInc))) - call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 + call utilities_updateCoords(reshape(F,shape(F_lastInc))) + call utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 reshape(F,shape(F_lastInc)), & ! target F 0.0_pReal) ! time increment call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! deassociate pointer @@ -231,19 +228,12 @@ end subroutine grid_mech_spectral_basic_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the basic scheme with internal iterations !-------------------------------------------------------------------------------------------------- -function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution) +function grid_mech_spectral_basic_solution(incInfoIn) result(solution) !-------------------------------------------------------------------------------------------------- ! input data for solution character(len=*), intent(in) :: & incInfoIn - real(pReal), intent(in) :: & - timeinc, & !< time increment of current solution - timeinc_old !< time increment of last successful increment - type(tBoundaryCondition), intent(in) :: & - stress_BC - type(rotation), intent(in) :: & - rotation_BC type(tSolutionState) :: & solution !-------------------------------------------------------------------------------------------------- @@ -255,17 +245,9 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_ !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) + S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg) if(num%update_gamma) call utilities_updateGamma(C_minMaxAvg) -!-------------------------------------------------------------------------------------------------- -! set module wide available data - params%stress_mask = stress_BC%maskFloat - params%stress_BC = stress_BC%values - params%rotation_BC = rotation_BC - params%timeinc = timeinc - params%timeincOld = timeinc_old - !-------------------------------------------------------------------------------------------------- ! solve BVP call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) @@ -303,7 +285,14 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo type(rotation), intent(in) :: & rotation_BC PetscErrorCode :: ierr - PetscScalar, dimension(:,:,:,:), pointer :: F + PetscScalar, pointer, dimension(:,:,:,:) :: F + +!-------------------------------------------------------------------------------------------------- +! set module wide available data + params%stress_mask = stress_BC%mask + params%rotation_BC = rotation_BC + params%timeinc = timeinc + params%timeincOld = timeinc_old call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) @@ -314,20 +303,20 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo C_volAvgLastInc = C_volAvg C_minMaxAvgLastInc = C_minMaxAvg - F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess) + F_aimDot = merge(merge((F_aim-F_aim_lastInc)/timeinc_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess) F_aim_lastInc = F_aim !----------------------------------------------------------------------------------------------- ! calculate rate for aim if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc) + F_aimDot = F_aimDot & + + merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask) elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * deformation_BC%values + F_aimDot = F_aimDot & + + merge(deformation_BC%values,.0_pReal,deformation_BC%mask) elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime + F_aimDot = F_aimDot & + + merge((deformation_BC%values - F_aim_lastInc)/loadCaseTime,.0_pReal,deformation_BC%mask) endif Fdot = utilities_calculateRate(guess, & @@ -341,7 +330,13 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo !-------------------------------------------------------------------------------------------------- ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * timeinc - F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average + if (stress_BC%myType=='p') then + P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask) + elseif (stress_BC%myType=='pdot') then !UNTESTED + P_aim = P_aim + merge(stress_BC%values*timeinc,.0_pReal,stress_BC%mask) + endif + + F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average rotation_BC%rotate(F_aim,active=.true.)),[9,grid(1),grid(2),grid3]) call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) @@ -375,7 +370,7 @@ subroutine grid_mech_spectral_basic_restartWrite call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) - print'(a)', ' writing solver data required for restart to file'; flush(IO_STDOUT) + print*, 'writing solver data required for restart to file'; flush(IO_STDOUT) write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5' fileHandle = HDF5_openFile(fileName,'w') @@ -469,6 +464,7 @@ subroutine formResidual(in, F, & 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 @@ -491,16 +487,16 @@ subroutine formResidual(in, F, & !-------------------------------------------------------------------------------------------------- ! stress BC handling - deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC) + deltaF_aim = math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc F_aim = F_aim - deltaF_aim - err_BC = maxval(abs(params%stress_mask * (P_av - params%stress_BC))) ! mask = 0.0 when no stress bc + err_BC = maxval(abs(merge(P_av - P_aim,.0_pReal,params%stress_mask))) !-------------------------------------------------------------------------------------------------- ! updated deformation gradient using fix point algorithm of basic scheme tensorField_real = 0.0_pReal tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residuum ! store fPK field for subsequent FFT forward transform call utilities_FFTtensorForward ! FFT forward of global "tensorField_real" - err_div = Utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use + err_div = utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use call utilities_fourierGammaConvolution(params%rotation_BC%rotate(deltaF_aim,active=.true.)) ! convolution of Gamma and tensorField_fourier call utilities_FFTtensorBackward ! FFT backward of global tensorField_fourier diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 5cc7d1602..f59b68d7a 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -15,7 +15,6 @@ module grid_mech_spectral_polarisation use DAMASK_interface use HDF5_utilities use math - use rotations use spectral_utilities use FEsolving use config @@ -25,8 +24,6 @@ module grid_mech_spectral_polarisation implicit none private -!-------------------------------------------------------------------------------------------------- -! derived types type(tSolutionParams) :: params type :: tNumerics @@ -48,7 +45,7 @@ module grid_mech_spectral_polarisation type(tNumerics) :: num ! numerics parameters. Better name? - logical, private :: debugRotation + logical :: debugRotation !-------------------------------------------------------------------------------------------------- ! PETSc data @@ -71,8 +68,8 @@ module grid_mech_spectral_polarisation F_aim = math_I3, & !< current prescribed deformation gradient F_aim_lastInc = math_I3, & !< previous average deformation gradient F_av = 0.0_pReal, & !< average incompatible def grad field - P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress - + P_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress + P_aim = 0.0_pReal character(len=:), allocatable :: incInfo !< time and increment information real(pReal), dimension(3,3,3,3) :: & C_volAvg = 0.0_pReal, & !< current volume average stiffness @@ -108,10 +105,6 @@ subroutine grid_mech_spectral_polarisation_init real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P real(pReal), dimension(3,3) :: & temp33_Real = 0.0_pReal - class (tNode), pointer :: & - num_grid, & - debug_grid - PetscErrorCode :: ierr PetscScalar, pointer, dimension(:,:,:,:) :: & FandF_tau, & ! overall pointer to solution data @@ -122,13 +115,16 @@ subroutine grid_mech_spectral_polarisation_init integer :: fileUnit character(len=pStringLen) :: & fileName + class (tNode), pointer :: & + num_grid, & + debug_grid print'(/,a)', ' <<<+- grid_mech_spectral_polarisation init -+>>>'; flush(IO_STDOUT) print*, 'Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006' -!------------------------------------------------------------------------------------------------ +!------------------------------------------------------------------------------------------------- ! debugging options debug_grid => config_debug%get('grid',defaultVal=emptyList) debugRotation = debug_grid%contains('rotation') @@ -136,17 +132,18 @@ subroutine grid_mech_spectral_polarisation_init !------------------------------------------------------------------------------------------------- ! read numerical parameters and do sanity checks num_grid => config_numerics%get('grid',defaultVal=emptyDict) - 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) + + 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=1.0e-3_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') @@ -228,8 +225,8 @@ subroutine grid_mech_spectral_polarisation_init endif restartRead materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent - call Utilities_updateCoords(reshape(F,shape(F_lastInc))) - call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 + call utilities_updateCoords(reshape(F,shape(F_lastInc))) + call utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 reshape(F,shape(F_lastInc)), & ! target F 0.0_pReal) ! time increment call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! deassociate pointer @@ -259,19 +256,12 @@ end subroutine grid_mech_spectral_polarisation_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the Polarisation scheme with internal iterations !-------------------------------------------------------------------------------------------------- -function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution) +function grid_mech_spectral_polarisation_solution(incInfoIn) result(solution) !-------------------------------------------------------------------------------------------------- ! input data for solution character(len=*), intent(in) :: & incInfoIn - real(pReal), intent(in) :: & - timeinc, & !< time increment of current solution - timeinc_old !< time increment of last successful increment - type(tBoundaryCondition), intent(in) :: & - stress_BC - type(rotation), intent(in) :: & - rotation_BC type(tSolutionState) :: & solution !-------------------------------------------------------------------------------------------------- @@ -283,21 +273,13 @@ function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old, !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) - if (num%update_gamma) then + S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg) + if(num%update_gamma) then call utilities_updateGamma(C_minMaxAvg) C_scale = C_minMaxAvg S_scale = math_invSym3333(C_minMaxAvg) endif -!-------------------------------------------------------------------------------------------------- -! set module wide available data - params%stress_mask = stress_BC%maskFloat - params%stress_BC = stress_BC%values - params%rotation_BC = rotation_BC - params%timeinc = timeinc - params%timeincOld = timeinc_old - !-------------------------------------------------------------------------------------------------- ! solve BVP call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) @@ -335,10 +317,17 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc type(rotation), intent(in) :: & rotation_BC PetscErrorCode :: ierr - PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau, F, F_tau + PetscScalar, pointer, dimension(:,:,:,:) :: FandF_tau, F, F_tau integer :: i, j, k real(pReal), dimension(3,3) :: F_lambda33 +!-------------------------------------------------------------------------------------------------- +! set module wide available data + params%stress_mask = stress_BC%mask + params%rotation_BC = rotation_BC + params%timeinc = timeinc + params%timeincOld = timeinc_old + call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) F => FandF_tau(0: 8,:,:,:) F_tau => FandF_tau(9:17,:,:,:) @@ -350,20 +339,20 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc C_volAvgLastInc = C_volAvg C_minMaxAvgLastInc = C_minMaxAvg - F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess) + F_aimDot = merge(merge((F_aim-F_aim_lastInc)/timeinc_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess) F_aim_lastInc = F_aim !----------------------------------------------------------------------------------------------- ! calculate rate for aim if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc) + F_aimDot = F_aimDot & + + merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask) elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * deformation_BC%values + F_aimDot = F_aimDot & + + merge(deformation_BC%values,.0_pReal,deformation_BC%mask) elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime + F_aimDot = F_aimDot & + + merge((deformation_BC%values - F_aim_lastInc)/loadCaseTime,.0_pReal,deformation_BC%mask) endif Fdot = utilities_calculateRate(guess, & @@ -381,6 +370,12 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc !-------------------------------------------------------------------------------------------------- ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * timeinc + if (stress_BC%myType=='p') then + P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask) + elseif (stress_BC%myType=='pdot') then !UNTESTED + P_aim = P_aim + merge(stress_BC%values*timeinc,.0_pReal,stress_BC%mask) + endif + F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average rotation_BC%rotate(F_aim,active=.true.)),& [9,grid(1),grid(2),grid3]) @@ -595,15 +590,15 @@ subroutine formResidual(in, FandF_tau, & !-------------------------------------------------------------------------------------------------- ! stress BC handling - F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc - err_BC = maxval(abs((1.0_pReal-params%stress_mask) * math_mul3333xx33(C_scale,F_aim & - -params%rotation_BC%rotate(F_av)) + & - params%stress_mask * (P_av-params%stress_BC))) ! mask = 0.0 for no bc + F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc + err_BC = maxval(abs(merge(P_av-P_aim, & + math_mul3333xx33(C_scale,F_aim-params%rotation_BC%rotate(F_av)),& + params%stress_mask))) ! calculate divergence tensorField_real = 0.0_pReal tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F !< stress field in disguise call utilities_FFTtensorForward - err_div = Utilities_divergenceRMS() !< root mean squared error in divergence of stress + err_div = utilities_divergenceRMS() !< root mean squared error in divergence of stress !-------------------------------------------------------------------------------------------------- ! constructing residual @@ -622,7 +617,7 @@ subroutine formResidual(in, FandF_tau, & tensorField_real = 0.0_pReal tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F call utilities_FFTtensorForward - err_curl = Utilities_curlRMS() + err_curl = utilities_curlRMS() end subroutine formResidual diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 49be5ad7e..b4f9acddb 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -148,11 +148,10 @@ end subroutine grid_thermal_spectral_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the spectral thermal scheme with internal iterations !-------------------------------------------------------------------------------------------------- -function grid_thermal_spectral_solution(timeinc,timeinc_old) result(solution) +function grid_thermal_spectral_solution(timeinc) result(solution) real(pReal), intent(in) :: & - timeinc, & !< increment in time for current solution - timeinc_old !< increment in time of last increment + timeinc !< increment in time for current solution integer :: i, j, k, cell type(tSolutionState) :: solution PetscInt :: devNull @@ -166,7 +165,6 @@ function grid_thermal_spectral_solution(timeinc,timeinc_old) result(solution) !-------------------------------------------------------------------------------------------------- ! set module wide availabe data params%timeinc = timeinc - params%timeincOld = timeinc_old call SNESSolve(thermal_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) call SNESGetConvergedReason(thermal_snes,reason,ierr); CHKERRQ(ierr) diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index dc6430c72..73aaa7789 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -47,15 +47,15 @@ module spectral_utilities complex(C_DOUBLE_COMPLEX),public, dimension(:,:,:,:), pointer :: vectorField_fourier !< vector field fourier representation for fftw real(C_DOUBLE), public, dimension(:,:,:), pointer :: scalarField_real !< scalar field real representation for fftw complex(C_DOUBLE_COMPLEX),public, dimension(:,:,:), pointer :: scalarField_fourier !< scalar field fourier representation for fftw - complex(pReal), private, dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat !< gamma operator (field) for spectral method - complex(pReal), private, dimension(:,:,:,:), allocatable :: xi1st !< wave vector field for first derivatives - complex(pReal), private, dimension(:,:,:,:), allocatable :: xi2nd !< wave vector field for second derivatives - real(pReal), private, dimension(3,3,3,3) :: C_ref !< mechanic reference stiffness + complex(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat !< gamma operator (field) for spectral method + complex(pReal), dimension(:,:,:,:), allocatable :: xi1st !< wave vector field for first derivatives + complex(pReal), dimension(:,:,:,:), allocatable :: xi2nd !< wave vector field for second derivatives + real(pReal), dimension(3,3,3,3) :: C_ref !< mechanic reference stiffness !-------------------------------------------------------------------------------------------------- ! plans for FFTW - type(C_PTR), private :: & + type(C_PTR) :: & planTensorForth, & !< FFTW MPI plan P(x) to P(k) planTensorBack, & !< FFTW MPI plan F(k) to F(x) planVectorForth, & !< FFTW MPI plan v(x) to v(k) @@ -65,7 +65,7 @@ module spectral_utilities !-------------------------------------------------------------------------------------------------- ! variables controlling debugging - logical, private :: & + logical :: & debugGeneral, & !< general debugging of spectral solver debugRotation, & !< also printing out results in lab frame debugPETSc !< use some in debug defined options for more verbose PETSc solution @@ -82,10 +82,9 @@ module spectral_utilities end type tSolutionState type, public :: tBoundaryCondition !< set of parameters defining a boundary condition - real(pReal), dimension(3,3) :: values = 0.0_pReal, & - maskFloat = 0.0_pReal - logical, dimension(3,3) :: maskLogical = .false. - character(len=pStringLen) :: myType = 'None' + real(pReal), dimension(3,3) :: values = 0.0_pReal + logical, dimension(3,3) :: mask = .false. + character(len=pStringLen) :: myType = 'None' end type tBoundaryCondition type, public :: tLoadCase @@ -101,21 +100,21 @@ module spectral_utilities integer(kind(FIELD_UNDEFINED_ID)), allocatable :: ID(:) end type tLoadCase - type, public :: tSolutionParams !< @todo use here the type definition for a full loadcase - real(pReal), dimension(3,3) :: stress_mask, stress_BC + type, public :: tSolutionParams + real(pReal), dimension(3,3) :: stress_BC + logical, dimension(3,3) :: stress_mask type(rotation) :: rotation_BC - real(pReal) :: timeinc - real(pReal) :: timeincOld + real(pReal) :: timeinc, timeincOld end type tSolutionParams - type, private :: tNumerics + type :: tNumerics integer :: & divergence_correction !< scale divergence/curl calculation: [0: no correction, 1: size scaled to 1, 2: size scaled to Npoints] logical :: & memory_efficient !< calculate gamma operator on the fly end type tNumerics - type(tNumerics), private :: num ! numerics parameters. Better name? + type(tNumerics) :: num ! numerics parameters. Better name? enum, bind(c); enumerator :: & DERIVATIVE_CONTINUOUS_ID, & diff --git a/src/material.f90 b/src/material.f90 index 6688a0ae5..742a29b9d 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -2,7 +2,7 @@ !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @brief Parses material config file, either solverJobName.materialConfig or material.config +!> @brief Defines phase and homogenization !-------------------------------------------------------------------------------------------------- module material use prec @@ -97,7 +97,7 @@ module material material_orientation0 !< initial orientation of each grain,IP,element integer, dimension(:), allocatable, private :: & - microstructure_Nconstituents !< number of constituents in each microstructure + material_Nconstituents !< number of constituents in each material @@ -180,8 +180,8 @@ subroutine material_init(restart) material_name_homogenization(myHomog) = trim(adjustl(sectionName))//material_homogenization%getKey(myHomog) enddo - call material_parseMicrostructure - print*, 'Microstructure parsed' + call material_parseMaterial + print*, 'Material parsed' call material_parseHomogenization print*, 'Homogenization parsed' @@ -317,16 +317,16 @@ end subroutine material_parseHomogenization !-------------------------------------------------------------------------------------------------- -!> @brief parses the microstructure part in the material configuration file +!> @brief parses the material part in the material configuration file !-------------------------------------------------------------------------------------------------- -subroutine material_parseMicrostructure +subroutine material_parseMaterial - class(tNode), pointer :: microstructures, & !> list of microstructures - microstructure, & !> microstructure definition + class(tNode), pointer :: materials, & !> list of materials + material, & !> material definition constituents, & !> list of constituents constituent, & !> constituent definition phases, & - homogenization + homogenizations integer, dimension(:), allocatable :: & counterPhase, & @@ -341,17 +341,17 @@ subroutine material_parseMicrostructure c, & maxNconstituents - microstructures => config_material%get('microstructure') - if(any(discretization_microstructureAt > microstructures%length)) & - call IO_error(155,ext_msg='More microstructures requested than found in material.yaml') + materials => config_material%get('material') + if(any(discretization_materialAt > materials%length)) & + call IO_error(155,ext_msg='More materials requested than found in material.yaml') - allocate(microstructure_Nconstituents(microstructures%length),source=0) - do m = 1, microstructures%length - microstructure => microstructures%get(m) - constituents => microstructure%get('constituents') - microstructure_Nconstituents(m) = constituents%length + allocate(material_Nconstituents(materials%length),source=0) + do m = 1, materials%length + material => materials%get(m) + constituents => material%get('constituents') + material_Nconstituents(m) = constituents%length enddo - maxNconstituents = maxval(microstructure_Nconstituents) + maxNconstituents = maxval(material_Nconstituents) allocate(material_homogenizationAt(discretization_nElem),source=0) allocate(material_homogenizationMemberAt(discretization_nIP,discretization_nElem),source=0) @@ -362,14 +362,14 @@ subroutine material_parseMicrostructure phases => config_material%get('phase') allocate(counterPhase(phases%length),source=0) - homogenization => config_material%get('homogenization') - allocate(counterHomogenization(homogenization%length),source=0) + homogenizations => config_material%get('homogenization') + allocate(counterHomogenization(homogenizations%length),source=0) do e = 1, discretization_nElem - microstructure => microstructures%get(discretization_microstructureAt(e)) - constituents => microstructure%get('constituents') + material => materials%get(discretization_materialAt(e)) + constituents => material%get('constituents') - material_homogenizationAt(e) = homogenization%getIndex(microstructure%get_asString('homogenization')) + material_homogenizationAt(e) = homogenizations%getIndex(material%get_asString('homogenization')) do i = 1, discretization_nIP counterHomogenization(material_homogenizationAt(e)) = counterHomogenization(material_homogenizationAt(e)) + 1 material_homogenizationMemberAt(i,e) = counterHomogenization(material_homogenizationAt(e)) @@ -385,7 +385,7 @@ subroutine material_parseMicrostructure counterPhase(material_phaseAt(c,e)) = counterPhase(material_phaseAt(c,e)) + 1 material_phaseMemberAt(c,i,e) = counterPhase(material_phaseAt(c,e)) - call material_orientation0(c,i,e)%fromQuaternion(constituent%get_asFloats('orientation',requiredSize=4)) + call material_orientation0(c,i,e)%fromQuaternion(constituent%get_asFloats('O',requiredSize=4)) enddo enddo @@ -393,7 +393,7 @@ subroutine material_parseMicrostructure enddo -end subroutine material_parseMicrostructure +end subroutine material_parseMaterial end module material diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index b57a4b332..bc96951a1 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -78,7 +78,7 @@ subroutine discretization_mesh_init(restart) IS :: faceSetIS PetscErrorCode :: ierr integer, dimension(:), allocatable :: & - microstructureAt + materialAt class(tNode), pointer :: & num_mesh integer :: integrationOrder !< order of quadrature rule required @@ -164,9 +164,9 @@ subroutine discretization_mesh_init(restart) call mesh_FEM_build_ipCoordinates(dimPlex,FEM_quadrature_points(dimPlex,integrationOrder)%p) call mesh_FEM_build_ipVolumes(dimPlex) - allocate(microstructureAt(mesh_NcpElems)) + allocate(materialAt(mesh_NcpElems)) do j = 1, mesh_NcpElems - call DMGetLabelValue(geomMesh,'material',j-1,microstructureAt(j),ierr) + call DMGetLabelValue(geomMesh,'material',j-1,materialAt(j),ierr) CHKERRQ(ierr) end do @@ -178,7 +178,7 @@ subroutine discretization_mesh_init(restart) allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal) - call discretization_init(microstructureAt,& + call discretization_init(materialAt,& reshape(mesh_ipCoordinates,[3,mesh_maxNips*mesh_NcpElems]), & mesh_node0) diff --git a/src/rotations.f90 b/src/rotations.f90 index e19513866..0b72c1dd5 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -1452,7 +1452,7 @@ subroutine selfTest contains - function quaternion_equal(qu1,qu2) result(ok) + pure recursive function quaternion_equal(qu1,qu2) result(ok) real(pReal), intent(in), dimension(4) :: qu1,qu2 logical :: ok