From 540428aab7868bec93ed90b5016cfb7aef3823c0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 29 Mar 2020 19:06:17 +0200 Subject: [PATCH 01/35] works also for ifort --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 0f488ec37..f6870137f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -148,7 +148,7 @@ else () set(OPENMP "${OPENMP}") endif () -# syntax check only (mainly for pre-receive hook, works only with gfortran) +# syntax check only (mainly for pre-receive hook) if (CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY") set (BUILDCMD_POST "${BUILDCMD_POST} -fsyntax-only") endif () From f1b4d81fb41f50f7911e3e1a29b41aa91d61a792 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 29 Mar 2020 19:11:36 +0200 Subject: [PATCH 02/35] simplified --- .../pre/geom_fromVoronoiTessellation.py | 34 +++++++++++-------- 1 file changed, 19 insertions(+), 15 deletions(-) diff --git a/processing/pre/geom_fromVoronoiTessellation.py b/processing/pre/geom_fromVoronoiTessellation.py index 3ecf7c099..edcbe83dc 100755 --- a/processing/pre/geom_fromVoronoiTessellation.py +++ b/processing/pre/geom_fromVoronoiTessellation.py @@ -21,36 +21,42 @@ def findClosestSeed(seeds, weights, point): return np.argmin(np.sum((np.broadcast_to(point,(len(seeds),3))-seeds)**2,axis=1) - weights) -def Laguerre_tessellation(grid, seeds, grains, size, periodic, weights, cpus): +def Laguerre_tessellation(grid, size, seeds, weights, origin = np.zeros(3), periodic = True, cpus = 2): if periodic: weights_p = np.tile(weights,27).flatten(order='F') # Laguerre weights (1,2,3,1,2,3,...,1,2,3) seeds_p = np.vstack((seeds -np.array([size[0],0.,0.]),seeds, seeds +np.array([size[0],0.,0.]))) seeds_p = np.vstack((seeds_p-np.array([0.,size[1],0.]),seeds_p,seeds_p+np.array([0.,size[1],0.]))) seeds_p = np.vstack((seeds_p-np.array([0.,0.,size[2]]),seeds_p,seeds_p+np.array([0.,0.,size[2]]))) + coords = damask.grid_filters.cell_coord0(grid*3,size*3,-origin-size).reshape(-1,3,order='F') else: weights_p = weights.flatten() seeds_p = seeds + coords = damask.grid_filters.cell_coord0(grid,size,-origin).reshape(-1,3,order='F') if cpus > 1: - default_args = partial(findClosestSeed,seeds_p,weights_p) - pool = multiprocessing.Pool(processes = cpus) # initialize workers - result = pool.map_async(default_args, [point for point in grid]) # evaluate function in parallel + pool = multiprocessing.Pool(processes = cpus) + result = pool.map_async(partial(findClosestSeed,seeds_p,weights_p), [coord for coord in coords]) pool.close() pool.join() - closestSeeds = np.array(result.get()).flatten() + closest_seed = np.array(result.get()) else: - closestSeeds= np.array([findClosestSeed(seeds_p,weights_p,point) for point in grid]) + closest_seed= np.array([findClosestSeed(seeds_p,weights_p,coord) for coord in coords]) - return grains[closestSeeds%seeds.shape[0]] + if periodic: + closest_seed = closest_seed.reshape(grid[2]*3,grid[1]*3,grid[0]*3) + return closest_seed[grid[2]:grid[2]*2,grid[1]:grid[1]*2,grid[0]:grid[0]*2]%seeds.shape[0] + else: + return closest_seed -def Voronoi_tessellation(grid, seeds, grains, size, periodic = True): +def Voronoi_tessellation(grid, size, seeds, origin = np.zeros(3), periodic = True): + coords = damask.grid_filters.cell_coord0(grid,size,-origin).reshape(-1,3,order='F') KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds) - devNull,closestSeeds = KDTree.query(grid) + devNull,closest_seed = KDTree.query(coords) - return grains[closestSeeds] + return closest_seed # -------------------------------------------------------------------- @@ -191,13 +197,11 @@ for name in filenames: if options.eulers in table.labels: eulers = table.get(options.eulers) - coords = damask.grid_filters.cell_coord0(grid,size,-origin).reshape(-1,3,order='F') - if options.laguerre: - indices = Laguerre_tessellation(coords,seeds,grains,size,options.periodic, - table.get(options.weight),options.cpus) + indices = grains[Laguerre_tessellation(grid,size,seeds,table.get(options.weight),origin, + options.periodic,options.cpus)] else: - indices = Voronoi_tessellation (coords,seeds,grains,size,options.periodic) + indices = grains[Voronoi_tessellation (grid,size,seeds,origin,options.periodic)] config_header = [] if options.config: From e61c1a027b2d326c595baf5d0de5d41a0936e00b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 29 Mar 2020 19:12:23 +0200 Subject: [PATCH 03/35] avoid detour via shell --- python/damask/_geom.py | 97 +++++++++++++++++++++++++++++++++++---- python/tests/test_Geom.py | 21 +++++++++ 2 files changed, 109 insertions(+), 9 deletions(-) diff --git a/python/damask/_geom.py b/python/damask/_geom.py index 6bc92e300..797b9dc81 100644 --- a/python/damask/_geom.py +++ b/python/damask/_geom.py @@ -1,11 +1,15 @@ import sys from io import StringIO +import multiprocessing +from functools import partial import numpy as np -from scipy import ndimage +from scipy import ndimage,spatial from . import VTK from . import util +from . import Environment +from . import grid_filters class Geom: @@ -50,7 +54,7 @@ class Geom: def update(self,microstructure=None,size=None,origin=None,rescale=False): """ - Updates microstructure and size. + Update microstructure and size. Parameters ---------- @@ -113,7 +117,7 @@ class Geom: def set_comments(self,comments): """ - Replaces all existing comments. + Replace all existing comments. Parameters ---------- @@ -127,7 +131,7 @@ class Geom: def add_comments(self,comments): """ - Appends comments to existing comments. + Append comments to existing comments. Parameters ---------- @@ -140,7 +144,7 @@ class Geom: def set_microstructure(self,microstructure): """ - Replaces the existing microstructure representation. + Replace the existing microstructure representation. Parameters ---------- @@ -159,7 +163,7 @@ class Geom: def set_size(self,size): """ - Replaces the existing size information. + Replace the existing size information. Parameters ---------- @@ -179,7 +183,7 @@ class Geom: def set_origin(self,origin): """ - Replaces the existing origin information. + Replace the existing origin information. Parameters ---------- @@ -196,7 +200,7 @@ class Geom: def set_homogenization(self,homogenization): """ - Replaces the existing homogenization index. + Replace the existing homogenization index. Parameters ---------- @@ -264,7 +268,7 @@ class Geom: @staticmethod def from_file(fname): """ - Reads a geom file. + Read a geom file. Parameters ---------- @@ -325,6 +329,81 @@ class Geom: return Geom(microstructure.reshape(grid),size,origin,homogenization,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) + @staticmethod + def from_Laguerre_tessellation(grid,size,seeds,weights,periodic=True): + """ + Generate geometry from Laguerre tessellation. + + Parameters + ---------- + grid : 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 microstructure in meter. + seeds : numpy.ndarray of shape (:,3) + position of the seed points in meter. All points need to lay within the box. + weights : numpy.ndarray of shape (seeds.shape[0]) + weights of the seeds. Setting all weights to 1.0 gives a standard Voronoi tessellation. + periodic : Boolean, optional + perform a periodic tessellation. Defaults to True. + + """ + if periodic: + weights_p = np.tile(weights,27).flatten(order='F') # Laguerre weights (1,2,3,1,2,3,...,1,2,3) + seeds_p = np.vstack((seeds -np.array([size[0],0.,0.]),seeds, seeds +np.array([size[0],0.,0.]))) + seeds_p = np.vstack((seeds_p-np.array([0.,size[1],0.]),seeds_p,seeds_p+np.array([0.,size[1],0.]))) + seeds_p = np.vstack((seeds_p-np.array([0.,0.,size[2]]),seeds_p,seeds_p+np.array([0.,0.,size[2]]))) + coords = grid_filters.cell_coord0(grid*3,size*3,-size).reshape(-1,3,order='F') + + else: + weights_p = weights.flatten() + seeds_p = seeds + coords = grid_filters.cell_coord0(grid,size).reshape(-1,3,order='F') + + pool = multiprocessing.Pool(processes = int(Environment().options['DAMASK_NUM_THREADS'])) + result = pool.map_async(partial(Geom._find_closest_seed,seeds_p,weights_p), [coord for coord in coords]) + pool.close() + pool.join() + microstructure = np.array(result.get()) + + if periodic: + microstructure = microstructure.reshape(grid[0]*3,grid[1]*3,grid[2]*3) + microstructure = microstructure[grid[0]:grid[0]*2,grid[1]:grid[1]*2,grid[2]:grid[2]*2]%seeds.shape[0] + else: + microstructure = microstructure.reshape(grid) + + #comments = 'geom.py:from_Laguerre_tessellation v{}'.format(version) + return Geom(microstructure+1,size,homogenization=1) + + + @staticmethod + def from_Voronoi_tessellation(grid,size,seeds,periodic=True): + """ + Generate geometry from Voronoi tessellation. + + Parameters + ---------- + grid : 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 microstructure in meter. + seeds : numpy.ndarray of shape (:,3) + position of the seed points in meter. All points need to lay within the box. + periodic : Boolean, optional + perform a periodic tessellation. Defaults to True. + + """ + coords = grid_filters.cell_coord0(grid,size).reshape(-1,3,order='F') + KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds) + devNull,microstructure = KDTree.query(coords) + + #comments = 'geom.py:from_Voronoi_tessellation v{}'.format(version) + return Geom(microstructure.reshape(grid)+1,size,homogenization=1) + + def to_file(self,fname,pack=None): """ Writes a geom file. diff --git a/python/tests/test_Geom.py b/python/tests/test_Geom.py index a6fd57cb6..1b1529e80 100644 --- a/python/tests/test_Geom.py +++ b/python/tests/test_Geom.py @@ -97,3 +97,24 @@ class TestGeom: reference = os.path.join(reference_dir,'scale_{}.geom'.format(tag)) if update: modified.to_file(reference) assert geom_equal(modified,Geom.from_file(reference)) + + @pytest.mark.parametrize('periodic',[(True),(False)]) + def test_tessellation(self,periodic): + grid = np.random.randint(10,20,3) + size = np.random.random(3) + 1.0 + N_seeds= np.random.randint(10,30) + seeds = np.random.rand(N_seeds,3) * np.broadcast_to(size,(N_seeds,3)) + Voronoi = Geom.from_Voronoi_tessellation( grid,size,seeds, periodic) + Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,np.ones(N_seeds),periodic) + assert geom_equal(Laguerre,Voronoi) + + def test_Laguerre(self): + grid = np.random.randint(10,20,3) + size = np.random.random(3) + 1.0 + N_seeds= np.random.randint(10,30) + seeds = np.random.rand(N_seeds,3) * np.broadcast_to(size,(N_seeds,3)) + weights= np.full((N_seeds),-np.inf) + ms = np.random.randint(1, N_seeds+1) + weights[ms-1] = np.random.random() + Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,weights,np.random.random()>0.5) + assert np.all(Laguerre.microstructure == ms) From 7b7ac294ca5d19f9c427b16a3d4181a7f1b61484 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 29 Mar 2020 19:13:29 +0200 Subject: [PATCH 04/35] volatile seems to make sense here the value can be changed surprisingly --- src/DAMASK_interface.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/DAMASK_interface.f90 b/src/DAMASK_interface.f90 index 6d5cbbdb6..04106d040 100644 --- a/src/DAMASK_interface.f90 +++ b/src/DAMASK_interface.f90 @@ -25,7 +25,7 @@ module DAMASK_interface implicit none private - logical, public, protected :: & + logical, volatile, public, protected :: & SIGTERM, & !< termination signal SIGUSR1, & !< 1. user-defined signal SIGUSR2 !< 2. user-defined signal From b6596a0310b97b6477ca956ec0fd144ac2f12928 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 29 Mar 2020 19:17:24 +0200 Subject: [PATCH 05/35] compiler can do the counting --- src/lattice.f90 | 22 +++++++++++----------- src/math.f90 | 16 +++++++--------- 2 files changed, 18 insertions(+), 20 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 31a73d7d6..120c58a15 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -18,16 +18,16 @@ module lattice !-------------------------------------------------------------------------------------------------- ! face centered cubic - integer, dimension(2), parameter :: & + integer, dimension(*), parameter :: & FCC_NSLIPSYSTEM = [12, 6] !< # of slip systems per family for fcc - integer, dimension(1), parameter :: & + integer, dimension(*), parameter :: & FCC_NTWINSYSTEM = [12] !< # of twin systems per family for fcc - integer, dimension(1), parameter :: & + integer, dimension(*), parameter :: & FCC_NTRANSSYSTEM = [12] !< # of transformation systems per family for fcc - integer, dimension(1), parameter :: & + integer, dimension(*), parameter :: & FCC_NCLEAVAGESYSTEM = [3] !< # of cleavage systems per family for fcc integer, parameter :: & @@ -109,13 +109,13 @@ module lattice !-------------------------------------------------------------------------------------------------- ! body centered cubic - integer, dimension(2), parameter :: & + integer, dimension(*), parameter :: & BCC_NSLIPSYSTEM = [12, 12] !< # of slip systems per family for bcc - integer, dimension(1), parameter :: & + integer, dimension(*), parameter :: & BCC_NTWINSYSTEM = [12] !< # of twin systems per family for bcc - integer, dimension(1), parameter :: & + integer, dimension(*), parameter :: & BCC_NCLEAVAGESYSTEM = [3] !< # of cleavage systems per family for bcc integer, parameter :: & @@ -187,10 +187,10 @@ module lattice !-------------------------------------------------------------------------------------------------- ! hexagonal - integer, dimension(6), parameter :: & + integer, dimension(*), parameter :: & HEX_NSLIPSYSTEM = [3, 3, 3, 6, 12, 6] !< # of slip systems per family for hex - integer, dimension(4), parameter :: & + integer, dimension(*), parameter :: & HEX_NTWINSYSTEM = [6, 6, 6, 6] !< # of slip systems per family for hex integer, parameter :: & @@ -280,7 +280,7 @@ module lattice !-------------------------------------------------------------------------------------------------- ! body centered tetragonal - integer, dimension(13), parameter :: & + integer, dimension(*), parameter :: & BCT_NSLIPSYSTEM = [2, 2, 2, 4, 2, 4, 2, 2, 4, 8, 4, 8, 8 ] !< # of slip systems per family for bct (Sn) Bieler J. Electr Mater 2009 integer, parameter :: & @@ -362,7 +362,7 @@ module lattice !-------------------------------------------------------------------------------------------------- ! orthorhombic - integer, dimension(3), parameter :: & + integer, dimension(*), parameter :: & ORT_NCLEAVAGESYSTEM = [1, 1, 1] !< # of cleavage systems per family for ortho integer, parameter :: & diff --git a/src/math.f90 b/src/math.f90 index a296dce31..b48037c65 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -30,14 +30,12 @@ module math 1.0_pReal,0.0_pReal,0.0_pReal, & 0.0_pReal,1.0_pReal,0.0_pReal, & 0.0_pReal,0.0_pReal,1.0_pReal & - ],[3,3]) !< 3x3 Identity + ],shape(math_I3)) !< 3x3 Identity - real(pReal), dimension(6), parameter, private :: & - NRMMANDEL = [& - 1.0_pReal, 1.0_pReal, 1.0_pReal, & - sqrt(2.0_pReal), sqrt(2.0_pReal), sqrt(2.0_pReal) ] !< forward weighting for Mandel notation + real(pReal), dimension(*), parameter, private :: & + NRMMANDEL = [1.0_pReal, 1.0_pReal,1.0_pReal, sqrt(2.0_pReal), sqrt(2.0_pReal), sqrt(2.0_pReal)] !< forward weighting for Mandel notation - real(pReal), dimension(6), parameter, private :: & + real(pReal), dimension(*), parameter, private :: & INVNRMMANDEL = 1.0_pReal/NRMMANDEL !< backward weighting for Mandel notation integer, dimension (2,6), parameter, private :: & @@ -48,7 +46,7 @@ module math 1,2, & 2,3, & 1,3 & - ],[2,6]) !< arrangement in Nye notation. + ],shape(MAPNYE)) !< arrangement in Nye notation. integer, dimension (2,6), parameter, private :: & MAPVOIGT = reshape([& @@ -58,7 +56,7 @@ module math 2,3, & 1,3, & 1,2 & - ],[2,6]) !< arrangement in Voigt notation + ],shape(MAPVOIGT)) !< arrangement in Voigt notation integer, dimension (2,9), parameter, private :: & MAPPLAIN = reshape([& @@ -71,7 +69,7 @@ module math 3,1, & 3,2, & 3,3 & - ],[2,9]) !< arrangement in Plain notation + ],shape(MAPPLAIN)) !< arrangement in Plain notation interface math_eye module procedure math_identity2nd From 08ad9d1d571f26d920c6d3b777a7b0e0150fdb4d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 29 Mar 2020 19:18:06 +0200 Subject: [PATCH 06/35] was mixed up --- src/grid/spectral_utilities.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 55375ee32..66426fc9a 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -1056,14 +1056,15 @@ subroutine utilities_updateCoords(F) call MPI_Irecv(IPfluct_padded(:,:,:,grid3+2),c,MPI_DOUBLE,rank_t,0,PETSC_COMM_WORLD,r,ierr) if(ierr /=0) call IO_error(894, ext_msg='update_IPcoords/MPI_Irecv') call MPI_Wait(r,s,ierr) + if(ierr /=0) call IO_error(894, ext_msg='update_IPcoords/MPI_Wait') ! send top layer to process above - if(ierr /=0) call IO_error(894, ext_msg='update_IPcoords/MPI_Wait') call MPI_Isend(IPfluct_padded(:,:,:,grid3+1),c,MPI_DOUBLE,rank_t,0,PETSC_COMM_WORLD,r,ierr) if(ierr /=0) call IO_error(894, ext_msg='update_IPcoords/MPI_Isend') call MPI_Irecv(IPfluct_padded(:,:,:,1), c,MPI_DOUBLE,rank_b,0,PETSC_COMM_WORLD,r,ierr) if(ierr /=0) call IO_error(894, ext_msg='update_IPcoords/MPI_Irecv') call MPI_Wait(r,s,ierr) + if(ierr /=0) call IO_error(894, ext_msg='update_IPcoords/MPI_Wait') !-------------------------------------------------------------------------------------------------- ! calculate nodal displacements From ba5538516c80737545ae3c87d56a036edd10a3c4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 29 Mar 2020 20:04:24 +0200 Subject: [PATCH 07/35] simplified --- src/math.f90 | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index b48037c65..b0852e8d4 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -485,21 +485,19 @@ function math_invSym3333(A) real(pReal),dimension(3,3,3,3),intent(in) :: A - integer :: ierr integer, dimension(6) :: ipiv6 real(pReal), dimension(6,6) :: temp66 real(pReal), dimension(6*(64+2)) :: work - logical :: error + integer :: ierr_i, ierr_f external :: & dgetrf, & dgetri temp66 = math_sym3333to66(A) - call dgetrf(6,6,temp66,6,ipiv6,ierr) - error = (ierr /= 0) - call dgetri(6,temp66,6,ipiv6,work,size(work,1),ierr) - error = error .or. (ierr /= 0) - if (error) then + call dgetrf(6,6,temp66,6,ipiv6,ierr_i) + call dgetri(6,temp66,6,ipiv6,work,size(work,1),ierr_f) + + if (ierr_i /= 0 .or. ierr_f /= 0) then call IO_error(400, ext_msg = 'math_invSym3333') else math_invSym3333 = math_66toSym3333(temp66) From 2a37acfe5eaa8ff83d61bac797bd15e36b0f48f3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 29 Mar 2020 20:04:51 +0200 Subject: [PATCH 08/35] store data where it is needed --- src/DAMASK_marc.f90 | 5 ++++- src/FEsolving.f90 | 10 +--------- src/homogenization.f90 | 2 ++ 3 files changed, 7 insertions(+), 10 deletions(-) diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index fe01301c9..9fd41459b 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -43,6 +43,9 @@ module DAMASK_interface logical, protected, public :: symmetricSolver character(len=*), parameter, public :: INPUTFILEEXTENSION = '.dat' + logical, dimension(:,:), public, allocatable :: & + calcMode !< calculate or collect (ping pong scheme) + public :: & DAMASK_interface_init, & getSolverJobName @@ -102,7 +105,6 @@ function getSolverJobName() character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forward and backward slash integer :: extPos - inputName='' inquire(5, name=inputName) ! determine inputfile extPos = len_trim(inputName)-4 getSolverJobName=inputName(scan(inputName,pathSep,back=.true.)+1:extPos) @@ -179,6 +181,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & use FEsolving use debug use discretization_marc + use homogenization use CPFEM implicit none diff --git a/src/FEsolving.f90 b/src/FEsolving.f90 index cce95a07e..3fc1482d3 100644 --- a/src/FEsolving.f90 +++ b/src/FEsolving.f90 @@ -4,20 +4,12 @@ !> @brief global variables for flow control !-------------------------------------------------------------------------------------------------- module FEsolving - + implicit none public - logical :: & - terminallyIll = .false. !< at least one material point is terminally ill - integer, dimension(2) :: & FEsolving_execElem, & !< for ping-pong scheme always whole range, otherwise one specific element FEsolving_execIP !< for ping-pong scheme always range to max IP, otherwise one specific IP - -#if defined(Marc4DAMASK) - logical, dimension(:,:), allocatable :: & - calcMode !< do calculation or simply collect when using ping pong scheme -#endif end module FEsolving diff --git a/src/homogenization.f90 b/src/homogenization.f90 index b657dfe3e..a9b173831 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -29,6 +29,8 @@ module homogenization !-------------------------------------------------------------------------------------------------- ! General variables for the homogenization at a material point + logical, public :: & + terminallyIll = .false. !< at least one material point is terminally ill real(pReal), dimension(:,:,:,:), allocatable, public :: & materialpoint_F0, & !< def grad of IP at start of FE increment materialpoint_F, & !< def grad of IP to be reached at end of FE increment From 9c90aa5acb622a38c5bf587648066819a2d5f8b4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 29 Mar 2020 20:07:09 +0200 Subject: [PATCH 09/35] polishing --- python/tests/test_Geom.py | 2 +- src/constitutive_plastic_nonlocal.f90 | 14 +++++++------- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/python/tests/test_Geom.py b/python/tests/test_Geom.py index 1b1529e80..c098603c6 100644 --- a/python/tests/test_Geom.py +++ b/python/tests/test_Geom.py @@ -28,7 +28,7 @@ def reference_dir(reference_dir_base): class TestGeom: - + def test_update(self,default): modified = copy.deepcopy(default) modified.update( diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index 58218ac00..f872e3846 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -16,15 +16,15 @@ submodule(constitutive) plastic_nonlocal kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin ! storage order of dislocation types - integer, dimension(8), parameter :: & + integer, dimension(*), parameter :: & sgl = [1,2,3,4,5,6,7,8] !< signed (single) - integer, dimension(5), parameter :: & + integer, dimension(*), parameter :: & edg = [1,2,5,6,9], & !< edge scr = [3,4,7,8,10] !< screw - integer, dimension(4), parameter :: & + integer, dimension(*), parameter :: & mob = [1,2,3,4], & !< mobile imm = [5,6,7,8] !< immobile (blocked) - integer, dimension(2), parameter :: & + integer, dimension(*), parameter :: & dip = [9,10], & !< dipole imm_edg = imm(1:2), & !< immobile edge imm_scr = imm(3:4) !< immobile screw @@ -1611,7 +1611,7 @@ end subroutine stateInit !-------------------------------------------------------------------------------------------------- !> @brief calculates kinetics !-------------------------------------------------------------------------------------------------- -subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Temperature, instance) +pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Temperature, instance) integer, intent(in) :: & c, & !< dislocation character (1:edge, 2:screw) @@ -1726,7 +1726,7 @@ end subroutine kinetics !> @brief returns copy of current dislocation densities from state !> @details raw values is rectified !-------------------------------------------------------------------------------------------------- -function getRho(instance,of,ip,el) +pure function getRho(instance,of,ip,el) integer, intent(in) :: instance, of,ip,el real(pReal), dimension(param(instance)%sum_N_sl,10) :: getRho @@ -1751,7 +1751,7 @@ end function getRho !> @brief returns copy of current dislocation densities from state !> @details raw values is rectified !-------------------------------------------------------------------------------------------------- -function getRho0(instance,of,ip,el) +pure function getRho0(instance,of,ip,el) integer, intent(in) :: instance, of,ip,el real(pReal), dimension(param(instance)%sum_N_sl,10) :: getRho0 From e334674a06bc5c8667264634f1d4efaf31c15ecb Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 29 Mar 2020 20:28:54 +0200 Subject: [PATCH 10/35] test bycristal tessellation --- python/tests/test_Geom.py | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/python/tests/test_Geom.py b/python/tests/test_Geom.py index c098603c6..f988179fa 100644 --- a/python/tests/test_Geom.py +++ b/python/tests/test_Geom.py @@ -99,7 +99,7 @@ class TestGeom: assert geom_equal(modified,Geom.from_file(reference)) @pytest.mark.parametrize('periodic',[(True),(False)]) - def test_tessellation(self,periodic): + def test_tessellation_approaches(self,periodic): grid = np.random.randint(10,20,3) size = np.random.random(3) + 1.0 N_seeds= np.random.randint(10,30) @@ -108,7 +108,7 @@ class TestGeom: Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,np.ones(N_seeds),periodic) assert geom_equal(Laguerre,Voronoi) - def test_Laguerre(self): + def test_Laguerre_weights(self): grid = np.random.randint(10,20,3) size = np.random.random(3) + 1.0 N_seeds= np.random.randint(10,30) @@ -118,3 +118,16 @@ class TestGeom: weights[ms-1] = np.random.random() Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,weights,np.random.random()>0.5) assert np.all(Laguerre.microstructure == ms) + + @pytest.mark.parametrize('approach',[('Laguerre'),('Voronoi')]) + def test_tessellate_bicrystal(self,approach): + grid = np.random.randint(5,10,3)*2 + size = grid.astype(np.float) + seeds = np.vstack((size*np.array([0.5,0.25,0.5]),size*np.array([0.5,0.75,0.5]))) + microstructure = np.ones(grid) + microstructure[:,grid[1]//2:,:] = 2 + if approach == 'Laguerre': + geom = Geom.from_Laguerre_tessellation(grid,size,seeds,np.ones(2),np.random.random()>0.5) + elif approach == 'Voronoi': + geom = Geom.from_Voronoi_tessellation(grid,size,seeds, np.random.random()>0.5) + assert np.all(geom.microstructure == microstructure) From 11e58bcc2f44b86894de7a6104b147336aa64be3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 29 Mar 2020 22:37:48 +0200 Subject: [PATCH 11/35] not needed why debugging allocate? --- src/homogenization.f90 | 16 ---------------- 1 file changed, 16 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index a9b173831..50331cbdd 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -173,22 +173,6 @@ subroutine homogenization_init write(6,'(/,a)') ' <<<+- homogenization init -+>>>'; flush(6) - if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0) then - write(6,'(a32,1x,7(i8,1x))') 'materialpoint_dPdF: ', shape(materialpoint_dPdF) - write(6,'(a32,1x,7(i8,1x))') 'materialpoint_F0: ', shape(materialpoint_F0) - write(6,'(a32,1x,7(i8,1x))') 'materialpoint_F: ', shape(materialpoint_F) - write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subF0: ', shape(materialpoint_subF0) - write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subF: ', shape(materialpoint_subF) - write(6,'(a32,1x,7(i8,1x))') 'materialpoint_P: ', shape(materialpoint_P) - write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subFrac: ', shape(materialpoint_subFrac) - write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subStep: ', shape(materialpoint_subStep) - write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subdt: ', shape(materialpoint_subdt) - write(6,'(a32,1x,7(i8,1x))') 'materialpoint_requested: ', shape(materialpoint_requested) - write(6,'(a32,1x,7(i8,1x))') 'materialpoint_converged: ', shape(materialpoint_converged) - write(6,'(a32,1x,7(i8,1x),/)') 'materialpoint_doneAndHappy: ', shape(materialpoint_doneAndHappy) - endif - flush(6) - if (debug_g < 1 .or. debug_g > homogenization_Ngrains(material_homogenizationAt(debug_e))) & call IO_error(602,ext_msg='constituent', el=debug_e, g=debug_g) From 7eeb5db15ff9c3cdadbd4285eaabb14553417919 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 29 Mar 2020 22:41:47 +0200 Subject: [PATCH 12/35] intention clearer slip/twin/trans happens in untwinned/untransformed volume, shear banding is independent of that --- src/constitutive_plastic_dislotwin.f90 | 34 ++++++++++++-------------- 1 file changed, 16 insertions(+), 18 deletions(-) diff --git a/src/constitutive_plastic_dislotwin.f90 b/src/constitutive_plastic_dislotwin.f90 index d36e08846..21fe555f2 100644 --- a/src/constitutive_plastic_dislotwin.f90 +++ b/src/constitutive_plastic_dislotwin.f90 @@ -568,7 +568,22 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) + ddot_gamma_dtau_slip(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i) enddo slipContribution - !ToDo: Why do this before shear banding? + call kinetics_twin(Mp,T,dot_gamma_sl,instance,of,dot_gamma_twin,ddot_gamma_dtau_twin) + twinContibution: do i = 1, prm%sum_N_tw + Lp = Lp + dot_gamma_twin(i)*prm%P_tw(1:3,1:3,i) + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + ddot_gamma_dtau_twin(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) + enddo twinContibution + + call kinetics_trans(Mp,T,dot_gamma_sl,instance,of,dot_gamma_tr,ddot_gamma_dtau_trans) + transContibution: do i = 1, prm%sum_N_tr + Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i) + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + ddot_gamma_dtau_trans(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i) + enddo transContibution + Lp = Lp * f_unrotated dLp_dMp = dLp_dMp * f_unrotated @@ -598,23 +613,6 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) endif shearBandingContribution - call kinetics_twin(Mp,T,dot_gamma_sl,instance,of,dot_gamma_twin,ddot_gamma_dtau_twin) - twinContibution: do i = 1, prm%sum_N_tw - Lp = Lp + dot_gamma_twin(i)*prm%P_tw(1:3,1:3,i) * f_unrotated - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + ddot_gamma_dtau_twin(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) * f_unrotated - enddo twinContibution - - call kinetics_trans(Mp,T,dot_gamma_sl,instance,of,dot_gamma_tr,ddot_gamma_dtau_trans) - transContibution: do i = 1, prm%sum_N_tr - Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i) * f_unrotated - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + ddot_gamma_dtau_trans(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i) * f_unrotated - enddo transContibution - - end associate end subroutine plastic_dislotwin_LpAndItsTangent From 396d428af73f4d1fb50c462700ef2deb8cfe5e44 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 29 Mar 2020 23:09:43 +0200 Subject: [PATCH 13/35] bugfix: works for all cuboids, not just cubes --- processing/pre/geom_fromVoronoiTessellation.py | 4 ++-- python/damask/_geom.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/processing/pre/geom_fromVoronoiTessellation.py b/processing/pre/geom_fromVoronoiTessellation.py index edcbe83dc..31e8caa7b 100755 --- a/processing/pre/geom_fromVoronoiTessellation.py +++ b/processing/pre/geom_fromVoronoiTessellation.py @@ -44,8 +44,8 @@ def Laguerre_tessellation(grid, size, seeds, weights, origin = np.zeros(3), peri closest_seed= np.array([findClosestSeed(seeds_p,weights_p,coord) for coord in coords]) if periodic: - closest_seed = closest_seed.reshape(grid[2]*3,grid[1]*3,grid[0]*3) - return closest_seed[grid[2]:grid[2]*2,grid[1]:grid[1]*2,grid[0]:grid[0]*2]%seeds.shape[0] + closest_seed = closest_seed.reshape(grid*3) + return closest_seed[grid[0]:grid[0]*2,grid[1]:grid[1]*2,grid[2]:grid[2]*2]%seeds.shape[0] else: return closest_seed diff --git a/python/damask/_geom.py b/python/damask/_geom.py index 797b9dc81..7288be6cb 100644 --- a/python/damask/_geom.py +++ b/python/damask/_geom.py @@ -370,7 +370,7 @@ class Geom: microstructure = np.array(result.get()) if periodic: - microstructure = microstructure.reshape(grid[0]*3,grid[1]*3,grid[2]*3) + microstructure = microstructure.reshape(grid*3) microstructure = microstructure[grid[0]:grid[0]*2,grid[1]:grid[1]*2,grid[2]:grid[2]*2]%seeds.shape[0] else: microstructure = microstructure.reshape(grid) From bcbdd87870229921eb42b785d8dd367bce6d293d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 29 Mar 2020 23:50:09 +0200 Subject: [PATCH 14/35] base substitution on original microstructure --- processing/pre/geom_translate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/processing/pre/geom_translate.py b/processing/pre/geom_translate.py index 2d4279821..1d08c23dc 100755 --- a/processing/pre/geom_translate.py +++ b/processing/pre/geom_translate.py @@ -52,7 +52,7 @@ for name in filenames: geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name) substituted = geom.get_microstructure() - for old,new in zip(sub[0::2],sub[1::2]): substituted[substituted==old] = new # substitute microstructure indices + for old,new in zip(sub[0::2],sub[1::2]): substituted[geom.microstructure==old] = new # substitute microstructure indices substituted += options.microstructure # constant shift damask.util.croak(geom.update(substituted,origin=geom.get_origin()+options.origin)) From 2f52165e0a604fcdfe7575dd41b41ad5e02bd028 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 30 Mar 2020 22:14:45 +0200 Subject: [PATCH 15/35] for test coverage reports using pytest --- python/.coveragerc | 2 ++ python/.gitignore | 1 + 2 files changed, 3 insertions(+) create mode 100644 python/.coveragerc diff --git a/python/.coveragerc b/python/.coveragerc new file mode 100644 index 000000000..c712d2595 --- /dev/null +++ b/python/.coveragerc @@ -0,0 +1,2 @@ +[run] +omit = tests/* diff --git a/python/.gitignore b/python/.gitignore index 1cb7eb8a5..7c60e1af8 100644 --- a/python/.gitignore +++ b/python/.gitignore @@ -1,3 +1,4 @@ *.pyc dist damask.egg-info +.coverage From 18f60a94a9680878930f073fc27efad79f5268a5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 31 Mar 2020 09:04:08 +0200 Subject: [PATCH 16/35] one variable suffices --- src/crystallite.f90 | 30 +++++++++++++----------------- 1 file changed, 13 insertions(+), 17 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 8a5a89d62..ece4db020 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1014,15 +1014,13 @@ subroutine integrateStateFPI sizeDotState real(pReal) :: & zeta - real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: & - residuum_plastic ! residuum for plastic state - real(pReal), dimension(constitutive_source_maxSizeDotState) :: & - residuum_source ! residuum for source state + real(pReal), dimension(max(constitutive_plasticity_maxSizeDotState,constitutive_source_maxSizeDotState)) :: & + r ! state residuum logical :: & nonlocalBroken nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE(sizeDotState,residuum_plastic,residuum_source,zeta,p,c) + !$OMP PARALLEL DO PRIVATE(sizeDotState,r,zeta,p,c) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) @@ -1095,13 +1093,12 @@ subroutine integrateStateFPI plasticState(p)%previousDotState2(:,c)) plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta & + plasticState(p)%previousDotState(:,c) * (1.0_pReal - zeta) - residuum_plastic(1:SizeDotState) = plasticState(p)%state (1:sizeDotState,c) & - - plasticState(p)%subState0(1:sizeDotState,c) & - - plasticState(p)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) + r(1:SizeDotState) = plasticState(p)%state (1:sizeDotState,c) & + - plasticState(p)%subState0(1:sizeDotState,c) & + - plasticState(p)%dotState (1:sizeDotState,c) * crystallite_subdt(g,i,e) plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%state(1:sizeDotState,c) & - - residuum_plastic(1:sizeDotState) - crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState), & + - r(1:sizeDotState) + crystallite_converged(g,i,e) = converged(r(1:sizeDotState), & plasticState(p)%state(1:sizeDotState,c), & plasticState(p)%atol(1:sizeDotState)) do s = 1, phase_Nsources(p) @@ -1111,14 +1108,13 @@ subroutine integrateStateFPI sourceState(p)%p(s)%previousDotState2(:,c)) sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) * zeta & + sourceState(p)%p(s)%previousDotState(:,c)* (1.0_pReal - zeta) - residuum_source(1: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) + r(1: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) sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%state(1:sizeDotState,c) & - - residuum_source(1:sizeDotState) + - r(1:sizeDotState) crystallite_converged(g,i,e) = & - crystallite_converged(g,i,e) .and. converged(residuum_source(1:sizeDotState), & + crystallite_converged(g,i,e) .and. converged(r(1:sizeDotState), & sourceState(p)%p(s)%state(1:sizeDotState,c), & sourceState(p)%p(s)%atol(1:sizeDotState)) enddo From 604bcd1229b2bdda91c89c998b457f3341471342 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 31 Mar 2020 09:38:08 +0200 Subject: [PATCH 17/35] memory efficient FPI state integrator --- src/crystallite.f90 | 30 +++++++++++++----------------- src/material.f90 | 8 -------- 2 files changed, 13 insertions(+), 25 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index ece4db020..26fa1b3f8 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1016,11 +1016,13 @@ subroutine integrateStateFPI zeta real(pReal), dimension(max(constitutive_plasticity_maxSizeDotState,constitutive_source_maxSizeDotState)) :: & r ! state residuum + real(pReal), dimension(:), allocatable :: plastic_dotState_p1, plastic_dotState_p2 + type(group_float), dimension(maxval(phase_Nsources)) :: source_dotState_p1, source_dotState_p2 logical :: & nonlocalBroken nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE(sizeDotState,r,zeta,p,c) + !$OMP PARALLEL DO PRIVATE(sizeDotState,r,zeta,p,c,plastic_dotState_p1, plastic_dotState_p2,source_dotState_p1, source_dotState_p2) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) @@ -1045,24 +1047,22 @@ subroutine integrateStateFPI plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & + plasticState(p)%dotState (1:sizeDotState,c) & * crystallite_subdt(g,i,e) + plastic_dotState_p2 = 0.0_pReal * plasticState(p)%dotState (1:sizeDotState,c) ! ToDo can be done smarter/clearer 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) + source_dotState_p2(p)%p = 0.0_pReal * sourceState(p)%p(s)%dotState (1:sizeDotState,c) ! ToDo can be done smarter/clearer enddo iteration: do NiterationState = 1, num%nState - plasticState(p)%previousDotState2(:,c) = merge(plasticState(p)%previousDotState(:,c),& - 0.0_pReal,& - NiterationState > 1) - plasticState(p)%previousDotState (:,c) = plasticState(p)%dotState(:,c) + if(nIterationState > 1) plastic_dotState_p2 = plastic_dotState_p1 + plastic_dotState_p1 = plasticState(p)%dotState(:,c) do s = 1, phase_Nsources(p) - sourceState(p)%p(s)%previousDotState2(:,c) = merge(sourceState(p)%p(s)%previousDotState(:,c),& - 0.0_pReal, & - NiterationState > 1) - sourceState(p)%p(s)%previousDotState (:,c) = sourceState(p)%p(s)%dotState(:,c) + if(nIterationState > 1) source_dotState_p2(p)%p = source_dotState_p1(p)%p + source_dotState_p1(p)%p = sourceState(p)%p(s)%dotState(:,c) enddo call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & @@ -1088,11 +1088,9 @@ subroutine integrateStateFPI if(.not. crystallite_todo(g,i,e)) exit iteration sizeDotState = plasticState(p)%sizeDotState - zeta = damper(plasticState(p)%dotState (:,c), & - plasticState(p)%previousDotState (:,c), & - plasticState(p)%previousDotState2(:,c)) + zeta = damper(plasticState(p)%dotState(:,c),plastic_dotState_p1,plastic_dotState_p2) plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta & - + plasticState(p)%previousDotState(:,c) * (1.0_pReal - zeta) + + plastic_dotState_p1 * (1.0_pReal - zeta) r(1:SizeDotState) = plasticState(p)%state (1:sizeDotState,c) & - plasticState(p)%subState0(1:sizeDotState,c) & - plasticState(p)%dotState (1:sizeDotState,c) * crystallite_subdt(g,i,e) @@ -1103,11 +1101,9 @@ subroutine integrateStateFPI plasticState(p)%atol(1:sizeDotState)) do s = 1, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - zeta = damper(sourceState(p)%p(s)%dotState (:,c), & - sourceState(p)%p(s)%previousDotState (:,c), & - sourceState(p)%p(s)%previousDotState2(:,c)) + zeta = damper(sourceState(p)%p(s)%dotState(:,c),source_dotState_p1(p)%p,source_dotState_p2(p)%p) sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) * zeta & - + sourceState(p)%p(s)%previousDotState(:,c)* (1.0_pReal - zeta) + + source_dotState_p1(p)%p* (1.0_pReal - zeta) r(1: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) diff --git a/src/material.f90 b/src/material.f90 index 4461ca958..14d4e4c19 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -724,10 +724,6 @@ subroutine material_allocatePlasticState(phase,NipcMyPhase,& allocate(plasticState(phase)%state (sizeState,NipcMyPhase), source=0.0_pReal) allocate(plasticState(phase)%dotState (sizeDotState,NipcMyPhase),source=0.0_pReal) - if (numerics_integrator == 1) then - allocate(plasticState(phase)%previousDotState (sizeDotState,NipcMyPhase),source=0.0_pReal) - allocate(plasticState(phase)%previousDotState2 (sizeDotState,NipcMyPhase),source=0.0_pReal) - endif if (numerics_integrator == 4) & allocate(plasticState(phase)%RK4dotState (4,sizeDotState,NipcMyPhase),source=0.0_pReal) if (numerics_integrator == 5) & @@ -762,10 +758,6 @@ subroutine material_allocateSourceState(phase,of,NipcMyPhase,& allocate(sourceState(phase)%p(of)%state (sizeState,NipcMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(of)%dotState (sizeDotState,NipcMyPhase),source=0.0_pReal) - if (numerics_integrator == 1) then - allocate(sourceState(phase)%p(of)%previousDotState (sizeDotState,NipcMyPhase),source=0.0_pReal) - allocate(sourceState(phase)%p(of)%previousDotState2 (sizeDotState,NipcMyPhase),source=0.0_pReal) - endif if (numerics_integrator == 4) & allocate(sourceState(phase)%p(of)%RK4dotState (4,sizeDotState,NipcMyPhase),source=0.0_pReal) if (numerics_integrator == 5) & From 28dadb4422d3898d0ba2420aa946d8853dcf330d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 31 Mar 2020 09:49:41 +0200 Subject: [PATCH 18/35] no need to check multiple times --- src/crystallite.f90 | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 26fa1b3f8..5c6510864 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1070,8 +1070,6 @@ subroutine integrateStateFPI g, i, e) crystallite_todo(g,i,e) = integrateStress(g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & - nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) exit iteration call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & @@ -1083,8 +1081,6 @@ subroutine integrateStateFPI do s = 1, phase_Nsources(p) crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c))) enddo - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & - nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) exit iteration sizeDotState = plasticState(p)%sizeDotState @@ -1117,12 +1113,12 @@ subroutine integrateStateFPI if(crystallite_converged(g,i,e)) then crystallite_todo(g,i,e) = stateJump(g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & - nonlocalBroken = .true. exit iteration endif enddo iteration + if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + nonlocalBroken = .true. endif enddo; enddo; enddo From 2cc0c746d31f3ed784c8664cff3389e45da89344 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 31 Mar 2020 10:07:01 +0200 Subject: [PATCH 19/35] no variable needed --- src/crystallite.f90 | 19 ++++--------------- 1 file changed, 4 insertions(+), 15 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 5c6510864..9c4eeb326 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1540,15 +1540,6 @@ subroutine integrateStateRKCK45 logical :: & nonlocalBroken - real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & - homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & - residuum_plastic - real(pReal), dimension(constitutive_source_maxSizeDotState, & - maxval(phase_Nsources), & - homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & - residuum_source - - nonlocalBroken = .false. !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c) do e = FEsolving_execElem(1),FEsolving_execElem(2) @@ -1629,13 +1620,12 @@ subroutine integrateStateRKCK45 sizeDotState = plasticState(p)%sizeDotState plasticState(p)%RKCK45dotState(6,:,c) = plasticState (p)%dotState(:,c) - residuum_plastic(1:sizeDotState,g,i,e) = matmul(DB,plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,c)) & - * crystallite_subdt(g,i,e) plasticState(p)%dotState(:,c) = matmul(B,plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,c)) plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & + plasticState(p)%dotState (1:sizeDotState,c) & * crystallite_subdt(g,i,e) - crystallite_todo(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), & + crystallite_todo(g,i,e) = converged(matmul(DB,plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,c)) & + * crystallite_subdt(g,i,e), & plasticState(p)%state(1:sizeDotState,c), & plasticState(p)%atol(1:sizeDotState)) @@ -1643,14 +1633,13 @@ subroutine integrateStateRKCK45 sizeDotState = sourceState(p)%p(s)%sizeDotState sourceState(p)%p(s)%RKCK45dotState(6,:,c) = sourceState(p)%p(s)%dotState(:,c) - residuum_source(1:sizeDotState,s,g,i,e) = matmul(DB,sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c)) & - * crystallite_subdt(g,i,e) sourceState(p)%p(s)%dotState(:,c) = matmul(B,sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c)) 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) crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. & - converged(residuum_source(1:sizeDotState,s,g,i,e), & + converged(matmul(DB,sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c)) & + * crystallite_subdt(g,i,e), & sourceState(p)%p(s)%state(1:sizeDotState,c), & sourceState(p)%p(s)%atol(1:sizeDotState)) enddo From 505c1432b1f61bc0718f7780110b722d7a4a2684 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 31 Mar 2020 10:12:25 +0200 Subject: [PATCH 20/35] no need for temp storage --- src/crystallite.f90 | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 9c4eeb326..52d102536 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -804,7 +804,6 @@ logical function integrateStress(ipc,ip,el,timeFraction) residuumLi_old, & ! last residuum of intermediate velocity gradient deltaLi, & ! direction of next guess Fe, & ! elastic deformation gradient - Fe_new, & S, & ! 2nd Piola-Kirchhoff Stress in plastic (lattice) configuration A, & B, & @@ -982,17 +981,15 @@ logical function integrateStress(ipc,ip,el,timeFraction) invFp_new = matmul(invFp_current,B) call math_invert33(Fp_new,devNull,error,invFp_new) if (error) return ! error - Fp_new = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal) ! regularize - Fe_new = matmul(matmul(F,invFp_new),invFi_new) integrateStress = .true. crystallite_P (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new))) crystallite_S (1:3,1:3,ipc,ip,el) = S crystallite_Lp (1:3,1:3,ipc,ip,el) = Lpguess crystallite_Li (1:3,1:3,ipc,ip,el) = Liguess - crystallite_Fp (1:3,1:3,ipc,ip,el) = Fp_new + crystallite_Fp (1:3,1:3,ipc,ip,el) = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal) ! regularize crystallite_Fi (1:3,1:3,ipc,ip,el) = Fi_new - crystallite_Fe (1:3,1:3,ipc,ip,el) = Fe_new + crystallite_Fe (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),invFi_new) end function integrateStress From 5b29af847349b8df64e40c1a5975be7092c004f8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 31 Mar 2020 10:16:14 +0200 Subject: [PATCH 21/35] clearer --- src/crystallite.f90 | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 52d102536..4df98c429 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -910,21 +910,19 @@ logical function integrateStress(ipc,ip,el,timeFraction) cycle LpLoop endif - !* calculate Jacobian for correction term - if (mod(jacoCounterLp, num%iJacoLpresiduum) == 0) then + calculateJacobiLi: if (mod(jacoCounterLp, num%iJacoLpresiduum) == 0) then jacoCounterLp = jacoCounterLp + 1 do o=1,3; do p=1,3 - dFe_dLp(o,1:3,p,1:3) = A(o,p)*transpose(invFi_new) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j) + dFe_dLp(o,1:3,p,1:3) = - dt * A(o,p)*transpose(invFi_new) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j) enddo; enddo - dFe_dLp = - dt * dFe_dLp dRLp_dLp = math_identity2nd(9) & - math_3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dS,dS_dFe),dFe_dLp)) temp_9 = math_33to9(residuumLp) call dgesv(9,1,dRLp_dLp,9,devNull_9,temp_9,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp if (ierr /= 0) return ! error deltaLp = - math_9to33(temp_9) - endif + endif calculateJacobiLi Lpguess = Lpguess & + deltaLp * steplengthLp @@ -952,8 +950,7 @@ logical function integrateStress(ipc,ip,el,timeFraction) cycle LiLoop endif - !* calculate Jacobian for correction term - if (mod(jacoCounterLi, num%iJacoLpresiduum) == 0) then + calculateJacobiLp: if (mod(jacoCounterLi, num%iJacoLpresiduum) == 0) then jacoCounterLi = jacoCounterLi + 1 temp_33 = matmul(matmul(A,B),invFi_current) @@ -972,7 +969,7 @@ logical function integrateStress(ipc,ip,el,timeFraction) call dgesv(9,1,dRLi_dLi,9,devNull_9,temp_9,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li if (ierr /= 0) return ! error deltaLi = - math_9to33(temp_9) - endif + endif calculateJacobiLp Liguess = Liguess & + deltaLi * steplengthLi From 6ef7410e5a899002ebaf06c6beb2dc8d284c0292 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 31 Mar 2020 11:04:06 +0200 Subject: [PATCH 22/35] testing VTK wrappers --- python/damask/_vtk.py | 2 +- python/tests/test_VTK.py | 48 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 49 insertions(+), 1 deletion(-) create mode 100644 python/tests/test_VTK.py diff --git a/python/damask/_vtk.py b/python/damask/_vtk.py index c11e07b4d..4505a5ebc 100644 --- a/python/damask/_vtk.py +++ b/python/damask/_vtk.py @@ -72,7 +72,7 @@ class VTK: connectivity : numpy.ndarray of np.dtype = int Cell connectivity (0-based), first dimension determines #Cells, second dimension determines #Nodes/Cell. cell_type : str - Name of the vtk.vtkCell subclass. Tested for TRIANGLE, QUAD, and HEXAHEDRON. + Name of the vtk.vtkCell subclass. Tested for TRIANGLE, QUAD, TETRA, and HEXAHEDRON. """ vtk_nodes = vtk.vtkPoints() diff --git a/python/tests/test_VTK.py b/python/tests/test_VTK.py new file mode 100644 index 000000000..ef4090da5 --- /dev/null +++ b/python/tests/test_VTK.py @@ -0,0 +1,48 @@ +import os +import sys + +import pytest +import numpy as np + +from damask import VTK + +@pytest.fixture +def reference_dir(reference_dir_base): + """Directory containing reference results.""" + return os.path.join(reference_dir_base,'Result') + +class TestVTK: + + def test_rectilinearGrid(self,tmp_path): + grid = np.random.randint(5,10,3)*2 + size = np.random.random(3) + 1.0 + origin = np.random.random(3) + v = VTK.from_rectilinearGrid(grid,size,origin) + s = v.__repr__() + v.write(os.path.join(tmp_path,'rectilinearGrid')) + v = VTK.from_file(os.path.join(tmp_path,'rectilinearGrid.vtr')) + assert(s == v.__repr__()) + + def test_polyData(self,tmp_path): + points = np.random.rand(3,100) + v = VTK.from_polyData(points) + s = v.__repr__() + v.write(os.path.join(tmp_path,'polyData')) + v = VTK.from_file(os.path.join(tmp_path,'polyData.vtp')) + assert(s == v.__repr__()) + + @pytest.mark.parametrize('cell_type,n',[ + ('VTK_hexahedron',8), + ('TETRA',4), + ('quad',4), + ('VTK_TRIANGLE',3) + ] + ) + def test_unstructuredGrid(self,tmp_path,cell_type,n): + nodes = np.random.rand(n,3) + connectivity = np.random.choice(np.arange(n),n,False).reshape(-1,n) + v = VTK.from_unstructuredGrid(nodes,connectivity,cell_type) + s = v.__repr__() + v.write(os.path.join(tmp_path,'unstructuredGrid')) + v = VTK.from_file(os.path.join(tmp_path,'unstructuredGrid.vtu')) + assert(s == v.__repr__()) From 2de4b87c6197b8565f83e425e2cf0ee78117cf2c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 31 Mar 2020 11:24:13 +0200 Subject: [PATCH 23/35] bugfix for FPI: loop over sources! no memory waste for adaptive Euler --- src/crystallite.f90 | 39 ++++++++++++++------------------------- 1 file changed, 14 insertions(+), 25 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 4df98c429..ae5091978 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1047,7 +1047,7 @@ subroutine integrateStateFPI 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) - source_dotState_p2(p)%p = 0.0_pReal * sourceState(p)%p(s)%dotState (1:sizeDotState,c) ! ToDo can be done smarter/clearer + source_dotState_p2(s)%p = 0.0_pReal * sourceState(p)%p(s)%dotState (1:sizeDotState,c) ! ToDo can be done smarter/clearer enddo iteration: do NiterationState = 1, num%nState @@ -1055,8 +1055,8 @@ subroutine integrateStateFPI if(nIterationState > 1) plastic_dotState_p2 = plastic_dotState_p1 plastic_dotState_p1 = plasticState(p)%dotState(:,c) do s = 1, phase_Nsources(p) - if(nIterationState > 1) source_dotState_p2(p)%p = source_dotState_p1(p)%p - source_dotState_p1(p)%p = sourceState(p)%p(s)%dotState(:,c) + if(nIterationState > 1) source_dotState_p2(s)%p = source_dotState_p1(s)%p + source_dotState_p1(s)%p = sourceState(p)%p(s)%dotState(:,c) enddo call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & @@ -1091,9 +1091,9 @@ subroutine integrateStateFPI plasticState(p)%atol(1:sizeDotState)) do s = 1, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - zeta = damper(sourceState(p)%p(s)%dotState(:,c),source_dotState_p1(p)%p,source_dotState_p2(p)%p) + zeta = damper(sourceState(p)%p(s)%dotState(:,c),source_dotState_p1(s)%p,source_dotState_p2(s)%p) sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) * zeta & - + source_dotState_p1(p)%p* (1.0_pReal - zeta) + + source_dotState_p1(s)%p* (1.0_pReal - zeta) r(1: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) @@ -1236,14 +1236,8 @@ subroutine integrateStateAdaptiveEuler logical :: & nonlocalBroken - real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & - homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & - residuum_plastic - real(pReal), dimension(constitutive_source_maxSizeDotState,& - maxval(phase_Nsources), & - homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & - residuum_source - + real(pReal), dimension(:), allocatable :: residuum_plastic + type(group_float), dimension(maxval(phase_Nsources)) :: residuum_source nonlocalBroken = .false. !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c) @@ -1269,15 +1263,14 @@ subroutine integrateStateAdaptiveEuler sizeDotState = plasticState(p)%sizeDotState - residuum_plastic(1:sizeDotState,g,i,e) = plasticState(p)%dotstate(1:sizeDotState,c) & - * (- 0.5_pReal * crystallite_subdt(g,i,e)) + residuum_plastic = - 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,g,i,e) = sourceState(p)%p(s)%dotstate(1:sizeDotState,c) & - * (- 0.5_pReal * crystallite_subdt(g,i,e)) + residuum_source(s)%p = - 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 @@ -1312,21 +1305,17 @@ subroutine integrateStateAdaptiveEuler sizeDotState = plasticState(p)%sizeDotState - residuum_plastic(1:sizeDotState,g,i,e) = residuum_plastic(1:sizeDotState,g,i,e) & - + 0.5_pReal * plasticState(p)%dotState(:,c) * crystallite_subdt(g,i,e) - - crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), & + crystallite_converged(g,i,e) = converged(residuum_plastic & + + 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 - residuum_source(1:sizeDotState,s,g,i,e) = & - residuum_source(1:sizeDotState,s,g,i,e) + 0.5_pReal * sourceState(p)%p(s)%dotState(:,c) * crystallite_subdt(g,i,e) - crystallite_converged(g,i,e) = & - crystallite_converged(g,i,e) .and. converged(residuum_source(1:sizeDotState,s,g,i,e), & + crystallite_converged(g,i,e) .and. converged(residuum_source(s)%p & + + 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 From 00be291fa02ce5af06641efa6adc506509ff580c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 31 Mar 2020 12:14:54 +0200 Subject: [PATCH 24/35] todo will be reset after state integration --- src/crystallite.f90 | 17 ++++++----------- 1 file changed, 6 insertions(+), 11 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index ae5091978..0555e28a6 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1118,8 +1118,7 @@ subroutine integrateStateFPI enddo; enddo; enddo !$OMP END PARALLEL DO - if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false. - if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck + if (nonlocalBroken) call nonlocalConvergenceCheck contains @@ -1214,8 +1213,7 @@ subroutine integrateStateEuler enddo; enddo; enddo !$OMP END PARALLEL DO - if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false. - if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck + if (nonlocalBroken) call nonlocalConvergenceCheck end subroutine integrateStateEuler @@ -1324,7 +1322,7 @@ subroutine integrateStateAdaptiveEuler enddo; enddo; enddo !$OMP END PARALLEL DO - if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck + if (nonlocalBroken) call nonlocalConvergenceCheck end subroutine integrateStateAdaptiveEuler @@ -1478,8 +1476,7 @@ subroutine integrateStateRK4 enddo; enddo; enddo !$OMP END PARALLEL DO - if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false. - if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck + if (nonlocalBroken) call nonlocalConvergenceCheck end subroutine integrateStateRK4 @@ -1652,8 +1649,7 @@ subroutine integrateStateRKCK45 enddo; enddo; enddo !$OMP END PARALLEL DO - if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false. - if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck + if (nonlocalBroken) call nonlocalConvergenceCheck end subroutine integrateStateRKCK45 @@ -1664,8 +1660,7 @@ end subroutine integrateStateRKCK45 !-------------------------------------------------------------------------------------------------- subroutine nonlocalConvergenceCheck - if (any(.not. crystallite_converged .and. .not. crystallite_localPlasticity)) & ! any non-local not yet converged (or broken)... - where( .not. crystallite_localPlasticity) crystallite_converged = .false. + where( .not. crystallite_localPlasticity) crystallite_converged = .false. end subroutine nonlocalConvergenceCheck From 01818cba80df0f80745cefcb0d394af9de56700b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 31 Mar 2020 11:05:25 +0200 Subject: [PATCH 25/35] tuples not needed for single arguments --- python/tests/test_Geom.py | 18 +++++++++--------- python/tests/test_VTK.py | 1 - 2 files changed, 9 insertions(+), 10 deletions(-) diff --git a/python/tests/test_Geom.py b/python/tests/test_Geom.py index f988179fa..0f8c10066 100644 --- a/python/tests/test_Geom.py +++ b/python/tests/test_Geom.py @@ -72,7 +72,7 @@ class TestGeom: if update: modified.to_file(reference) assert geom_equal(modified,Geom.from_file(reference)) - @pytest.mark.parametrize('stencil',[(1),(2),(3),(4)]) + @pytest.mark.parametrize('stencil',[1,2,3,4]) def test_clean(self,default,update,reference_dir,stencil): modified = copy.deepcopy(default) modified.clean(stencil) @@ -82,12 +82,12 @@ class TestGeom: assert geom_equal(modified,Geom.from_file(reference)) @pytest.mark.parametrize('grid',[ - ((10,11,10)), - ([10,13,10]), - (np.array((10,10,10))), - (np.array((8, 10,12))), - (np.array((5, 4, 20))), - (np.array((10,20,2)) ) + (10,11,10), + [10,13,10], + np.array((10,10,10)), + np.array((8, 10,12)), + np.array((5, 4, 20)), + np.array((10,20,2)) ] ) def test_scale(self,default,update,reference_dir,grid): @@ -98,7 +98,7 @@ class TestGeom: if update: modified.to_file(reference) assert geom_equal(modified,Geom.from_file(reference)) - @pytest.mark.parametrize('periodic',[(True),(False)]) + @pytest.mark.parametrize('periodic',[True,False]) def test_tessellation_approaches(self,periodic): grid = np.random.randint(10,20,3) size = np.random.random(3) + 1.0 @@ -119,7 +119,7 @@ class TestGeom: Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,weights,np.random.random()>0.5) assert np.all(Laguerre.microstructure == ms) - @pytest.mark.parametrize('approach',[('Laguerre'),('Voronoi')]) + @pytest.mark.parametrize('approach',['Laguerre','Voronoi']) def test_tessellate_bicrystal(self,approach): grid = np.random.randint(5,10,3)*2 size = grid.astype(np.float) diff --git a/python/tests/test_VTK.py b/python/tests/test_VTK.py index ef4090da5..627ea4007 100644 --- a/python/tests/test_VTK.py +++ b/python/tests/test_VTK.py @@ -1,5 +1,4 @@ import os -import sys import pytest import numpy as np From 6254063b4b1154f80fcbaf5b7ac54cc3fe8d6d3d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 31 Mar 2020 12:52:46 +0200 Subject: [PATCH 26/35] not used --- processing/pre/geom_fromVoronoiTessellation.py | 1 - 1 file changed, 1 deletion(-) diff --git a/processing/pre/geom_fromVoronoiTessellation.py b/processing/pre/geom_fromVoronoiTessellation.py index 31e8caa7b..df40176d8 100755 --- a/processing/pre/geom_fromVoronoiTessellation.py +++ b/processing/pre/geom_fromVoronoiTessellation.py @@ -192,7 +192,6 @@ for name in filenames: grains = table.get(options.microstructure) if options.microstructure in table.labels else np.arange(len(seeds))+1 grainIDs = np.unique(grains).astype('i') - NgrainIDs = len(grainIDs) if options.eulers in table.labels: eulers = table.get(options.eulers) From e818dfdb3e5c6d9a4bd44a05291b951ffa02cfee Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 31 Mar 2020 20:42:51 +0200 Subject: [PATCH 27/35] not used anymore --- src/prec.f90 | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/prec.f90 b/src/prec.f90 index a8d7d9e33..974c63c2f 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -30,10 +30,6 @@ module prec real(pReal), dimension(:), pointer :: p end type group_float - type :: group_int - integer, dimension(:), pointer :: p - end type group_int - ! http://stackoverflow.com/questions/3948210/can-i-have-a-pointer-to-an-item-in-an-allocatable-array type :: tState integer :: & From 9c95ce36f4b6ef1a783fce150481501f889891e0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 31 Mar 2020 20:57:09 +0200 Subject: [PATCH 28/35] automatic LHS (re)-allocation does not work for pointers group_float has pointers, not allocatables --- src/crystallite.f90 | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 5c53e9165..d547a6a01 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1011,7 +1011,7 @@ subroutine integrateStateFPI real(pReal), dimension(max(constitutive_plasticity_maxSizeDotState,constitutive_source_maxSizeDotState)) :: & r ! state residuum real(pReal), dimension(:), allocatable :: plastic_dotState_p1, plastic_dotState_p2 - type(group_float), dimension(maxval(phase_Nsources)) :: source_dotState_p1, source_dotState_p2 + real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState logical :: & nonlocalBroken @@ -1047,7 +1047,7 @@ subroutine integrateStateFPI 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) - source_dotState_p2(s)%p = 0.0_pReal * sourceState(p)%p(s)%dotState (1:sizeDotState,c) ! ToDo can be done smarter/clearer + source_dotState(1:sizeDotState,2,s) = 0.0_pReal enddo iteration: do NiterationState = 1, num%nState @@ -1055,8 +1055,9 @@ subroutine integrateStateFPI if(nIterationState > 1) plastic_dotState_p2 = plastic_dotState_p1 plastic_dotState_p1 = plasticState(p)%dotState(:,c) do s = 1, phase_Nsources(p) - if(nIterationState > 1) source_dotState_p2(s)%p = source_dotState_p1(s)%p - source_dotState_p1(s)%p = sourceState(p)%p(s)%dotState(:,c) + sizeDotState = sourceState(p)%p(s)%sizeDotState + if(nIterationState > 1) source_dotState(1:sizeDotState,2,s) = source_dotState(1:sizeDotState,1,s) + source_dotState(1:sizeDotState,1,s) = sourceState(p)%p(s)%dotState(:,c) enddo call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & @@ -1091,9 +1092,11 @@ subroutine integrateStateFPI plasticState(p)%atol(1:sizeDotState)) do s = 1, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - zeta = damper(sourceState(p)%p(s)%dotState(:,c),source_dotState_p1(s)%p,source_dotState_p2(s)%p) + zeta = damper(sourceState(p)%p(s)%dotState(:,c), & + source_dotState(1:sizeDotState,1,s),& + source_dotState(1:sizeDotState,2,s)) sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) * zeta & - + source_dotState_p1(s)%p* (1.0_pReal - zeta) + + source_dotState(1:sizeDotState,1,s)* (1.0_pReal - zeta) r(1: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) From e46220cd8a471039438d967fae0dd698e97a727e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 31 Mar 2020 22:02:53 +0200 Subject: [PATCH 29/35] OMP bugfix for FPI integrator, memory-efficient RK4 --- src/crystallite.f90 | 31 ++++++++++++++++++------------- src/material.f90 | 4 ---- src/prec.f90 | 5 +---- 3 files changed, 19 insertions(+), 21 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index d547a6a01..9617cc283 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1016,7 +1016,7 @@ subroutine integrateStateFPI nonlocalBroken nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE(sizeDotState,r,zeta,p,c,plastic_dotState_p1, plastic_dotState_p2,source_dotState_p1, source_dotState_p2) + !$OMP PARALLEL DO PRIVATE(sizeDotState,r,zeta,p,c,plastic_dotState_p1, plastic_dotState_p2,source_dotState) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) @@ -1359,8 +1359,10 @@ subroutine integrateStateRK4 logical :: & nonlocalBroken + real(pReal), dimension(constitutive_source_maxSizeDotState,4,maxval(phase_Nsources)) :: source_RK4dotState + real(pReal), dimension(constitutive_plasticity_maxSizeDotState,4) :: plastic_RK4dotState nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c) + !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,source_RK4dotState,plastic_RK4dotState) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) @@ -1382,20 +1384,23 @@ subroutine integrateStateRK4 if(.not. crystallite_todo(g,i,e)) cycle do stage = 1,3 - - plasticState(p)%RK4dotState(stage,:,c) = plasticState(p)%dotState(:,c) - plasticState(p)%dotState(:,c) = A(1,stage) * plasticState(p)%RK4dotState(1,:,c) + sizeDotState = plasticState(p)%sizeDotState + plastic_RK4dotState(1:sizeDotState,stage) = plasticState(p)%dotState(:,c) + plasticState(p)%dotState(:,c) = A(1,stage) * plastic_RK4dotState(1:sizeDotState,1) do s = 1, phase_Nsources(p) - sourceState(p)%p(s)%RK4dotState(stage,:,c) = sourceState(p)%p(s)%dotState(:,c) - sourceState(p)%p(s)%dotState(:,c) = A(1,stage) * sourceState(p)%p(s)%RK4dotState(1,:,c) + sizeDotState = sourceState(p)%p(s)%sizeDotState + source_RK4dotState(1:sizeDotState,stage,s) = sourceState(p)%p(s)%dotState(:,c) + sourceState(p)%p(s)%dotState(:,c) = A(1,stage) * source_RK4dotState(1:sizeDotState,1,s) enddo do n = 2, stage + sizeDotState = plasticState(p)%sizeDotState plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) & - + A(n,stage) * plasticState(p)%RK4dotState(n,:,c) + + A(n,stage) * plastic_RK4dotState(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) * sourceState(p)%p(s)%RK4dotState(n,:,c) + + A(n,stage) * source_RK4dotState(1:sizeDotState,n,s) enddo enddo @@ -1438,9 +1443,9 @@ subroutine integrateStateRK4 sizeDotState = plasticState(p)%sizeDotState - plasticState(p)%RK4dotState(4,:,c) = plasticState (p)%dotState(:,c) + plastic_RK4dotState(1:sizeDotState,4) = plasticState (p)%dotState(:,c) - plasticState(p)%dotState(:,c) = matmul(B,plasticState(p)%RK4dotState(1:4,1:sizeDotState,c)) + plasticState(p)%dotState(:,c) = matmul(plastic_RK4dotState(1:sizeDotState,1:4),B) plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & + plasticState(p)%dotState (1:sizeDotState,c) & * crystallite_subdt(g,i,e) @@ -1448,9 +1453,9 @@ subroutine integrateStateRK4 do s = 1, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - sourceState(p)%p(s)%RK4dotState(4,:,c) = sourceState(p)%p(s)%dotState(:,c) + source_RK4dotState(1:sizeDotState,4,s) = sourceState(p)%p(s)%dotState(:,c) - sourceState(p)%p(s)%dotState(:,c) = matmul(B,sourceState(p)%p(s)%RK4dotState(1:4,1:sizeDotState,c)) + sourceState(p)%p(s)%dotState(:,c) = matmul(source_RK4dotState(1:sizeDotState,1:4,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) diff --git a/src/material.f90 b/src/material.f90 index 14d4e4c19..fc55bb17c 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -724,8 +724,6 @@ subroutine material_allocatePlasticState(phase,NipcMyPhase,& allocate(plasticState(phase)%state (sizeState,NipcMyPhase), source=0.0_pReal) allocate(plasticState(phase)%dotState (sizeDotState,NipcMyPhase),source=0.0_pReal) - if (numerics_integrator == 4) & - allocate(plasticState(phase)%RK4dotState (4,sizeDotState,NipcMyPhase),source=0.0_pReal) if (numerics_integrator == 5) & allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NipcMyPhase),source=0.0_pReal) @@ -758,8 +756,6 @@ subroutine material_allocateSourceState(phase,of,NipcMyPhase,& allocate(sourceState(phase)%p(of)%state (sizeState,NipcMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(of)%dotState (sizeDotState,NipcMyPhase),source=0.0_pReal) - if (numerics_integrator == 4) & - allocate(sourceState(phase)%p(of)%RK4dotState (4,sizeDotState,NipcMyPhase),source=0.0_pReal) if (numerics_integrator == 5) & allocate(sourceState(phase)%p(of)%RKCK45dotState (6,sizeDotState,NipcMyPhase),source=0.0_pReal) diff --git a/src/prec.f90 b/src/prec.f90 index 974c63c2f..0afaeb536 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -46,11 +46,8 @@ module prec deltaState !< increment of state change real(pReal), allocatable, dimension(:,:) :: & partionedState0, & - subState0, & - previousDotState, & - previousDotState2 + subState0 real(pReal), allocatable, dimension(:,:,:) :: & - RK4dotState, & RKCK45dotState end type From 6bfa51c30777d6a9ccab1898af32a27ec32f6c20 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 07:02:23 +0200 Subject: [PATCH 30/35] LHS allocation does not work for pointers --- src/crystallite.f90 | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 9617cc283..5884e2cef 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1237,11 +1237,11 @@ subroutine integrateStateAdaptiveEuler logical :: & nonlocalBroken - real(pReal), dimension(:), allocatable :: residuum_plastic - type(group_float), dimension(maxval(phase_Nsources)) :: residuum_source + 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) + !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,residuum_plastic,residuum_source) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) @@ -1264,14 +1264,14 @@ subroutine integrateStateAdaptiveEuler sizeDotState = plasticState(p)%sizeDotState - residuum_plastic = - plasticState(p)%dotstate(1:sizeDotState,c) * 0.5_pReal * crystallite_subdt(g,i,e) + 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(s)%p = - sourceState(p)%p(s)%dotstate(1:sizeDotState,c) & - * 0.5_pReal * crystallite_subdt(g,i,e) + 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 @@ -1306,7 +1306,7 @@ subroutine integrateStateAdaptiveEuler sizeDotState = plasticState(p)%sizeDotState - crystallite_converged(g,i,e) = converged(residuum_plastic & + 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)) @@ -1315,7 +1315,7 @@ subroutine integrateStateAdaptiveEuler sizeDotState = sourceState(p)%p(s)%sizeDotState crystallite_converged(g,i,e) = & - crystallite_converged(g,i,e) .and. converged(residuum_source(s)%p & + 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)) From 54e3455bd47a31a9e82b63b7b7e6be2bcb04e313 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 07:09:33 +0200 Subject: [PATCH 31/35] RKCK45 store intermediate state only per point --- src/crystallite.f90 | 33 +++++++++++++++++++-------------- src/material.f90 | 4 ---- src/prec.f90 | 2 -- 3 files changed, 19 insertions(+), 20 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 5884e2cef..fb12c439b 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1526,9 +1526,11 @@ subroutine integrateStateRKCK45 sizeDotState logical :: & nonlocalBroken + real(pReal), dimension(constitutive_source_maxSizeDotState,6,maxval(phase_Nsources)) :: source_RKdotState + real(pReal), dimension(constitutive_plasticity_maxSizeDotState,6) :: plastic_RKdotState nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c) + !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,plastic_RKdotState,source_RKdotState) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) @@ -1550,20 +1552,23 @@ subroutine integrateStateRKCK45 if(.not. crystallite_todo(g,i,e)) cycle do stage = 1,5 - - plasticState(p)%RKCK45dotState(stage,:,c) = plasticState(p)%dotState(:,c) - plasticState(p)%dotState(:,c) = A(1,stage) * plasticState(p)%RKCK45dotState(1,:,c) + 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) - sourceState(p)%p(s)%RKCK45dotState(stage,:,c) = sourceState(p)%p(s)%dotState(:,c) - sourceState(p)%p(s)%dotState(:,c) = A(1,stage) * sourceState(p)%p(s)%RKCK45dotState(1,:,c) + 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) * plasticState(p)%RKCK45dotState(n,:,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) * sourceState(p)%p(s)%RKCK45dotState(n,:,c) + + A(n,stage) * source_RKdotState(1:sizeDotState,n,s) enddo enddo @@ -1606,12 +1611,12 @@ subroutine integrateStateRKCK45 sizeDotState = plasticState(p)%sizeDotState - plasticState(p)%RKCK45dotState(6,:,c) = plasticState (p)%dotState(:,c) - plasticState(p)%dotState(:,c) = matmul(B,plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,c)) + plastic_RKdotState(1:sizeDotState,6) = plasticState (p)%dotState(:,c) + plasticState(p)%dotState(:,c) = matmul(plastic_RKdotState(1:sizeDotState,1:6),B) plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & + plasticState(p)%dotState (1:sizeDotState,c) & * crystallite_subdt(g,i,e) - crystallite_todo(g,i,e) = converged(matmul(DB,plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,c)) & + crystallite_todo(g,i,e) = converged( matmul(plastic_RKdotState(1:sizeDotState,1:6),DB) & * crystallite_subdt(g,i,e), & plasticState(p)%state(1:sizeDotState,c), & plasticState(p)%atol(1:sizeDotState)) @@ -1619,13 +1624,13 @@ subroutine integrateStateRKCK45 do s = 1, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - sourceState(p)%p(s)%RKCK45dotState(6,:,c) = sourceState(p)%p(s)%dotState(:,c) - sourceState(p)%p(s)%dotState(:,c) = matmul(B,sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c)) + source_RKdotState(1:sizeDotState,6,s) = sourceState(p)%p(s)%dotState(:,c) + sourceState(p)%p(s)%dotState(:,c) = matmul(source_RKdotState(1:sizeDotState,1:6,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) crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. & - converged(matmul(DB,sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c)) & + converged(matmul(source_RKdotState(1:sizeDotState,1:6,s),DB) & * crystallite_subdt(g,i,e), & sourceState(p)%p(s)%state(1:sizeDotState,c), & sourceState(p)%p(s)%atol(1:sizeDotState)) diff --git a/src/material.f90 b/src/material.f90 index fc55bb17c..aefe44878 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -724,8 +724,6 @@ subroutine material_allocatePlasticState(phase,NipcMyPhase,& allocate(plasticState(phase)%state (sizeState,NipcMyPhase), source=0.0_pReal) allocate(plasticState(phase)%dotState (sizeDotState,NipcMyPhase),source=0.0_pReal) - if (numerics_integrator == 5) & - allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NipcMyPhase),source=0.0_pReal) allocate(plasticState(phase)%deltaState (sizeDeltaState,NipcMyPhase),source=0.0_pReal) @@ -756,8 +754,6 @@ subroutine material_allocateSourceState(phase,of,NipcMyPhase,& allocate(sourceState(phase)%p(of)%state (sizeState,NipcMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(of)%dotState (sizeDotState,NipcMyPhase),source=0.0_pReal) - if (numerics_integrator == 5) & - allocate(sourceState(phase)%p(of)%RKCK45dotState (6,sizeDotState,NipcMyPhase),source=0.0_pReal) allocate(sourceState(phase)%p(of)%deltaState (sizeDeltaState,NipcMyPhase),source=0.0_pReal) diff --git a/src/prec.f90 b/src/prec.f90 index 0afaeb536..dd39ddc05 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -47,8 +47,6 @@ module prec real(pReal), allocatable, dimension(:,:) :: & partionedState0, & subState0 - real(pReal), allocatable, dimension(:,:,:) :: & - RKCK45dotState end type type, extends(tState) :: tPlasticState From a25e163ab97c38506c846b5cc9aa0ff8b702f478 Mon Sep 17 00:00:00 2001 From: Test User Date: Thu, 2 Apr 2020 14:37:53 +0200 Subject: [PATCH 32/35] [skip ci] updated version information after successful test of v2.0.3-2206-gf174dd6a --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index d111f6358..43d4d275b 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-2204-gdb1f2151 +v2.0.3-2206-gf174dd6a From 5e9ff7947bfbd8ef1d122cc412b517fef7df9b75 Mon Sep 17 00:00:00 2001 From: Franz Roters Date: Fri, 3 Apr 2020 14:01:35 +0200 Subject: [PATCH 33/35] [skip ci] plastic_dotstate always before source_dotstate --- src/crystallite.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index fb12c439b..efa066696 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1359,10 +1359,10 @@ subroutine integrateStateRK4 logical :: & nonlocalBroken - real(pReal), dimension(constitutive_source_maxSizeDotState,4,maxval(phase_Nsources)) :: source_RK4dotState real(pReal), dimension(constitutive_plasticity_maxSizeDotState,4) :: plastic_RK4dotState + real(pReal), dimension(constitutive_source_maxSizeDotState,4,maxval(phase_Nsources)) :: source_RK4dotState nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,source_RK4dotState,plastic_RK4dotState) + !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,plastic_RK4dotState,source_RK4dotState) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) @@ -1526,8 +1526,8 @@ subroutine integrateStateRKCK45 sizeDotState logical :: & nonlocalBroken - real(pReal), dimension(constitutive_source_maxSizeDotState,6,maxval(phase_Nsources)) :: source_RKdotState real(pReal), dimension(constitutive_plasticity_maxSizeDotState,6) :: plastic_RKdotState + real(pReal), dimension(constitutive_source_maxSizeDotState,6,maxval(phase_Nsources)) :: source_RKdotState nonlocalBroken = .false. !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,plastic_RKdotState,source_RKdotState) From 407fb12f8c269050e7a7e33b2c957da9221175ee Mon Sep 17 00:00:00 2001 From: Test User Date: Fri, 3 Apr 2020 18:28:16 +0200 Subject: [PATCH 34/35] [skip ci] updated version information after successful test of v2.0.3-2223-g8631653f --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 43d4d275b..d2968346e 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-2206-gf174dd6a +v2.0.3-2223-g8631653f From a86f00d82752dc0bf1139ab1428e3c9ab0bb87c1 Mon Sep 17 00:00:00 2001 From: Test User Date: Sun, 5 Apr 2020 02:20:26 +0200 Subject: [PATCH 35/35] [skip ci] updated version information after successful test of v2.0.3-2243-gbb03483b --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index d2968346e..2d8aa9597 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-2223-g8631653f +v2.0.3-2243-gbb03483b