diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index ce645dcf5..a56b108d5 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -104,7 +104,7 @@ checkout: - release ################################################################################################### -Pytest: +Pytest_python: stage: python script: - cd $DAMASKROOT/python @@ -385,6 +385,15 @@ Phenopowerlaw_singleSlip: - master - release +Pytest_grid: + stage: grid + script: + - cd pytest + - pytest + except: + - master + - release + ################################################################################################### Marc_compileIfort: stage: compileMarc diff --git a/PRIVATE b/PRIVATE index 6db5f4666..9573ce7bd 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 6db5f4666fc651b4de3b44ceaed3f2b848170ac9 +Subproject commit 9573ce7bd2c1a7188c1aac5b83aa76d480c2bdb0 diff --git a/examples/ConfigFiles/Phase_Phenopowerlaw_multiField.config b/examples/ConfigFiles/Phase_Phenopowerlaw_multiField.config deleted file mode 100644 index adad3710e..000000000 --- a/examples/ConfigFiles/Phase_Phenopowerlaw_multiField.config +++ /dev/null @@ -1,49 +0,0 @@ -[Aluminum] -elasticity hooke -plasticity phenopowerlaw - -(output) resistance_slip -(output) shearrate_slip -(output) resolvedstress_slip -(output) accumulated_shear_slip -(output) resistance_twin -(output) shearrate_twin -(output) resolvedstress_twin -(output) accumulated_shear_twin - -lattice_structure fcc -Nslip 12 # per family -Ntwin 0 # per family - -c11 106.75e9 -c12 60.41e9 -c44 28.34e9 - -gdot0_slip 0.001 -n_slip 20 -tau0_slip 31e6 # per family -tausat_slip 63e6 # per family -a_slip 2.25 -gdot0_twin 0.001 -n_twin 20 -tau0_twin 31e6 # per family -h0_slipslip 75e6 -interaction_slipslip 1 1 1.4 1.4 1.4 1.4 -atol_resistance 1 - -(stiffness_degradation) damage -(stiffness_degradation) porosity -{./Phase_Damage.config} -{./Phase_Thermal.config} -{./Phase_Vacancy.config} -{./Phase_Porosity.config} -{./Phase_Hydrogen.config} -{./Source_Damage_IsoBrittle.config} -{./Source_Thermal_Dissipation.config} -{./Source_Vacancy_PhenoPlasticity.config} -{./Source_Vacancy_Irradiation.config} -{./Kinematics_Thermal_Expansion.config} -{./Kinematics_Vacancy_Strain.config} -{./Kinematics_Hydrogen_Strain.config} - - diff --git a/processing/post/permuteData.py b/processing/post/permuteData.py index a589cb6da..81af71adb 100755 --- a/processing/post/permuteData.py +++ b/processing/post/permuteData.py @@ -83,7 +83,7 @@ for name in filenames: # ------------------------------------------ assemble header --------------------------------------- - randomSeed = int(os.urandom(4).encode('hex'), 16) if options.randomSeed is None else options.randomSeed # random seed per file + randomSeed = int(os.urandom(4).hex(), 16) if options.randomSeed is None else options.randomSeed # random seed per file np.random.seed(randomSeed) table.info_append([scriptID + '\t' + ' '.join(sys.argv[1:]), diff --git a/processing/pre/seeds_fromDistribution.py b/processing/pre/seeds_fromDistribution.py index 2e8936f27..7bd97db60 100755 --- a/processing/pre/seeds_fromDistribution.py +++ b/processing/pre/seeds_fromDistribution.py @@ -1,11 +1,9 @@ #!/usr/bin/env python3 -# -*- coding: UTF-8 no BOM -*- import threading,time,os,sys,random import numpy as np from optparse import OptionParser from io import StringIO -import binascii import damask scriptName = os.path.splitext(os.path.basename(__file__))[0] @@ -16,9 +14,10 @@ currentSeedsName = None #--------------------------------------------------------------------------------------------------- class myThread (threading.Thread): - """perturbes seed in seed file, performes Voronoi tessellation, evaluates, and updates best match""" + """Perturb seed in seed file, performes Voronoi tessellation, evaluates, and updates best match.""" def __init__(self, threadID): + """Threading class with thread ID.""" threading.Thread.__init__(self) self.threadID = threadID @@ -215,7 +214,7 @@ options = parser.parse_args()[0] damask.util.report(scriptName,options.seedFile) if options.randomSeed is None: - options.randomSeed = int(binascii.hexlify(os.urandom(4)),16) + options.randomSeed = int(os.urandom(4).hex(),16) damask.util.croak(options.randomSeed) delta = (options.scale/options.grid[0],options.scale/options.grid[1],options.scale/options.grid[2]) baseFile=os.path.splitext(os.path.basename(options.seedFile))[0] diff --git a/src/IO.f90 b/src/IO.f90 index c31ebdfab..6ab87715c 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -400,6 +400,8 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) msg = 'number of values does not match' case (147) msg = 'not supported anymore' + case (148) + msg = 'Nconstituents mismatch between homogenization and microstructure' !-------------------------------------------------------------------------------------------------- ! material error messages and related messages in mesh diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 0678f9631..1a8252cce 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -729,7 +729,7 @@ subroutine crystallite_results real(pReal), allocatable, dimension(:,:,:) :: select_tensors integer :: e,i,c,j - allocate(select_tensors(3,3,count(material_phaseAt==instance)*homogenization_maxNgrains*discretization_nIP)) + allocate(select_tensors(3,3,count(material_phaseAt==instance)*discretization_nIP)) j=0 do e = 1, size(material_phaseAt,2) diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index b51ebb56a..57a27ebe8 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -21,10 +21,10 @@ module grid_mech_FEM use discretization use mesh_grid use debug - + implicit none private - + !-------------------------------------------------------------------------------------------------- ! derived types type(tSolutionParams), private :: params @@ -52,20 +52,20 @@ module grid_mech_FEM 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=pStringLen), private :: incInfo !< time and increment information - + real(pReal), private, dimension(3,3,3,3) :: & - C_volAvg = 0.0_pReal, & !< current volume average stiffness + 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 :: & err_BC !< deviation from stress BC - + integer, private :: & totalIter = 0 !< total iteration in current increment - + public :: & grid_mech_FEM_init, & grid_mech_FEM_solution, & @@ -79,7 +79,7 @@ contains !> @brief allocates all necessary fields and fills them with data, potentially from restart info !-------------------------------------------------------------------------------------------------- subroutine grid_mech_FEM_init - + real(pReal) :: HGCoeff = 0.0e-2_pReal PetscInt, dimension(0:worldsize-1) :: localK real(pReal), dimension(3,3) :: & @@ -99,9 +99,9 @@ subroutine grid_mech_FEM_init real(pReal), dimension(3,3,3,3) :: devNull PetscScalar, pointer, dimension(:,:,:,:) :: & u_current,u_lastInc - + write(6,'(/,a)') ' <<<+- grid_mech_FEM init -+>>>'; flush(6) - + !-------------------------------------------------------------------------------------------------- ! set default and user defined options for PETSc call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type newtonls -mech_ksp_type fgmres & @@ -115,11 +115,11 @@ subroutine grid_mech_FEM_init allocate(F (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) allocate(P_current (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) allocate(F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - + !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc call SNESCreate(PETSC_COMM_WORLD,mech_snes,ierr); CHKERRQ(ierr) - call SNESSetOptionsPrefix(mech_snes,'mech_',ierr);CHKERRQ(ierr) + call SNESSetOptionsPrefix(mech_snes,'mech_',ierr);CHKERRQ(ierr) localK = 0 localK(worldrank) = grid3 call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) @@ -141,12 +141,12 @@ subroutine grid_mech_FEM_init call DMCreateGlobalVector(mech_grid,solution_lastInc,ierr); CHKERRQ(ierr) call DMCreateGlobalVector(mech_grid,solution_rate ,ierr); CHKERRQ(ierr) call DMSNESSetFunctionLocal(mech_grid,formResidual,PETSC_NULL_SNES,ierr) - CHKERRQ(ierr) + CHKERRQ(ierr) call DMSNESSetJacobianLocal(mech_grid,formJacobian,PETSC_NULL_SNES,ierr) CHKERRQ(ierr) call SNESSetConvergenceTest(mech_snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged" CHKERRQ(ierr) - call SNESSetMaxLinearSolveFailures(mech_snes, huge(1), ierr); CHKERRQ(ierr) ! ignore linear solve failures + call SNESSetMaxLinearSolveFailures(mech_snes, huge(1), ierr); CHKERRQ(ierr) ! ignore linear solve failures call SNESSetFromOptions(mech_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional cli arguments !-------------------------------------------------------------------------------------------------- @@ -156,15 +156,15 @@ subroutine grid_mech_FEM_init call VecSet(solution_rate ,0.0_pReal,ierr);CHKERRQ(ierr) call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) - + call DMDAGetCorners(mech_grid,xstart,ystart,zstart,xend,yend,zend,ierr) ! local grid extent - CHKERRQ(ierr) + CHKERRQ(ierr) xend = xstart+xend-1 yend = ystart+yend-1 zend = zstart+zend-1 delta = geomSize/real(grid,pReal) ! grid spacing detJ = product(delta) ! cell volume - + BMat = reshape(real([-1.0_pReal/delta(1),-1.0_pReal/delta(2),-1.0_pReal/delta(3), & 1.0_pReal/delta(1),-1.0_pReal/delta(2),-1.0_pReal/delta(3), & -1.0_pReal/delta(1), 1.0_pReal/delta(2),-1.0_pReal/delta(3), & @@ -173,7 +173,7 @@ subroutine grid_mech_FEM_init 1.0_pReal/delta(1),-1.0_pReal/delta(2), 1.0_pReal/delta(3), & -1.0_pReal/delta(1), 1.0_pReal/delta(2), 1.0_pReal/delta(3), & 1.0_pReal/delta(1), 1.0_pReal/delta(2), 1.0_pReal/delta(3)],pReal), [3,8])/4.0_pReal ! shape function derivative matrix - + HGMat = matmul(transpose(HGcomp),HGcomp) & * HGCoeff*(delta(1)*delta(2) + delta(2)*delta(3) + delta(3)*delta(1))/16.0_pReal ! hourglass stabilization matrix @@ -181,11 +181,11 @@ subroutine grid_mech_FEM_init ! init fields restartRead: if (interface_restartInc > 0) then write(6,'(/,a,i0,a)') ' reading restart data of increment ', interface_restartInc, ' from file' - + write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5' fileHandle = HDF5_openFile(fileName) groupHandle = HDF5_openGroup(fileHandle,'solver') - + call HDF5_read(groupHandle,F_aim, 'F_aim') call HDF5_read(groupHandle,F_aim_lastInc,'F_aim_lastInc') call HDF5_read(groupHandle,F_aimDot, 'F_aimDot') @@ -193,7 +193,7 @@ subroutine grid_mech_FEM_init call HDF5_read(groupHandle,F_lastInc, 'F_lastInc') call HDF5_read(groupHandle,u_current, 'u') call HDF5_read(groupHandle,u_lastInc, 'u_lastInc') - + elseif (interface_restartInc == 0) then restartRead 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) @@ -207,12 +207,12 @@ subroutine grid_mech_FEM_init CHKERRQ(ierr) call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr) CHKERRQ(ierr) - + restartRead2: if (interface_restartInc > 0) then write(6,'(/,a,i0,a)') ' reading more restart data of increment ', interface_restartInc, ' from file' call HDF5_read(groupHandle,C_volAvg, 'C_volAvg') call HDF5_read(groupHandle,C_volAvgLastInc,'C_volAvgLastInc') - + call HDF5_closeGroup(groupHandle) call HDF5_closeFile(fileHandle) @@ -243,14 +243,14 @@ function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation ! PETSc Data PetscErrorCode :: ierr SNESConvergedReason :: reason - + incInfo = incInfoIn !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) !-------------------------------------------------------------------------------------------------- -! set module wide available data +! set module wide available data params%stress_mask = stress_BC%maskFloat params%stress_BC = stress_BC%values params%rotation_BC = rotation_BC @@ -258,13 +258,13 @@ function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation params%timeincOld = timeinc_old !-------------------------------------------------------------------------------------------------- -! solve BVP +! solve BVP call SNESsolve(mech_snes,PETSC_NULL_VEC,solution_current,ierr);CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! check convergence call SNESGetConvergedReason(mech_snes,reason,ierr);CHKERRQ(ierr) - + solution%converged = reason > 0 solution%iterationsNeeded = totalIter solution%termIll = terminallyIll @@ -296,15 +296,15 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime, PetscErrorCode :: ierr PetscScalar, pointer, dimension(:,:,:,:) :: & u_current,u_lastInc - + call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) - + if (cutBack) then C_volAvg = C_volAvgLastInc else C_volAvgLastInc = C_volAvg - + F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess) F_aim_lastInc = F_aim @@ -329,10 +329,10 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime, call VecSet(solution_rate,0.0_pReal,ierr); CHKERRQ(ierr) endif call VecCopy(solution_current,solution_lastInc,ierr); CHKERRQ(ierr) - + F_lastInc = F - - materialpoint_F0 = reshape(F, [3,3,1,product(grid(1:2))*grid3]) + + materialpoint_F0 = reshape(F, [3,3,1,product(grid(1:2))*grid3]) endif !-------------------------------------------------------------------------------------------------- @@ -365,12 +365,12 @@ subroutine grid_mech_FEM_restartWrite integer(HID_T) :: fileHandle, groupHandle PetscScalar, dimension(:,:,:,:), pointer :: u_current,u_lastInc character(len=pStringLen) :: fileName - + call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) write(6,'(a)') ' writing solver data required for restart to file'; flush(6) - + write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5' fileHandle = HDF5_openFile(fileName,'w') groupHandle = HDF5_addGroup(fileHandle,'solver') @@ -388,7 +388,7 @@ subroutine grid_mech_FEM_restartWrite call HDF5_closeGroup(groupHandle) call HDF5_closeFile(fileHandle) - + call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr);CHKERRQ(ierr) call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr);CHKERRQ(ierr) @@ -405,7 +405,7 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,i PetscReal, intent(in) :: & devNull1, & devNull2, & - fnorm + fnorm SNESConvergedReason :: reason PetscObject :: dummy PetscErrorCode :: ierr @@ -421,7 +421,7 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,i if ((totalIter >= itmin .and. & all([ err_div/divTol, & err_BC /BCTol ] < 1.0_pReal)) & - .or. terminallyIll) then + .or. terminallyIll) then reason = 1 elseif (totalIter >= itmax) then reason = -1 @@ -435,10 +435,10 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,i write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & err_div/divTol, ' (',err_div,' / m, tol = ',divTol,')' write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', & - err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' + err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' write(6,'(/,a)') ' ===========================================================================' flush(6) - + end subroutine converged @@ -475,7 +475,7 @@ subroutine formResidual(da_local,x_local, & write(6,'(1x,a,3(a,i0))') trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter+1, '≤', itmax if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotTensor2(F_aim,active=.true.)) + ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & ' deformation gradient aim =', transpose(F_aim) flush(6) @@ -491,7 +491,7 @@ subroutine formResidual(da_local,x_local, & x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk) enddo; enddo; enddo ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1 - F(1:3,1:3,ii,jj,kk) = params%rotation_BC%rotTensor2(F_aim,active=.true.) + transpose(matmul(BMat,x_elem)) + F(1:3,1:3,ii,jj,kk) = params%rotation_BC%rotate(F_aim,active=.true.) + transpose(matmul(BMat,x_elem)) enddo; enddo; enddo call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr) @@ -501,7 +501,7 @@ subroutine formResidual(da_local,x_local, & P_av,C_volAvg,devNull, & F,params%timeinc,params%rotation_BC) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) - + !-------------------------------------------------------------------------------------------------- ! stress BC handling F_aim_lastIter = F_aim @@ -535,24 +535,24 @@ subroutine formResidual(da_local,x_local, & enddo; enddo; enddo call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr) call DMDAVecRestoreArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr) - + !-------------------------------------------------------------------------------------------------- ! applying boundary conditions call DMDAVecGetArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr) if (zstart == 0) then - f_scal(0:2,xstart,ystart,zstart) = 0.0 - f_scal(0:2,xend+1,ystart,zstart) = 0.0 - f_scal(0:2,xstart,yend+1,zstart) = 0.0 - f_scal(0:2,xend+1,yend+1,zstart) = 0.0 + f_scal(0:2,xstart,ystart,zstart) = 0.0 + f_scal(0:2,xend+1,ystart,zstart) = 0.0 + f_scal(0:2,xstart,yend+1,zstart) = 0.0 + f_scal(0:2,xend+1,yend+1,zstart) = 0.0 endif if (zend + 1 == grid(3)) then - f_scal(0:2,xstart,ystart,zend+1) = 0.0 - f_scal(0:2,xend+1,ystart,zend+1) = 0.0 - f_scal(0:2,xstart,yend+1,zend+1) = 0.0 - f_scal(0:2,xend+1,yend+1,zend+1) = 0.0 + f_scal(0:2,xstart,ystart,zend+1) = 0.0 + f_scal(0:2,xend+1,ystart,zend+1) = 0.0 + f_scal(0:2,xstart,yend+1,zend+1) = 0.0 + f_scal(0:2,xend+1,yend+1,zend+1) = 0.0 endif - call DMDAVecRestoreArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr) - + call DMDAVecRestoreArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr) + end subroutine formResidual @@ -574,7 +574,7 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,ierr) PetscObject :: dummy MatNullSpace :: matnull PetscErrorCode :: ierr - + BMatFull = 0.0 BMatFull(1:3,1 :8 ) = BMat BMatFull(4:6,9 :16) = BMat @@ -623,7 +623,7 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,ierr) call MatAssemblyEnd(Jac,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr) call MatAssemblyBegin(Jac_pre,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr) call MatAssemblyEnd(Jac_pre,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr) - + !-------------------------------------------------------------------------------------------------- ! applying boundary conditions diag = (C_volAvg(1,1,1,1)/delta(1)**2.0_pReal + & diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index af8cbc377..25ac9cb63 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -28,11 +28,11 @@ module grid_mech_spectral_basic !-------------------------------------------------------------------------------------------------- ! derived types type(tSolutionParams), private :: params - + type, private :: tNumerics logical :: update_gamma !< update gamma operator with current stiffness end type tNumerics - + type(tNumerics) :: num ! numerics parameters. Better name? !-------------------------------------------------------------------------------------------------- @@ -55,21 +55,21 @@ module grid_mech_spectral_basic F_aim_lastInc = math_I3, & !< previous average deformation gradient P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress - character(len=pStringLen), private :: incInfo !< time and increment information + character(len=pStringLen), private :: incInfo !< time and increment information real(pReal), private, dimension(3,3,3,3) :: & - C_volAvg = 0.0_pReal, & !< current volume average stiffness + 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 C_minMaxAvgLastInc = 0.0_pReal, & !< previous (min+max)/2 stiffness S = 0.0_pReal !< current compliance (filled up with zeros) - + real(pReal), private :: & err_BC, & !< deviation from stress BC err_div !< RMS of div of P - + integer, private :: & totalIter = 0 !< total iteration in current increment - + public :: & grid_mech_spectral_basic_init, & grid_mech_spectral_basic_solution, & @@ -83,27 +83,27 @@ contains !> @brief allocates all necessary fields and fills them with data, potentially from restart info !-------------------------------------------------------------------------------------------------- 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 - + PetscErrorCode :: ierr PetscScalar, pointer, dimension(:,:,:,:) :: & F ! pointer to solution data - PetscInt, dimension(worldsize) :: localK + PetscInt, dimension(worldsize) :: localK integer(HID_T) :: fileHandle, groupHandle integer :: fileUnit character(len=pStringLen) :: fileName - + write(6,'(/,a)') ' <<<+- grid_mech_spectral_basic init -+>>>'; flush(6) - + write(6,'(/,a)') ' Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013' write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2012.09.012' - + write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006' - + num%update_gamma = config_numerics%getInt('update_gamma',defaultVal=0) > 0 !-------------------------------------------------------------------------------------------------- @@ -117,11 +117,11 @@ subroutine grid_mech_spectral_basic_init ! allocate global fields allocate (F_lastInc(3,3,grid(1),grid(2),grid3),source = 0.0_pReal) allocate (Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - + !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) - call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) + call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) localK = 0 localK(worldrank+1) = grid3 call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) @@ -139,45 +139,45 @@ subroutine grid_mech_spectral_basic_init call DMsetUp(da,ierr); CHKERRQ(ierr) call DMcreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 9, i.e. every def grad tensor) call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector - CHKERRQ(ierr) + CHKERRQ(ierr) call SNESsetConvergenceTest(snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "converged" CHKERRQ(ierr) call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments !-------------------------------------------------------------------------------------------------- -! init fields +! init fields call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! places pointer on PETSc data - + restartRead: if (interface_restartInc > 0) then write(6,'(/,a,i0,a)') ' reading restart data of increment ', interface_restartInc, ' from file' write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5' fileHandle = HDF5_openFile(fileName) groupHandle = HDF5_openGroup(fileHandle,'solver') - + call HDF5_read(groupHandle,F_aim, 'F_aim') call HDF5_read(groupHandle,F_aim_lastInc,'F_aim_lastInc') call HDF5_read(groupHandle,F_aimDot, 'F_aimDot') call HDF5_read(groupHandle,F, 'F') call HDF5_read(groupHandle,F_lastInc, 'F_lastInc') - + elseif (interface_restartInc == 0) then restartRead F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity F = reshape(F_lastInc,[9,grid(1),grid(2),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(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 - + restartRead2: if (interface_restartInc > 0) then write(6,'(/,a,i0,a)') ' reading more restart data of increment ', interface_restartInc, ' from file' call HDF5_read(groupHandle,C_volAvg, 'C_volAvg') call HDF5_read(groupHandle,C_volAvgLastInc,'C_volAvgLastInc') - + call HDF5_closeGroup(groupHandle) call HDF5_closeFile(fileHandle) @@ -197,7 +197,7 @@ 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) - + !-------------------------------------------------------------------------------------------------- ! input data for solution character(len=*), intent(in) :: & @@ -215,7 +215,7 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_ ! PETSc Data PetscErrorCode :: ierr SNESConvergedReason :: reason - + incInfo = incInfoIn !-------------------------------------------------------------------------------------------------- @@ -224,7 +224,7 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_ if(num%update_gamma) call utilities_updateGamma(C_minMaxAvg) !-------------------------------------------------------------------------------------------------- -! set module wide available data +! set module wide available data params%stress_mask = stress_BC%maskFloat params%stress_BC = stress_BC%values params%rotation_BC = rotation_BC @@ -232,13 +232,13 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_ params%timeincOld = timeinc_old !-------------------------------------------------------------------------------------------------- -! solve BVP +! solve BVP call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! check convergence call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr) - + solution%converged = reason > 0 solution%iterationsNeeded = totalIter solution%termIll = terminallyIll @@ -271,14 +271,14 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo PetscScalar, dimension(:,:,:,:), pointer :: F call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) - + if (cutBack) then C_volAvg = C_volAvgLastInc C_minMaxAvg = C_minMaxAvgLastInc else 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_aim_lastInc = F_aim @@ -297,9 +297,9 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo Fdot = utilities_calculateRate(guess, & F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, & - rotation_BC%rotTensor2(F_aimDot,active=.true.)) + rotation_BC%rotate(F_aimDot,active=.true.)) F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid3]) - + materialpoint_F0 = reshape(F, [3,3,1,product(grid(1:2))*grid3]) endif @@ -307,9 +307,9 @@ 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 - rotation_BC%rotTensor2(F_aim,active=.true.)),[9,grid(1),grid(2),grid3]) + rotation_BC%rotate(F_aim,active=.true.)),[9,grid(1),grid(2),grid3]) call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) - + end subroutine grid_mech_spectral_basic_forward @@ -341,11 +341,11 @@ subroutine grid_mech_spectral_basic_restartWrite call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) write(6,'(a)') ' writing solver data required for restart to file'; flush(6) - + write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5' fileHandle = HDF5_openFile(fileName,'w') groupHandle = HDF5_addGroup(fileHandle,'solver') - + call HDF5_write(groupHandle,F_aim, 'F_aim') call HDF5_write(groupHandle,F_aim_lastInc,'F_aim_lastInc') call HDF5_write(groupHandle,F_aimDot, 'F_aimDot') @@ -358,7 +358,7 @@ subroutine grid_mech_spectral_basic_restartWrite call HDF5_closeGroup(groupHandle) call HDF5_closeFile(fileHandle) - + if (num%update_gamma) call utilities_saveReferenceStiffness call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) @@ -376,7 +376,7 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm PetscReal, intent(in) :: & devNull1, & devNull2, & - devNull3 + devNull3 SNESConvergedReason :: reason PetscObject :: dummy PetscErrorCode :: ierr @@ -390,7 +390,7 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm if ((totalIter >= itmin .and. & all([ err_div/divTol, & err_BC /BCTol ] < 1.0_pReal)) & - .or. terminallyIll) then + .or. terminallyIll) then reason = 1 elseif (totalIter >= itmax) then reason = -1 @@ -404,10 +404,10 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & err_div/divTol, ' (',err_div,' / m, tol = ',divTol,')' write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', & - err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' + err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' write(6,'(/,a)') ' ===========================================================================' - flush(6) - + flush(6) + end subroutine converged @@ -441,7 +441,7 @@ subroutine formResidual(in, F, & write(6,'(1x,a,3(a,i0))') trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotTensor2(F_aim,active=.true.)) + ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & ' deformation gradient aim =', transpose(F_aim) flush(6) @@ -453,7 +453,7 @@ subroutine formResidual(in, F, & P_av,C_volAvg,C_minMaxAvg, & F,params%timeinc,params%rotation_BC) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) - + !-------------------------------------------------------------------------------------------------- ! stress BC handling deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC) @@ -466,9 +466,9 @@ subroutine formResidual(in, F, & 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 - call utilities_fourierGammaConvolution(params%rotation_BC%rotTensor2(deltaF_aim,active=.true.)) ! convolution of Gamma and tensorField_fourier + 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 - + !-------------------------------------------------------------------------------------------------- ! constructing residual residuum = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index bdc65a8c5..ab880e2a6 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -22,20 +22,20 @@ module grid_mech_spectral_polarisation use homogenization use mesh_grid use debug - + implicit none private - + !-------------------------------------------------------------------------------------------------- ! derived types type(tSolutionParams), private :: params - + type, private :: tNumerics logical :: update_gamma !< update gamma operator with current stiffness end type tNumerics - + type(tNumerics) :: num ! numerics parameters. Better name? - + !-------------------------------------------------------------------------------------------------- ! PETSc data DM, private :: da @@ -46,7 +46,7 @@ module grid_mech_spectral_polarisation ! common pointwise data real(pReal), private, dimension(:,:,:,:,:), allocatable :: & F_lastInc, & !< field of previous compatible deformation gradients - F_tau_lastInc, & !< field of previous incompatible deformation gradient + F_tau_lastInc, & !< field of previous incompatible deformation gradient Fdot, & !< field of assumed rate of compatible deformation gradient F_tauDot !< field of assumed rate of incopatible deformation gradient @@ -58,25 +58,25 @@ module grid_mech_spectral_polarisation 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 - + character(len=pStringLen), private :: incInfo !< time and increment information real(pReal), private, dimension(3,3,3,3) :: & - C_volAvg = 0.0_pReal, & !< current volume average stiffness + 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 C_minMaxAvgLastInc = 0.0_pReal, & !< previous (min+max)/2 stiffness S = 0.0_pReal, & !< current compliance (filled up with zeros) C_scale = 0.0_pReal, & S_scale = 0.0_pReal - + real(pReal), private :: & err_BC, & !< deviation from stress BC err_curl, & !< RMS of curl of F err_div !< RMS of div of P - + integer, private :: & totalIter = 0 !< total iteration in current increment - + public :: & grid_mech_spectral_polarisation_init, & grid_mech_spectral_polarisation_solution, & @@ -90,26 +90,26 @@ contains !> @brief allocates all necessary fields and fills them with data, potentially from restart info !-------------------------------------------------------------------------------------------------- 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 - + PetscErrorCode :: ierr PetscScalar, pointer, dimension(:,:,:,:) :: & FandF_tau, & ! overall pointer to solution data F, & ! specific (sub)pointer F_tau ! specific (sub)pointer - PetscInt, dimension(0:worldsize-1) :: localK + PetscInt, dimension(0:worldsize-1) :: localK integer(HID_T) :: fileHandle, groupHandle integer :: fileUnit character(len=pStringLen) :: fileName - + write(6,'(/,a)') ' <<<+- grid_mech_spectral_polarisation init -+>>>'; flush(6) - + write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006' - + num%update_gamma = config_numerics%getInt('update_gamma',defaultVal=0) > 0 !-------------------------------------------------------------------------------------------------- @@ -125,11 +125,11 @@ subroutine grid_mech_spectral_polarisation_init allocate(Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) allocate(F_tau_lastInc(3,3,grid(1),grid(2),grid3),source = 0.0_pReal) allocate(F_tauDot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - + !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) - call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) + call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) localK = 0 localK(worldrank) = grid3 call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) @@ -147,24 +147,24 @@ subroutine grid_mech_spectral_polarisation_init call DMsetUp(da,ierr); CHKERRQ(ierr) call DMcreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 18, i.e. every def grad tensor) call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector - CHKERRQ(ierr) + CHKERRQ(ierr) call SNESsetConvergenceTest(snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "converged" CHKERRQ(ierr) call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments !-------------------------------------------------------------------------------------------------- -! init fields +! init fields call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! places pointer on PETSc data F => FandF_tau(0: 8,:,:,:) F_tau => FandF_tau(9:17,:,:,:) - + restartRead: if (interface_restartInc > 0) then write(6,'(/,a,i0,a)') ' reading restart data of increment ', interface_restartInc, ' from file' write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5' fileHandle = HDF5_openFile(fileName) groupHandle = HDF5_openGroup(fileHandle,'solver') - + call HDF5_read(groupHandle,F_aim, 'F_aim') call HDF5_read(groupHandle,F_aim_lastInc,'F_aim_lastInc') call HDF5_read(groupHandle,F_aimDot, 'F_aimDot') @@ -172,21 +172,21 @@ subroutine grid_mech_spectral_polarisation_init call HDF5_read(groupHandle,F_lastInc, 'F_lastInc') call HDF5_read(groupHandle,F_tau, 'F_tau') call HDF5_read(groupHandle,F_tau_lastInc,'F_tau_lastInc') - + elseif (interface_restartInc == 0) then restartRead F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) F_tau = 2.0_pReal*F F_tau_lastInc = 2.0_pReal*F_lastInc 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 reshape(F,shape(F_lastInc)), & ! target F 0.0_pReal) ! time increment call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! deassociate pointer - + restartRead2: if (interface_restartInc > 0) then write(6,'(/,a,i0,a)') ' reading more restart data of increment ', interface_restartInc, ' from file' call HDF5_read(groupHandle,C_volAvg, 'C_volAvg') @@ -200,12 +200,12 @@ subroutine grid_mech_spectral_polarisation_init call MPI_File_read(fileUnit,C_minMaxAvg,81,MPI_DOUBLE,MPI_STATUS_IGNORE,ierr) call MPI_File_close(fileUnit,ierr) endif restartRead2 - + call utilities_updateGamma(C_minMaxAvg) call utilities_saveReferenceStiffness C_scale = C_minMaxAvg S_scale = math_invSym3333(C_minMaxAvg) - + end subroutine grid_mech_spectral_polarisation_init @@ -229,9 +229,9 @@ function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old, solution !-------------------------------------------------------------------------------------------------- ! PETSc Data - PetscErrorCode :: ierr + PetscErrorCode :: ierr SNESConvergedReason :: reason - + incInfo = incInfoIn !-------------------------------------------------------------------------------------------------- @@ -241,10 +241,10 @@ function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old, call utilities_updateGamma(C_minMaxAvg) C_scale = C_minMaxAvg S_scale = math_invSym3333(C_minMaxAvg) - endif + endif !-------------------------------------------------------------------------------------------------- -! set module wide available data +! set module wide available data params%stress_mask = stress_BC%maskFloat params%stress_BC = stress_BC%values params%rotation_BC = rotation_BC @@ -252,13 +252,13 @@ function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old, params%timeincOld = timeinc_old !-------------------------------------------------------------------------------------------------- -! solve BVP +! solve BVP call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! check convergence call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr) - + solution%converged = reason > 0 solution%iterationsNeeded = totalIter solution%termIll = terminallyIll @@ -299,7 +299,7 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc if (cutBack) then C_volAvg = C_volAvgLastInc C_minMaxAvg = C_minMaxAvgLastInc - else + else C_volAvgLastInc = C_volAvg C_minMaxAvgLastInc = C_minMaxAvg @@ -321,13 +321,13 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc Fdot = utilities_calculateRate(guess, & F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, & - rotation_BC%rotTensor2(F_aimDot,active=.true.)) + rotation_BC%rotate(F_aimDot,active=.true.)) F_tauDot = utilities_calculateRate(guess, & F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid3]), timeinc_old, & - rotation_BC%rotTensor2(F_aimDot,active=.true.)) + rotation_BC%rotate(F_aimDot,active=.true.)) F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3]) - + materialpoint_F0 = reshape(F,[3,3,1,product(grid(1:2))*grid3]) endif @@ -335,7 +335,7 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc ! 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 - rotation_BC%rotTensor2(F_aim,active=.true.)),& + rotation_BC%rotate(F_aim,active=.true.)),& [9,grid(1),grid(2),grid3]) if (guess) then F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), & @@ -351,7 +351,7 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(1:9,i,j,k) enddo; enddo; enddo endif - + call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) end subroutine grid_mech_spectral_polarisation_forward @@ -364,7 +364,7 @@ subroutine grid_mech_spectral_polarisation_updateCoords PetscErrorCode :: ierr PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau - + call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) call utilities_updateCoords(FandF_tau(0:8,:,:,:)) call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) @@ -391,7 +391,7 @@ subroutine grid_mech_spectral_polarisation_restartWrite write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5' fileHandle = HDF5_openFile(fileName,'w') groupHandle = HDF5_addGroup(fileHandle,'solver') - + call HDF5_write(groupHandle,F_aim, 'F_aim') call HDF5_write(groupHandle,F_aim_lastInc,'F_aim_lastInc') call HDF5_write(groupHandle,F_aimDot, 'F_aimDot') @@ -405,9 +405,9 @@ subroutine grid_mech_spectral_polarisation_restartWrite call HDF5_closeGroup(groupHandle) call HDF5_closeFile(fileHandle) - + if(num%update_gamma) call utilities_saveReferenceStiffness - + call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) end subroutine grid_mech_spectral_polarisation_restartWrite @@ -417,13 +417,13 @@ end subroutine grid_mech_spectral_polarisation_restartWrite !> @brief convergence check !-------------------------------------------------------------------------------------------------- subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dummy,ierr) - + SNES :: snes_local PetscInt, intent(in) :: PETScIter PetscReal, intent(in) :: & devNull1, & devNull2, & - devNull3 + devNull3 SNESConvergedReason :: reason PetscObject :: dummy PetscErrorCode :: ierr @@ -431,11 +431,11 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm curlTol, & divTol, & BCTol - + curlTol = max(maxval(abs(F_aim-math_I3))*err_curl_tolRel ,err_curl_tolAbs) divTol = max(maxval(abs(P_av)) *err_div_tolRel ,err_div_tolAbs) BCTol = max(maxval(abs(P_av)) *err_stress_tolRel,err_stress_tolAbs) - + if ((totalIter >= itmin .and. & all([ err_div /divTol, & err_curl/curlTol, & @@ -456,9 +456,9 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', & err_curl/curlTol,' (',err_curl,' -, tol = ',curlTol,')' write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', & - err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' + err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' write(6,'(/,a)') ' ===========================================================================' - flush(6) + flush(6) end subroutine converged @@ -498,7 +498,7 @@ subroutine formResidual(in, FandF_tau, & F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - + call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) @@ -510,14 +510,14 @@ subroutine formResidual(in, FandF_tau, & write(6,'(1x,a,3(a,i0))') trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotTensor2(F_aim,active=.true.)) + ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & ' deformation gradient aim =', transpose(F_aim) flush(6) endif newIteration !-------------------------------------------------------------------------------------------------- -! +! tensorField_real = 0.0_pReal do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) tensorField_real(1:3,1:3,i,j,k) = & @@ -525,15 +525,15 @@ subroutine formResidual(in, FandF_tau, & polarAlpha*matmul(F(1:3,1:3,i,j,k), & math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3)) enddo; enddo; enddo - + !-------------------------------------------------------------------------------------------------- -! doing convolution in Fourier space +! doing convolution in Fourier space call utilities_FFTtensorForward - call utilities_fourierGammaConvolution(params%rotation_BC%rotTensor2(polarBeta*F_aim,active=.true.)) + call utilities_fourierGammaConvolution(params%rotation_BC%rotate(polarBeta*F_aim,active=.true.)) call utilities_FFTtensorBackward !-------------------------------------------------------------------------------------------------- -! constructing residual +! constructing residual residual_F_tau = polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) !-------------------------------------------------------------------------------------------------- @@ -541,13 +541,13 @@ subroutine formResidual(in, FandF_tau, & call utilities_constitutiveResponse(residual_F, & ! "residuum" gets field of first PK stress (to save memory) P_av,C_volAvg,C_minMaxAvg, & F - residual_F_tau/polarBeta,params%timeinc,params%rotation_BC) - call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) - + call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) + !-------------------------------------------------------------------------------------------------- ! 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%rotTensor2(F_av)) + & + -params%rotation_BC%rotate(F_av)) + & params%stress_mask * (P_av-params%stress_BC))) ! mask = 0.0 for no bc ! calculate divergence tensorField_real = 0.0_pReal @@ -566,7 +566,7 @@ subroutine formResidual(in, FandF_tau, & math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) & + residual_F_tau(1:3,1:3,i,j,k) enddo; enddo; enddo - + !-------------------------------------------------------------------------------------------------- ! calculating curl tensorField_real = 0.0_pReal diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index cd064b1f8..87e906727 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -18,12 +18,12 @@ module spectral_utilities use config use discretization use homogenization - + implicit none private - + include 'fftw3-mpi.f03' - + !-------------------------------------------------------------------------------------------------- ! field labels information enum, bind(c) @@ -109,8 +109,8 @@ module spectral_utilities real(pReal) :: timeinc real(pReal) :: timeincOld end type tSolutionParams - - type, private :: tNumerics + + type, private :: tNumerics real(pReal) :: & FFTW_timelimit !< timelimit for FFTW plan creation, see www.fftw.org integer :: & @@ -122,7 +122,7 @@ module spectral_utilities FFTW_plan_mode, & !< FFTW plan mode, see www.fftw.org PETSc_options end type tNumerics - + type(tNumerics) :: num ! numerics parameters. Better name? enum, bind(c) @@ -189,18 +189,18 @@ subroutine utilities_init scalarSize = 1_C_INTPTR_T, & vecSize = 3_C_INTPTR_T, & tensorSize = 9_C_INTPTR_T - + write(6,'(/,a)') ' <<<+- spectral_utilities init -+>>>' - + write(6,'(/,a)') ' Diehl, Diploma Thesis TU München, 2010' write(6,'(a)') ' https://doi.org/10.13140/2.1.3234.3840' - + write(6,'(/,a)') ' Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013' write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2012.09.012' - + write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006' - + write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, 2019' write(6,'(a)') ' https://doi.org/10.1007/978-981-10-6855-3_80' @@ -209,34 +209,34 @@ subroutine utilities_init debugGeneral = iand(debug_level(debug_SPECTRAL),debug_LEVELBASIC) /= 0 debugRotation = iand(debug_level(debug_SPECTRAL),debug_SPECTRALROTATION) /= 0 debugPETSc = iand(debug_level(debug_SPECTRAL),debug_SPECTRALPETSC) /= 0 - + if(debugPETSc) write(6,'(3(/,a),/)') & ' Initializing PETSc with debug options: ', & trim(PETScDebug), & ' add more using the PETSc_Options keyword in numerics.config '; flush(6) - + call PETScOptionsClear(PETSC_NULL_OPTIONS,ierr) CHKERRQ(ierr) if(debugPETSc) call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),ierr) CHKERRQ(ierr) call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) CHKERRQ(ierr) - + grid1Red = grid(1)/2 + 1 wgt = 1.0/real(product(grid),pReal) - + write(6,'(/,a,3(i12 ))') ' grid a b c: ', grid write(6,'(a,3(es12.5))') ' size x y z: ', geomSize - + num%memory_efficient = config_numerics%getInt ('memory_efficient', defaultVal=1) > 0 num%FFTW_timelimit = config_numerics%getFloat ('fftw_timelimit', defaultVal=-1.0_pReal) num%divergence_correction = config_numerics%getInt ('divergence_correction', defaultVal=2) num%spectral_derivative = config_numerics%getString('spectral_derivative', defaultVal='continuous') num%FFTW_plan_mode = config_numerics%getString('fftw_plan_mode', defaultVal='FFTW_MEASURE') - + if (num%divergence_correction < 0 .or. num%divergence_correction > 2) & call IO_error(301,ext_msg='divergence_correction') - + select case (num%spectral_derivative) case ('continuous') spectral_derivative_ID = DERIVATIVE_CONTINUOUS_ID @@ -265,8 +265,8 @@ subroutine utilities_init else scaledGeomSize = geomSize endif - - + + select case(IO_lc(num%FFTW_plan_mode)) ! setting parameters for the plan creation of FFTW. Basically a translation from fftw3.f case('fftw_estimate') ! ordered from slow execution (but fast plan creation) to fast execution FFTW_planner_flag = FFTW_ESTIMATE @@ -285,7 +285,7 @@ subroutine utilities_init ! general initialization of FFTW (see manual on fftw.org for more details) if (pReal /= C_DOUBLE .or. kind(1) /= C_INT) call IO_error(0,ext_msg='Fortran to C') ! check for correct precision in C call fftw_set_timelimit(num%FFTW_timelimit) ! set timelimit for plan creation - + if (debugGeneral) write(6,'(/,a)') ' FFTW initialized'; flush(6) !-------------------------------------------------------------------------------------------------- @@ -295,19 +295,19 @@ subroutine utilities_init PETSC_COMM_WORLD, local_K, local_K_offset) allocate (xi1st (3,grid1Red,grid(2),grid3),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for first derivatives, only half the size for first dimension allocate (xi2nd (3,grid1Red,grid(2),grid3),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for second derivatives, only half the size for first dimension - + tensorField = fftw_alloc_complex(tensorSize*alloc_local) call c_f_pointer(tensorField, tensorField_real, [3_C_INTPTR_T,3_C_INTPTR_T, & 2_C_INTPTR_T*(gridFFTW(1)/2_C_INTPTR_T + 1_C_INTPTR_T),gridFFTW(2),local_K]) ! place a pointer for a real tensor representation call c_f_pointer(tensorField, tensorField_fourier, [3_C_INTPTR_T,3_C_INTPTR_T, & gridFFTW(1)/2_C_INTPTR_T + 1_C_INTPTR_T , gridFFTW(2),local_K]) ! place a pointer for a fourier tensor representation - + vectorField = fftw_alloc_complex(vecSize*alloc_local) call c_f_pointer(vectorField, vectorField_real, [3_C_INTPTR_T,& 2_C_INTPTR_T*(gridFFTW(1)/2_C_INTPTR_T + 1_C_INTPTR_T),gridFFTW(2),local_K]) ! place a pointer for a real vector representation call c_f_pointer(vectorField, vectorField_fourier,[3_C_INTPTR_T,& gridFFTW(1)/2_C_INTPTR_T + 1_C_INTPTR_T, gridFFTW(2),local_K]) ! place a pointer for a fourier vector representation - + scalarField = fftw_alloc_complex(scalarSize*alloc_local) ! allocate data for real representation (no in place transform) call c_f_pointer(scalarField, scalarField_real, & [2_C_INTPTR_T*(gridFFTW(1)/2_C_INTPTR_T + 1),gridFFTW(2),local_K]) ! place a pointer for a real scalar representation @@ -371,7 +371,7 @@ subroutine utilities_init xi1st(1:3,i,j,k-grid3Offset) = xi2nd(1:3,i,j,k-grid3Offset) endwhere enddo; enddo; enddo - + if(num%memory_efficient) then ! allocate just single fourth order tensor allocate (gamma_hat(3,3,3,3,1,1,1), source = cmplx(0.0_pReal,0.0_pReal,pReal)) else ! precalculation of gamma_hat field @@ -388,7 +388,7 @@ end subroutine utilities_init !> In case of an on-the-fly calculation, only the reference stiffness is updated. !--------------------------------------------------------------------------------------------------- subroutine utilities_updateGamma(C) - + real(pReal), intent(in), dimension(3,3,3,3) :: C !< input stiffness to store as reference stiffness complex(pReal), dimension(3,3) :: temp33_complex, xiDyad_cmplx real(pReal), dimension(6,6) :: A, A_inv @@ -396,9 +396,9 @@ subroutine utilities_updateGamma(C) i, j, k, & l, m, n, o logical :: err - + C_ref = C - + if(.not. num%memory_efficient) then gamma_hat = cmplx(0.0_pReal,0.0_pReal,pReal) ! for the singular point and any non invertible A do k = grid3Offset+1, grid3Offset+grid3; do j = 1, grid(2); do i = 1, grid1Red @@ -419,7 +419,7 @@ subroutine utilities_updateGamma(C) endif enddo; enddo; enddo endif - + end subroutine utilities_updateGamma @@ -501,17 +501,17 @@ end subroutine utilities_FFTvectorBackward !> @brief doing convolution gamma_hat * field_real, ensuring that average value = fieldAim !-------------------------------------------------------------------------------------------------- subroutine utilities_fourierGammaConvolution(fieldAim) - + real(pReal), intent(in), dimension(3,3) :: fieldAim !< desired average value of the field after convolution complex(pReal), dimension(3,3) :: temp33_complex, xiDyad_cmplx real(pReal), dimension(6,6) :: A, A_inv - + integer :: & i, j, k, & l, m, n, o logical :: err - - + + write(6,'(/,a)') ' ... doing gamma convolution ...............................................' flush(6) @@ -531,7 +531,7 @@ subroutine utilities_fourierGammaConvolution(fieldAim) temp33_complex = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal) forall(l=1:3, m=1:3, n=1:3, o=1:3) & gamma_hat(l,m,n,o,1,1,1) = temp33_complex(l,n)*conjg(-xi1st(o,i,j,k))*xi1st(m,i,j,k) - else + else gamma_hat(1:3,1:3,1:3,1:3,1,1,1) = cmplx(0.0_pReal,0.0_pReal,pReal) endif forall(l = 1:3, m = 1:3) & @@ -546,7 +546,7 @@ subroutine utilities_fourierGammaConvolution(fieldAim) tensorField_fourier(1:3,1:3,i,j,k) = temp33_Complex enddo; enddo; enddo endif memoryEfficient - + if (grid3Offset == 0) tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim/wgt,0.0_pReal,pReal) end subroutine utilities_fourierGammaConvolution @@ -561,7 +561,7 @@ subroutine utilities_fourierGreenConvolution(D_ref, mobility_ref, deltaT) real(pReal), intent(in) :: mobility_ref, deltaT complex(pReal) :: GreenOp_hat integer :: i, j, k - + !-------------------------------------------------------------------------------------------------- ! do the actual spectral method calculation do k = 1, grid3; do j = 1, grid(2) ;do i = 1, grid1Red @@ -625,16 +625,16 @@ real(pReal) function utilities_curlRMS() integer :: i, j, k, l, ierr complex(pReal), dimension(3,3) :: curl_fourier complex(pReal), dimension(3) :: rescaledGeom - + write(6,'(/,a)') ' ... calculating curl ......................................................' flush(6) - + rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal) - + !-------------------------------------------------------------------------------------------------- ! calculating max curl criterion in Fourier space utilities_curlRMS = 0.0_pReal - + do k = 1, grid3; do j = 1, grid(2); do i = 2, grid1Red - 1 do l = 1, 3 @@ -669,7 +669,7 @@ real(pReal) function utilities_curlRMS() utilities_curlRMS = utilities_curlRMS & + sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal) ! this layer (Nyquist) does not have a conjugate complex counterpart (if grid(1) /= 1) enddo; enddo - + call MPI_Allreduce(MPI_IN_PLACE,utilities_curlRMS,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) if(ierr /=0) call IO_error(894, ext_msg='utilities_curlRMS') utilities_curlRMS = sqrt(utilities_curlRMS) * wgt @@ -688,9 +688,10 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) type(rotation), intent(in) :: rot_BC !< rotation of load frame logical, intent(in), dimension(3,3) :: mask_stress !< mask of stress BC - integer :: j, k, m, n - logical, dimension(9) :: mask_stressVector - real(pReal), dimension(9,9) :: temp99_Real + integer :: i, j + logical, dimension(9) :: mask_stressVector + logical, dimension(9,9) :: mask + real(pReal), dimension(9,9) :: temp99_real integer :: size_reduced = 0 real(pReal), dimension(:,:), allocatable :: & s_reduced, & !< reduced compliance matrix (depending on number of stress BC) @@ -698,57 +699,33 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) sTimesC !< temp variable to check inversion logical :: errmatinv character(len=pStringLen):: formatString - + mask_stressVector = reshape(transpose(mask_stress), [9]) size_reduced = count(mask_stressVector) - if(size_reduced > 0 )then - allocate (c_reduced(size_reduced,size_reduced), source =0.0_pReal) - allocate (s_reduced(size_reduced,size_reduced), source =0.0_pReal) - allocate (sTimesC(size_reduced,size_reduced), source =0.0_pReal) - temp99_Real = math_3333to99(rot_BC%rotTensor4(C)) - + if(size_reduced > 0) then + temp99_real = math_3333to99(rot_BC%rotate(C)) + if(debugGeneral) then write(6,'(/,a)') ' ... updating masked compliance ............................................' write(6,'(/,a,/,9(9(2x,f12.7,1x)/))',advance='no') ' Stiffness C (load) / GPa =',& transpose(temp99_Real)*1.0e-9_pReal flush(6) endif - k = 0 ! calculate reduced stiffness - do n = 1,9 - if(mask_stressVector(n)) then - k = k + 1 - j = 0 - do m = 1,9 - if(mask_stressVector(m)) then - j = j + 1 - c_reduced(k,j) = temp99_Real(n,m) - endif; enddo; endif; enddo - + + do i = 1,9; do j = 1,9 + mask(i,j) = mask_stressVector(i) .and. mask_stressVector(j) + enddo; enddo + c_reduced = reshape(pack(temp99_Real,mask),[size_reduced,size_reduced]) + + allocate(s_reduced,mold = c_reduced) call math_invert(s_reduced, errmatinv, c_reduced) ! invert reduced stiffness if (any(IEEE_is_NaN(s_reduced))) errmatinv = .true. if (errmatinv) call IO_error(error_ID=400,ext_msg='utilities_maskedCompliance') - temp99_Real = 0.0_pReal ! fill up compliance with zeros - k = 0 - do n = 1,9 - if(mask_stressVector(n)) then - k = k + 1 - j = 0 - do m = 1,9 - if(mask_stressVector(m)) then - j = j + 1 - temp99_Real(n,m) = s_reduced(k,j) - endif; enddo; endif; enddo !-------------------------------------------------------------------------------------------------- ! check if inversion was successful sTimesC = matmul(c_reduced,s_reduced) - do m=1, size_reduced - do n=1, size_reduced - errmatinv = errmatinv & - .or. (m==n .and. abs(sTimesC(m,n)-1.0_pReal) > 1.0e-12_pReal) & ! diagonal elements of S*C should be 1 - .or. (m/=n .and. abs(sTimesC(m,n)) > 1.0e-12_pReal) ! off-diagonal elements of S*C should be 0 - enddo - enddo + errmatinv = errmatinv .or. any(dNeq(sTimesC,math_identity2nd(size_reduced),1.0e-12_pReal)) if (debugGeneral .or. errmatinv) then write(formatString, '(i2)') size_reduced formatString = '(/,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))' @@ -757,15 +734,18 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) write(6,trim(formatString),advance='no') ' S (load) ', transpose(s_reduced) if(errmatinv) call IO_error(error_ID=400,ext_msg='utilities_maskedCompliance') endif + temp99_real = reshape(unpack(reshape(s_reduced,[size_reduced**2]),reshape(mask,[81]),0.0_pReal),[9,9]) else temp99_real = 0.0_pReal endif + + utilities_maskedCompliance = math_99to3333(temp99_Real) + if(debugGeneral) then write(6,'(/,a,/,9(9(2x,f10.5,1x)/),/)',advance='no') & ' Masked Compliance (load) * GPa =', transpose(temp99_Real)*1.0e9_pReal flush(6) endif - utilities_maskedCompliance = math_99to3333(temp99_Real) end function utilities_maskedCompliance @@ -774,9 +754,9 @@ end function utilities_maskedCompliance !> @brief calculate scalar gradient in fourier field !-------------------------------------------------------------------------------------------------- subroutine utilities_fourierScalarGradient() - + integer :: i, j, k - + do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red vectorField_fourier(1:3,i,j,k) = scalarField_fourier(i,j,k)*xi1st(1:3,i,j,k) ! ToDo: no -conjg? enddo; enddo; enddo @@ -788,9 +768,9 @@ end subroutine utilities_fourierScalarGradient !> @brief calculate vector divergence in fourier field !-------------------------------------------------------------------------------------------------- subroutine utilities_fourierVectorDivergence() - + integer :: i, j, k - + do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red scalarField_fourier(i,j,k) = sum(vectorField_fourier(1:3,i,j,k)*conjg(-xi1st(1:3,i,j,k))) enddo; enddo; enddo @@ -802,9 +782,9 @@ end subroutine utilities_fourierVectorDivergence !> @brief calculate vector gradient in fourier field !-------------------------------------------------------------------------------------------------- subroutine utilities_fourierVectorGradient() - + integer :: i, j, k, m, n - + do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red do m = 1, 3; do n = 1, 3 tensorField_fourier(m,n,i,j,k) = vectorField_fourier(m,i,j,k)*xi1st(n,i,j,k) @@ -820,7 +800,7 @@ end subroutine utilities_fourierVectorGradient subroutine utilities_fourierTensorDivergence() integer :: i, j, k - + do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red vectorField_fourier(:,i,j,k) = matmul(tensorField_fourier(:,:,i,j,k),conjg(-xi1st(:,i,j,k))) enddo; enddo; enddo @@ -833,28 +813,28 @@ end subroutine utilities_fourierTensorDivergence !-------------------------------------------------------------------------------------------------- subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& F,timeinc,rotation_BC) - + real(pReal), intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness real(pReal), intent(out), dimension(3,3) :: P_av !< average PK stress real(pReal), intent(out), dimension(3,3,grid(1),grid(2),grid3) :: P !< PK stress real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: F !< deformation gradient target real(pReal), intent(in) :: timeinc !< loading time type(rotation), intent(in), optional :: rotation_BC !< rotation of load frame - - + + integer :: & i,ierr real(pReal), dimension(3,3,3,3) :: dPdF_max, dPdF_min real(pReal) :: dPdF_norm_max, dPdF_norm_min real(pReal), dimension(2) :: valueAndRank !< pair of min/max norm of dPdF to synchronize min/max of dPdF - + write(6,'(/,a)') ' ... evaluating constitutive response ......................................' flush(6) - + materialpoint_F = reshape(F,[3,3,1,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field - + call materialpoint_stressAndItsTangent(.true.,timeinc) ! calculate P field - + P = reshape(materialpoint_P, [3,3,grid(1),grid(2),grid3]) P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) @@ -862,11 +842,11 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress (lab) / MPa =',& transpose(P_av)*1.e-6_pReal if(present(rotation_BC)) & - P_av = rotation_BC%rotTensor2(P_av) + P_av = rotation_BC%rotate(P_av) write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',& transpose(P_av)*1.e-6_pReal flush(6) - + dPdF_max = 0.0_pReal dPdF_norm_max = 0.0_pReal dPdF_min = huge(1.0_pReal) @@ -881,21 +861,21 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& dPdF_norm_min = sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal) endif end do - + valueAndRank = [dPdF_norm_max,real(worldrank,pReal)] call MPI_Allreduce(MPI_IN_PLACE,valueAndRank,1, MPI_2DOUBLE_PRECISION, MPI_MAXLOC, PETSC_COMM_WORLD, ierr) if (ierr /= 0) call IO_error(894, ext_msg='MPI_Allreduce max') call MPI_Bcast(dPdF_max,81,MPI_DOUBLE,int(valueAndRank(2)),PETSC_COMM_WORLD, ierr) if (ierr /= 0) call IO_error(894, ext_msg='MPI_Bcast max') - + valueAndRank = [dPdF_norm_min,real(worldrank,pReal)] call MPI_Allreduce(MPI_IN_PLACE,valueAndRank,1, MPI_2DOUBLE_PRECISION, MPI_MINLOC, PETSC_COMM_WORLD, ierr) if (ierr /= 0) call IO_error(894, ext_msg='MPI_Allreduce min') call MPI_Bcast(dPdF_min,81,MPI_DOUBLE,int(valueAndRank(2)),PETSC_COMM_WORLD, ierr) if (ierr /= 0) call IO_error(894, ext_msg='MPI_Bcast min') - + C_minmaxAvg = 0.5_pReal*(dPdF_max + dPdF_min) - + C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) C_volAvg = C_volAvg * wgt @@ -908,7 +888,7 @@ end subroutine utilities_constitutiveResponse !> @brief calculates forward rate, either guessing or just add delta/timeinc !-------------------------------------------------------------------------------------------------- pure function utilities_calculateRate(heterogeneous,field0,field,dt,avRate) - + real(pReal), intent(in), dimension(3,3) :: & avRate !< homogeneous addon real(pReal), intent(in) :: & @@ -920,7 +900,7 @@ pure function utilities_calculateRate(heterogeneous,field0,field,dt,avRate) field !< data of current step real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: & utilities_calculateRate - + if (heterogeneous) then utilities_calculateRate = (field-field0) / dt else @@ -971,14 +951,14 @@ pure function utilities_getFreqDerivative(k_s) complex(pReal), dimension(3) :: utilities_getFreqDerivative select case (spectral_derivative_ID) - case (DERIVATIVE_CONTINUOUS_ID) + case (DERIVATIVE_CONTINUOUS_ID) utilities_getFreqDerivative = cmplx(0.0_pReal, 2.0_pReal*PI*real(k_s,pReal)/geomSize,pReal) - case (DERIVATIVE_CENTRAL_DIFF_ID) + case (DERIVATIVE_CENTRAL_DIFF_ID) utilities_getFreqDerivative = cmplx(0.0_pReal, sin(2.0_pReal*PI*real(k_s,pReal)/real(grid,pReal)), pReal)/ & - cmplx(2.0_pReal*geomSize/real(grid,pReal), 0.0_pReal, pReal) - - case (DERIVATIVE_FWBW_DIFF_ID) + cmplx(2.0_pReal*geomSize/real(grid,pReal), 0.0_pReal, pReal) + + case (DERIVATIVE_FWBW_DIFF_ID) utilities_getFreqDerivative(1) = & cmplx(cos(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)) - 1.0_pReal, & sin(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)), pReal)* & @@ -986,7 +966,7 @@ pure function utilities_getFreqDerivative(k_s) sin(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)), pReal)* & cmplx(cos(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)) + 1.0_pReal, & sin(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)), pReal)/ & - cmplx(4.0_pReal*geomSize(1)/real(grid(1),pReal), 0.0_pReal, pReal) + cmplx(4.0_pReal*geomSize(1)/real(grid(1),pReal), 0.0_pReal, pReal) utilities_getFreqDerivative(2) = & cmplx(cos(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)) + 1.0_pReal, & sin(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)), pReal)* & @@ -994,7 +974,7 @@ pure function utilities_getFreqDerivative(k_s) sin(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)), pReal)* & cmplx(cos(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)) + 1.0_pReal, & sin(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)), pReal)/ & - cmplx(4.0_pReal*geomSize(2)/real(grid(2),pReal), 0.0_pReal, pReal) + cmplx(4.0_pReal*geomSize(2)/real(grid(2),pReal), 0.0_pReal, pReal) utilities_getFreqDerivative(3) = & cmplx(cos(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)) + 1.0_pReal, & sin(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)), pReal)* & @@ -1002,7 +982,7 @@ pure function utilities_getFreqDerivative(k_s) sin(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)), pReal)* & cmplx(cos(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)) - 1.0_pReal, & sin(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)), pReal)/ & - cmplx(4.0_pReal*geomSize(3)/real(grid(3),pReal), 0.0_pReal, pReal) + cmplx(4.0_pReal*geomSize(3)/real(grid(3),pReal), 0.0_pReal, pReal) end select end function utilities_getFreqDerivative @@ -1014,7 +994,7 @@ end function utilities_getFreqDerivative ! convolution !-------------------------------------------------------------------------------------------------- subroutine utilities_updateCoords(F) - + real(pReal), dimension(3,3,grid(1),grid(2),grid3), intent(in) :: F real(pReal), dimension(3, grid(1),grid(2),grid3) :: IPcoords real(pReal), dimension(3, grid(1),grid(2),grid3+2) :: IPfluct_padded ! Fluctuations of cell center displacement (padded along z for MPI) @@ -1040,7 +1020,7 @@ subroutine utilities_updateCoords(F) 1, 0, 1, & 1, 1, 1, & 0, 1, 1 ], [3,8]) - + step = geomSize/real(grid, pReal) !-------------------------------------------------------------------------------------------------- ! integration in Fourier space to get fluctuations of cell center discplacements @@ -1057,27 +1037,27 @@ subroutine utilities_updateCoords(F) enddo; enddo; enddo call fftw_mpi_execute_dft_c2r(planVectorBack,vectorField_fourier,vectorField_real) - + !-------------------------------------------------------------------------------------------------- ! average F if (grid3Offset == 0) Favg = real(tensorField_fourier(1:3,1:3,1,1,1),pReal)*wgt call MPI_Bcast(Favg,9,MPI_DOUBLE,0,PETSC_COMM_WORLD,ierr) if(ierr /=0) call IO_error(894, ext_msg='update_IPcoords/MPI_Bcast') - + !-------------------------------------------------------------------------------------------------- ! pad cell center fluctuations along z-direction (needed when running MPI simulation) IPfluct_padded(1:3,1:grid(1),1:grid(2),2:grid3+1) = vectorField_real(1:3,1:grid(1),1:grid(2),1:grid3) c = product(shape(IPfluct_padded(:,:,:,1))) !< amount of data to transfer rank_t = modulo(worldrank+1,worldsize) rank_b = modulo(worldrank-1,worldsize) - + ! send bottom layer to process below call MPI_Isend(IPfluct_padded(:,:,:,2), c,MPI_DOUBLE,rank_b,0,PETSC_COMM_WORLD,r,ierr) if(ierr /=0) call IO_error(894, ext_msg='update_IPcoords/MPI_Isend') 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) - + ! 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) @@ -1085,9 +1065,9 @@ subroutine utilities_updateCoords(F) 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) - + !-------------------------------------------------------------------------------------------------- - ! calculate nodal displacements + ! calculate nodal displacements nodeCoords = 0.0_pReal do k = 0,grid3; do j = 0,grid(2); do i = 0,grid(1) nodeCoords(1:3,i+1,j+1,k+1) = matmul(Favg,step*(real([i,j,k+grid3Offset],pReal))) @@ -1097,17 +1077,17 @@ subroutine utilities_updateCoords(F) + IPfluct_padded(1:3,modulo(me(1)-1,grid(1))+1,modulo(me(2)-1,grid(2))+1,me(3)+1)*0.125_pReal enddo averageFluct enddo; enddo; enddo - + !-------------------------------------------------------------------------------------------------- ! calculate cell center displacements do k = 1,grid3; do j = 1,grid(2); do i = 1,grid(1) IPcoords(1:3,i,j,k) = vectorField_real(1:3,i,j,k) & + matmul(Favg,step*real([i,j,k+grid3Offset]-0.5_pReal,pReal)) enddo; enddo; enddo - + call discretization_setNodeCoords(reshape(NodeCoords,[3,(grid(1)+1)*(grid(2)+1)*(grid3+1)])) call discretization_setIPcoords (reshape(IPcoords, [3,grid(1)*grid(2)*grid3])) - + end subroutine utilities_updateCoords @@ -1115,7 +1095,7 @@ end subroutine utilities_updateCoords !> @brief Write out the current reference stiffness for restart. !--------------------------------------------------------------------------------------------------- subroutine utilities_saveReferenceStiffness - + integer :: & fileUnit @@ -1125,7 +1105,7 @@ subroutine utilities_saveReferenceStiffness write(fileUnit) C_ref close(fileUnit) endif - + end subroutine utilities_saveReferenceStiffness end module spectral_utilities diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 08e02290b..9d935866c 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -146,9 +146,7 @@ subroutine homogenization_init !-------------------------------------------------------------------------------------------------- ! allocate and initialize global variables allocate(materialpoint_dPdF(3,3,3,3,discretization_nIP,discretization_nElem), source=0.0_pReal) - allocate(materialpoint_F0(3,3,discretization_nIP,discretization_nElem), source=0.0_pReal) materialpoint_F0 = spread(spread(math_I3,3,discretization_nIP),4,discretization_nElem) ! initialize to identity - allocate(materialpoint_F(3,3,discretization_nIP,discretization_nElem), source=0.0_pReal) materialpoint_F = materialpoint_F0 ! initialize to identity allocate(materialpoint_subF0(3,3,discretization_nIP,discretization_nElem), source=0.0_pReal) allocate(materialpoint_subF(3,3,discretization_nIP,discretization_nElem), source=0.0_pReal) @@ -333,12 +331,10 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) !$OMP FLUSH(terminallyIll) if (.not. terminallyIll) then ! so first signals terminally ill... !$OMP CRITICAL (write2out) - write(6,*) 'Integration point ', i,' at element ', e, ' terminally ill' + write(6,*) 'Integration point ', i,' at element ', e, ' terminally ill' !$OMP END CRITICAL (write2out) endif - !$OMP CRITICAL (setTerminallyIll) - terminallyIll = .true. ! ...and kills all others - !$OMP END CRITICAL (setTerminallyIll) + terminallyIll = .true. ! ...and kills all others else ! cutback makes sense materialpoint_subStep(i,e) = subStepSizeHomog * materialpoint_subStep(i,e) ! crystallite had severe trouble, so do a significant cutback diff --git a/src/lattice.f90 b/src/lattice.f90 index 3d624ba48..ab2f43284 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -15,14 +15,14 @@ module lattice implicit none private - + ! BEGIN DEPRECATED integer, parameter, public :: & LATTICE_maxNcleavageFamily = 3 !< max # of transformation system families over lattice structures - + integer, allocatable, dimension(:,:), protected, public :: & lattice_NcleavageSystem !< total # of transformation systems in each family - + real(pReal), allocatable, dimension(:,:,:,:,:), protected, public :: & lattice_Scleavage !< Schmid matrices for cleavage systems ! END DEPRECATED @@ -32,16 +32,16 @@ module lattice ! face centered cubic integer, dimension(2), parameter :: & LATTICE_FCC_NSLIPSYSTEM = [12, 6] !< # of slip systems per family for fcc - + integer, dimension(1), parameter :: & LATTICE_FCC_NTWINSYSTEM = [12] !< # of twin systems per family for fcc - + integer, dimension(1), parameter :: & LATTICE_FCC_NTRANSSYSTEM = [12] !< # of transformation systems per family for fcc - + integer, dimension(2), parameter :: & LATTICE_FCC_NCLEAVAGESYSTEM = [3, 4] !< # of cleavage systems per family for fcc - + integer, parameter :: & #ifndef __PGI LATTICE_FCC_NSLIP = sum(LATTICE_FCC_NSLIPSYSTEM), & !< total # of slip systems for fcc @@ -78,7 +78,7 @@ module lattice 0, 1, 1, 0, 1,-1, & 0, 1,-1, 0, 1, 1 & ],pReal),shape(LATTICE_FCC_SYSTEMSLIP)) !< Slip system <110>{111} directions. Sorted according to Eisenlohr & Hantcherli - + real(pReal), dimension(3+3,LATTICE_FCC_NTWIN), parameter :: & LATTICE_FCC_SYSTEMTWIN = reshape(real( [& -2, 1, 1, 1, 1, 1, & @@ -94,7 +94,7 @@ module lattice -1,-2,-1, -1, 1,-1, & -1, 1, 2, -1, 1,-1 & ],pReal),shape(LATTICE_FCC_SYSTEMTWIN)) !< Twin system <112>{111} directions. Sorted according to Eisenlohr & Hantcherli - + integer, dimension(2,LATTICE_FCC_NTWIN), parameter, public :: & LATTICE_FCC_TWINNUCLEATIONSLIPPAIR = reshape( [& 2,3, & @@ -110,7 +110,7 @@ module lattice 10,12, & 10,11 & ],shape(LATTICE_FCC_TWINNUCLEATIONSLIPPAIR)) - + real(pReal), dimension(3+3,LATTICE_FCC_NCLEAVAGE), parameter :: & LATTICE_FCC_SYSTEMCLEAVAGE = reshape(real([& ! Cleavage direction Plane normal @@ -122,18 +122,18 @@ module lattice -1, 0,-1, 1,-1,-1, & 0, 1, 1, -1, 1,-1 & ],pReal),shape(LATTICE_FCC_SYSTEMCLEAVAGE)) - + !-------------------------------------------------------------------------------------------------- ! body centered cubic integer, dimension(2), parameter :: & LATTICE_BCC_NSLIPSYSTEM = [12, 12] !< # of slip systems per family for bcc - + integer, dimension(1), parameter :: & LATTICE_BCC_NTWINSYSTEM = [12] !< # of twin systems per family for bcc - + integer, dimension(2), parameter :: & LATTICE_BCC_NCLEAVAGESYSTEM = [3, 6] !< # of cleavage systems per family for bcc - + integer, parameter :: & #ifndef __PGI LATTICE_BCC_NSLIP = sum(LATTICE_BCC_NSLIPSYSTEM), & !< total # of slip systems for bcc @@ -175,7 +175,7 @@ module lattice -1, 1, 1, 1,-1, 2, & 1, 1, 1, 1, 1,-2 & ],pReal),shape(LATTICE_BCC_SYSTEMSLIP)) - + real(pReal), dimension(3+3,LATTICE_BCC_NTWIN), parameter :: & LATTICE_BCC_SYSTEMTWIN = reshape(real([& ! Twin system <111>{112} @@ -191,8 +191,8 @@ module lattice 1,-1, 1, -1, 1, 2, & -1, 1, 1, 1,-1, 2, & 1, 1, 1, 1, 1,-2 & - ],pReal),shape(LATTICE_BCC_SYSTEMTWIN)) - + ],pReal),shape(LATTICE_BCC_SYSTEMTWIN)) + real(pReal), dimension(3+3,LATTICE_BCC_NCLEAVAGE), parameter :: & LATTICE_BCC_SYSTEMCLEAVAGE = reshape(real([& ! Cleavage direction Plane normal @@ -206,18 +206,18 @@ module lattice -1, 1, 1, 1, 1, 0, & 1, 1, 1, -1, 1, 0 & ],pReal),shape(LATTICE_BCC_SYSTEMCLEAVAGE)) - + !-------------------------------------------------------------------------------------------------- ! hexagonal integer, dimension(6), parameter :: & LATTICE_HEX_NSLIPSYSTEM = [3, 3, 3, 6, 12, 6] !< # of slip systems per family for hex - + integer, dimension(4), parameter :: & LATTICE_HEX_NTWINSYSTEM = [6, 6, 6, 6] !< # of slip systems per family for hex - + integer, dimension(1), parameter :: & LATTICE_HEX_NCLEAVAGESYSTEM = [3] !< # of cleavage systems per family for hex - + integer, parameter :: & #ifndef __PGI LATTICE_HEX_NSLIP = sum(LATTICE_HEX_NSLIPSYSTEM), & !< total # of slip systems for hex @@ -265,14 +265,14 @@ module lattice -1, 2, -1, 3, 1, -1, 0, 1, & -2, 1, 1, 3, 1, -1, 0, 1, & ! pyramidal system: c+a slip <11.3>{-1-1.2} -- as for hexagonal ice (Castelnau et al. 1996, similar to twin system found below) - -1, -1, 2, 3, 1, 1, -2, 2, & ! <11.3>{-1-1.2} shear = 2((c/a)^2-2)/(3 c/a) + -1, -1, 2, 3, 1, 1, -2, 2, & ! <11.3>{-1-1.2} shear = 2((c/a)^2-2)/(3 c/a) 1, -2, 1, 3, -1, 2, -1, 2, & 2, -1, -1, 3, -2, 1, 1, 2, & 1, 1, -2, 3, -1, -1, 2, 2, & -1, 2, -1, 3, 1, -2, 1, 2, & -2, 1, 1, 3, 2, -1, -1, 2 & ],pReal),shape(LATTICE_HEX_SYSTEMSLIP)) !< slip systems for hex, sorted by P. Eisenlohr CCW around starting next to a_1 axis - + real(pReal), dimension(4+4,LATTICE_HEX_NTWIN), parameter :: & LATTICE_HEX_SYSTEMTWIN = reshape(real([& ! Compression or Tension = f(twinning shear=f(c/a)) for each metal ! (according to Yoo 1981) @@ -304,7 +304,7 @@ module lattice 1, -2, 1, -3, 1, -2, 1, 2, & 2, -1, -1, -3, 2, -1, -1, 2 & ],pReal),shape(LATTICE_HEX_SYSTEMTWIN)) !< twin systems for hex, sorted by P. Eisenlohr CCW around starting next to a_1 axis - + real(pReal), dimension(4+4,LATTICE_HEX_NCLEAVAGE), parameter :: & LATTICE_HEX_SYSTEMCLEAVAGE = reshape(real([& ! Cleavage direction Plane normal @@ -312,20 +312,20 @@ module lattice 0, 0, 0, 1, 2,-1,-1, 0, & 0, 0, 0, 1, 0, 1,-1, 0 & ],pReal),shape(LATTICE_HEX_SYSTEMCLEAVAGE)) - - + + !-------------------------------------------------------------------------------------------------- ! body centered tetragonal integer, dimension(13), parameter :: & LATTICE_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 :: & #ifndef __PGI LATTICE_BCT_NSLIP = sum(LATTICE_BCT_NSLIPSYSTEM) !< total # of slip systems for bct #else LATTICE_BCT_NSLIP = 52 #endif - + real(pReal), dimension(3+3,LATTICE_BCT_NSLIP), parameter :: & LATTICE_BCT_SYSTEMSLIP = reshape(real([& ! Slip direction Plane normal @@ -395,12 +395,12 @@ module lattice -1, 1, 1, -1,-2, 1, & 1, 1, 1, 1,-2, 1 & ],pReal),shape(LATTICE_BCT_SYSTEMSLIP)) !< slip systems for bct sorted by Bieler - + !-------------------------------------------------------------------------------------------------- ! isotropic integer, dimension(1), parameter :: & LATTICE_ISO_NCLEAVAGESYSTEM = [3] !< # of cleavage systems per family for iso - + integer, parameter :: & #ifndef __PGI LATTICE_ISO_NCLEAVAGE = sum(LATTICE_ISO_NCLEAVAGESYSTEM) !< total # of cleavage systems for iso @@ -415,13 +415,13 @@ module lattice 0, 0, 1, 0, 1, 0, & 1, 0, 0, 0, 0, 1 & ],pReal),shape(LATTICE_ISO_SYSTEMCLEAVAGE)) - - + + !-------------------------------------------------------------------------------------------------- ! orthorhombic integer, dimension(3), parameter :: & LATTICE_ORT_NCLEAVAGESYSTEM = [1, 1, 1] !< # of cleavage systems per family for ortho - + integer, parameter :: & #ifndef __PGI LATTICE_ORT_NCLEAVAGE = sum(LATTICE_ORT_NCLEAVAGESYSTEM) !< total # of cleavage systems for ortho @@ -436,23 +436,23 @@ module lattice 0, 0, 1, 0, 1, 0, & 1, 0, 0, 0, 0, 1 & ],pReal),shape(LATTICE_ORT_SYSTEMCLEAVAGE)) - - - + + + ! BEGIN DEPRECATED integer, parameter, public :: & LATTICE_maxNcleavage = max(LATTICE_fcc_Ncleavage,LATTICE_bcc_Ncleavage, & LATTICE_hex_Ncleavage, & LATTICE_iso_Ncleavage,LATTICE_ort_Ncleavage) ! END DEPRECATED - + real(pReal), dimension(:,:,:), allocatable, public, protected :: & lattice_C66 real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: & lattice_C3333 real(pReal), dimension(:), allocatable, public, protected :: & lattice_mu, lattice_nu - + ! SHOULD NOT BE PART OF LATTICE BEGIN real(pReal), dimension(:,:,:,:), allocatable, public, protected :: & ! with higher-order parameters (e.g. temperature-dependent) lattice_thermalExpansion33 @@ -465,7 +465,7 @@ module lattice lattice_specificHeat, & lattice_referenceTemperature ! SHOULD NOT BE PART OF LATTICE END - + enum, bind(c) enumerator :: LATTICE_undefined_ID, & LATTICE_iso_ID, & @@ -475,18 +475,18 @@ module lattice LATTICE_bct_ID, & LATTICE_ort_ID end enum - + integer(kind(LATTICE_undefined_ID)), dimension(:), allocatable, public, protected :: & lattice_structure - + interface lattice_forestProjection_edge module procedure slipProjection_transverse end interface lattice_forestProjection_edge - + interface lattice_forestProjection_screw module procedure slipProjection_direction end interface lattice_forestProjection_screw - + public :: & lattice_init, & LATTICE_iso_ID, & @@ -516,29 +516,29 @@ module lattice lattice_slip_transverse, & lattice_labels_slip, & lattice_labels_twin - + contains !-------------------------------------------------------------------------------------------------- !> @brief Module initialization !-------------------------------------------------------------------------------------------------- subroutine lattice_init - - integer :: Nphases, p + + integer :: Nphases, p,i character(len=pStringLen) :: & tag = '' real(pReal) :: CoverA real(pReal), dimension(:), allocatable :: & temp - + write(6,'(/,a)') ' <<<+- lattice init -+>>>' - + Nphases = size(config_phase) - + allocate(lattice_structure(Nphases),source = LATTICE_undefined_ID) allocate(lattice_C66(6,6,Nphases), source=0.0_pReal) allocate(lattice_C3333(3,3,3,3,Nphases), source=0.0_pReal) - + allocate(lattice_thermalExpansion33 (3,3,3,Nphases), source=0.0_pReal) ! constant, linear, quadratic coefficients allocate(lattice_thermalConductivity33 (3,3,Nphases), source=0.0_pReal) allocate(lattice_damageDiffusion33 (3,3,Nphases), source=0.0_pReal) @@ -546,32 +546,16 @@ subroutine lattice_init allocate(lattice_massDensity ( Nphases), source=0.0_pReal) allocate(lattice_specificHeat ( Nphases), source=0.0_pReal) allocate(lattice_referenceTemperature ( Nphases), source=300.0_pReal) - + allocate(lattice_mu(Nphases), source=0.0_pReal) allocate(lattice_nu(Nphases), source=0.0_pReal) - - + allocate(lattice_Scleavage(3,3,3,lattice_maxNcleavage,Nphases),source=0.0_pReal) allocate(lattice_NcleavageSystem(lattice_maxNcleavageFamily,Nphases),source=0) - + + do p = 1, size(config_phase) - tag = config_phase(p)%getString('lattice_structure') - select case(trim(tag(1:3))) - case('iso') - lattice_structure(p) = LATTICE_iso_ID - case('fcc') - lattice_structure(p) = LATTICE_fcc_ID - case('bcc') - lattice_structure(p) = LATTICE_bcc_ID - case('hex') - lattice_structure(p) = LATTICE_hex_ID - case('bct') - lattice_structure(p) = LATTICE_bct_ID - case('ort') - lattice_structure(p) = LATTICE_ort_ID - end select - - + lattice_C66(1,1,p) = config_phase(p)%getFloat('c11',defaultVal=0.0_pReal) lattice_C66(1,2,p) = config_phase(p)%getFloat('c12',defaultVal=0.0_pReal) lattice_C66(1,3,p) = config_phase(p)%getFloat('c13',defaultVal=0.0_pReal) @@ -581,21 +565,56 @@ subroutine lattice_init lattice_C66(4,4,p) = config_phase(p)%getFloat('c44',defaultVal=0.0_pReal) lattice_C66(5,5,p) = config_phase(p)%getFloat('c55',defaultVal=0.0_pReal) lattice_C66(6,6,p) = config_phase(p)%getFloat('c66',defaultVal=0.0_pReal) - - + CoverA = config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal) - + + tag = config_phase(p)%getString('lattice_structure') + select case(tag(1:3)) + case('iso') + lattice_structure(p) = LATTICE_iso_ID + case('fcc') + lattice_structure(p) = LATTICE_fcc_ID + case('bcc') + lattice_structure(p) = LATTICE_bcc_ID + case('hex') + if(CoverA < 1.0_pReal .or. CoverA > 2.0_pReal) call IO_error(131,el=p) + lattice_structure(p) = LATTICE_hex_ID + case('bct') + if(CoverA > 2.0_pReal) call IO_error(131,el=p) + lattice_structure(p) = LATTICE_bct_ID + case('ort') + lattice_structure(p) = LATTICE_ort_ID + case default + call IO_error(130,ext_msg='lattice_init') + end select + + lattice_C66(1:6,1:6,p) = lattice_symmetrizeC66(lattice_structure(p),lattice_C66(1:6,1:6,p)) + + lattice_mu(p) = 0.2_pReal *(lattice_C66(1,1,p) -lattice_C66(1,2,p) +3.0_pReal*lattice_C66(4,4,p)) ! (C11iso-C12iso)/2 with C11iso=(3*C11+2*C12+4*C44)/5 and C12iso=(C11+4*C12-2*C44)/5 + lattice_nu(p) = ( lattice_C66(1,1,p) +4.0_pReal*lattice_C66(1,2,p) -2.0_pReal*lattice_C66(4,4,p)) & + / (4.0_pReal*lattice_C66(1,1,p) +6.0_pReal*lattice_C66(1,2,p) +2.0_pReal*lattice_C66(4,4,p))! C12iso/(C11iso+C12iso) with C11iso=(3*C11+2*C12+4*C44)/5 and C12iso=(C11+4*C12-2*C44)/5 + + lattice_C3333(1:3,1:3,1:3,1:3,p) = math_Voigt66to3333(lattice_C66(1:6,1:6,p)) ! Literature data is Voigt + lattice_C66(1:6,1:6,p) = math_sym3333to66(lattice_C3333(1:3,1:3,1:3,1:3,p)) ! DAMASK uses Mandel-weighting + + do i = 1, 6 + if (abs(lattice_C66(i,i,p)) 2.0_pReal) & - .and. lattice_structure(p) == LATTICE_hex_ID) call IO_error(131,el=p) ! checking physical significance of c/a - if ((CoverA > 2.0_pReal) & - .and. lattice_structure(p) == LATTICE_bct_ID) call IO_error(131,el=p) ! checking physical significance of c/a + call lattice_initializeStructure(p, CoverA) enddo - + end subroutine lattice_init - - + + !-------------------------------------------------------------------------------------------------- !> @brief !!!!!!!DEPRECTATED!!!!!! !-------------------------------------------------------------------------------------------------- @@ -622,102 +637,65 @@ subroutine lattice_initializeStructure(myPhase,CoverA) integer, intent(in) :: myPhase real(pReal), intent(in) :: & CoverA - + integer :: & - i, & - myNcleavage - - lattice_C66(1:6,1:6,myPhase) = lattice_symmetrizeC66(lattice_structure(myPhase),& - lattice_C66(1:6,1:6,myPhase)) - - lattice_mu(myPhase) = 0.2_pReal *( lattice_C66(1,1,myPhase) & - - lattice_C66(1,2,myPhase) & - + 3.0_pReal*lattice_C66(4,4,myPhase)) ! (C11iso-C12iso)/2 with C11iso=(3*C11+2*C12+4*C44)/5 and C12iso=(C11+4*C12-2*C44)/5 - lattice_nu(myPhase) = ( lattice_C66(1,1,myPhase) & - + 4.0_pReal*lattice_C66(1,2,myPhase) & - - 2.0_pReal*lattice_C66(4,4,myPhase)) & - /( 4.0_pReal*lattice_C66(1,1,myPhase) & - + 6.0_pReal*lattice_C66(1,2,myPhase) & - + 2.0_pReal*lattice_C66(4,4,myPhase))! C12iso/(C11iso+C12iso) with C11iso=(3*C11+2*C12+4*C44)/5 and C12iso=(C11+4*C12-2*C44)/5 - lattice_C3333(1:3,1:3,1:3,1:3,myPhase) = math_Voigt66to3333(lattice_C66(1:6,1:6,myPhase)) ! Literature data is Voigt - lattice_C66(1:6,1:6,myPhase) = math_sym3333to66(lattice_C3333(1:3,1:3,1:3,1:3,myPhase)) ! DAMASK uses Mandel-weighting - do i = 1, 6 - if (abs(lattice_C66(i,i,myPhase)) @brief Symmetrizes stiffness matrix according to lattice type !> @details J. A. Rayne and B. S. Chandrasekhar Phys. Rev. 120, 1658 Erratum Phys. Rev. 122, 1962 !-------------------------------------------------------------------------------------------------- pure function lattice_symmetrizeC66(struct,C66) - + integer(kind(LATTICE_undefined_ID)), intent(in) :: struct real(pReal), dimension(6,6), intent(in) :: C66 real(pReal), dimension(6,6) :: lattice_symmetrizeC66 integer :: j,k - + lattice_symmetrizeC66 = 0.0_pReal - + select case(struct) case (LATTICE_iso_ID) do k=1,3 @@ -777,22 +755,22 @@ pure function lattice_symmetrizeC66(struct,C66) case default lattice_symmetrizeC66 = C66 end select - + end function lattice_symmetrizeC66 - - + + !-------------------------------------------------------------------------------------------------- !> @brief Symmetrizes 2nd order tensor according to lattice type !-------------------------------------------------------------------------------------------------- pure function lattice_symmetrize33(struct,T33) - + integer(kind(LATTICE_undefined_ID)), intent(in) :: struct real(pReal), dimension(3,3), intent(in) :: T33 real(pReal), dimension(3,3) :: lattice_symmetrize33 integer :: k - + lattice_symmetrize33 = 0.0_pReal - + select case(struct) case (LATTICE_iso_ID,LATTICE_fcc_ID,LATTICE_bcc_ID) do k=1,3 @@ -809,26 +787,26 @@ pure function lattice_symmetrize33(struct,T33) case default lattice_symmetrize33 = T33 end select - + end function lattice_symmetrize33 - - + + !-------------------------------------------------------------------------------------------------- !> @brief Characteristic shear for twinning !-------------------------------------------------------------------------------------------------- function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(characteristicShear) - + integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family character(len=*), intent(in) :: structure !< lattice structure real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), dimension(sum(Ntwin)) :: characteristicShear - + integer :: & a, & !< index of active system p, & !< index in potential system list f, & !< index of my family s !< index of my system in current family - + integer, dimension(LATTICE_HEX_NTWIN), parameter :: & HEX_SHEARTWIN = reshape( [& 1, & ! <-10.1>{10.2} @@ -856,10 +834,10 @@ function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(charact 4, & 4 & ],[LATTICE_HEX_NTWIN]) ! indicator to formulas below - + if (len_trim(structure) /= 3) & call IO_error(137,ext_msg='lattice_characteristicShear_Twin: '//trim(structure)) - + a = 0 myFamilies: do f = 1,size(Ntwin,1) mySystems: do s = 1,Ntwin(f) @@ -886,28 +864,28 @@ function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(charact end select enddo mySystems enddo myFamilies - + end function lattice_characteristicShear_Twin - - + + !-------------------------------------------------------------------------------------------------- !> @brief Rotated elasticity matrices for twinning in 66-vector notation !-------------------------------------------------------------------------------------------------- function lattice_C66_twin(Ntwin,C66,structure,CoverA) - + integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family character(len=*), intent(in) :: structure !< lattice structure real(pReal), dimension(6,6), intent(in) :: C66 !< unrotated parent stiffness matrix real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), dimension(6,6,sum(Ntwin)) :: lattice_C66_twin - + real(pReal), dimension(3,3,sum(Ntwin)):: coordinateSystem type(rotation) :: R integer :: i - + if (len_trim(structure) /= 3) & call IO_error(137,ext_msg='lattice_C66_twin: '//trim(structure)) - + select case(structure(1:3)) case('fcc') coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_FCC_NSLIPSYSTEM,LATTICE_FCC_SYSTEMTWIN,& @@ -921,35 +899,35 @@ function lattice_C66_twin(Ntwin,C66,structure,CoverA) case default call IO_error(137,ext_msg='lattice_C66_twin: '//trim(structure)) end select - + do i = 1, sum(Ntwin) call R%fromAxisAngle([coordinateSystem(1:3,2,i),PI],P=1) ! ToDo: Why always 180 deg? lattice_C66_twin(1:6,1:6,i) = R%rotTensor4sym(C66) enddo end function lattice_C66_twin - - + + !-------------------------------------------------------------------------------------------------- !> @brief Rotated elasticity matrices for transformation in 66-vector notation !-------------------------------------------------------------------------------------------------- function lattice_C66_trans(Ntrans,C_parent66,structure_target, & cOverA_trans,a_bcc,a_fcc) - + integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family character(len=*), intent(in) :: structure_target !< lattice structure real(pReal), dimension(6,6), intent(in) :: C_parent66 real(pReal), dimension(6,6,sum(Ntrans)) :: lattice_C66_trans - + real(pReal), dimension(6,6) :: C_bar66, C_target_unrotated66 real(pReal), dimension(3,3,sum(Ntrans)) :: Q,S type(rotation) :: R real(pReal) :: a_bcc, a_fcc, cOverA_trans integer :: i - + if (len_trim(structure_target) /= 3) & call IO_error(137,ext_msg='lattice_C66_trans (target): '//trim(structure_target)) - + !-------------------------------------------------------------------------------------------------- ! elasticity matrix of the target phase in cube orientation if (structure_target(1:3) == 'hex') then @@ -961,7 +939,7 @@ function lattice_C66_trans(Ntrans,C_parent66,structure_target, & C_bar66(1,3) = (C_parent66(1,1) + 2.0_pReal*C_parent66(1,2) - 2.0_pReal*C_parent66(4,4))/3.0_pReal C_bar66(4,4) = (C_parent66(1,1) - C_parent66(1,2) + C_parent66(4,4))/3.0_pReal C_bar66(1,4) = (C_parent66(1,1) - C_parent66(1,2) - 2.0_pReal*C_parent66(4,4)) /(3.0_pReal*sqrt(2.0_pReal)) - + C_target_unrotated66 = 0.0_pReal C_target_unrotated66(1,1) = C_bar66(1,1) - C_bar66(1,4)**2.0_pReal/C_bar66(4,4) C_target_unrotated66(1,2) = C_bar66(1,2) + C_bar66(1,4)**2.0_pReal/C_bar66(4,4) @@ -976,22 +954,22 @@ function lattice_C66_trans(Ntrans,C_parent66,structure_target, & else call IO_error(137,ext_msg='lattice_C66_trans : '//trim(structure_target)) endif - + do i = 1, 6 if (abs(C_target_unrotated66(i,i)) @brief Non-schmid projections for bcc with up to 6 coefficients ! Koester et al. 2012, Acta Materialia 60 (2012) 3894–3901, eq. (17) @@ -1003,19 +981,19 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc real(pReal), dimension(:), intent(in) :: nonSchmidCoefficients !< non-Schmid coefficients for projections integer, intent(in) :: sense !< sense (-1,+1) real(pReal), dimension(1:3,1:3,sum(Nslip)) :: nonSchmidMatrix - + real(pReal), dimension(1:3,1:3,sum(Nslip)) :: coordinateSystem !< coordinate system of slip system real(pReal), dimension(3) :: direction, normal, np type(rotation) :: R integer :: i - + if (abs(sense) /= 1) call IO_error(0,ext_msg='lattice_nonSchmidMatrix') - + coordinateSystem = buildCoordinateSystem(Nslip,LATTICE_BCC_NSLIPSYSTEM,LATTICE_BCC_SYSTEMSLIP,& 'bcc',0.0_pReal) coordinateSystem(1:3,1,1:sum(Nslip)) = coordinateSystem(1:3,1,1:sum(Nslip)) *real(sense,pReal) ! convert unidirectional coordinate system nonSchmidMatrix = lattice_SchmidMatrix_slip(Nslip,'bcc',0.0_pReal) ! Schmid contribution - + do i = 1,sum(Nslip) direction = coordinateSystem(1:3,1,i) normal = coordinateSystem(1:3,2,i) @@ -1038,8 +1016,8 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc enddo end function lattice_nonSchmidMatrix - - + + !-------------------------------------------------------------------------------------------------- !> @brief Slip-slip interaction matrix !> details only active slip systems are considered @@ -1050,10 +1028,10 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul real(pReal), dimension(:), intent(in) :: interactionValues !< values for slip-slip interaction character(len=*), intent(in) :: structure !< lattice structure real(pReal), dimension(sum(Nslip),sum(Nslip)) :: interactionMatrix - + integer, dimension(:), allocatable :: NslipMax integer, dimension(:,:), allocatable :: interactionTypes - + integer, dimension(LATTICE_FCC_NSLIP,LATTICE_FCC_NSLIP), parameter :: & FCC_INTERACTIONSLIPSLIP = reshape( [& 1, 2, 2, 4, 6, 5, 3, 5, 5, 4, 5, 6, 9,10, 9,10,11,12, & ! -----> acting @@ -1068,7 +1046,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul 4, 5, 6, 3, 5, 5, 4, 6, 5, 1, 2, 2, 10, 9, 9,10,12,11, & 5, 3, 5, 5, 4, 6, 6, 4, 5, 2, 1, 2, 10, 9,11,12,10, 9, & 6, 5, 4, 5, 6, 4, 5, 5, 3, 2, 2, 1, 12,11, 9,10,10, 9, & - + 9, 9,11, 9, 9,11,10,10,12,10,10,12, 1, 7, 8, 8, 8, 8, & 10,10,12,10,10,12, 9, 9,11, 9, 9,11, 7, 1, 8, 8, 8, 8, & 9,11, 9,10,12,10,10,12,10, 9,11, 9, 8, 8, 1, 7, 8, 8, & @@ -1088,7 +1066,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul !<10: similar to glissile junctions in <110>{111} btw one {110} and one {111} plane !<11: crossing btw one {110} and one {111} plane !<12: collinear btw one {110} and one {111} plane - + integer, dimension(LATTICE_BCC_NSLIP,LATTICE_BCC_NSLIP), parameter :: & BCC_INTERACTIONSLIPSLIP = reshape( [& 1,2,6,6,5,4,4,3,4,3,5,4, 6,6,4,3,3,4,6,6,4,3,6,6, & ! -----> acting @@ -1123,7 +1101,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul !< 4: mixed-asymmetrical junction !< 5: mixed-symmetrical junction !< 6: edge junction - + integer, dimension(LATTICE_HEX_NSLIP,LATTICE_HEX_NSLIP), parameter :: & HEX_INTERACTIONSLIPSLIP = reshape( [& 1, 2, 2, 3, 3, 3, 7, 7, 7, 13,13,13,13,13,13, 21,21,21,21,21,21,21,21,21,21,21,21, 31,31,31,31,31,31, & ! -----> acting @@ -1165,7 +1143,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul 42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 37,37,37,37,36,37, & 42,42,42, 41,41,41, 40,40,40, 39,39,39,39,39,39, 38,38,38,38,38,38,38,38,38,38,38,38, 37,37,37,37,37,36 & ],shape(HEX_INTERACTIONSLIPSLIP)) !< Slip--slip interaction types for hex (onion peel naming scheme) - + integer, dimension(LATTICE_BCT_NSLIP,LATTICE_BCT_NSLIP), parameter :: & BCT_INTERACTIONSLIPSLIP = reshape( [& 1, 2, 3, 3, 7, 7, 13, 13, 13, 13, 21, 21, 31, 31, 31, 31, 43, 43, 57, 57, 73, 73, 73, 73, 91, 91, 91, 91, 91, 91, 91, 91, 111, 111, 111, 111, 133,133,133,133,133,133,133,133, 157,157,157,157,157,157,157,157, & ! -----> acting @@ -1233,11 +1211,11 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul 182,182, 181,181, 180,180, 179,179,179,179, 178,178, 177,177,177,177, 176,176, 175,175, 174,174,174,174, 173,173,173,173,173,173,173,173, 172, 172, 172, 172, 171,171,171,171,171,171,171,171, 169,170,170,170,170,170,169,170, & 182,182, 181,181, 180,180, 179,179,179,179, 178,178, 177,177,177,177, 176,176, 175,175, 174,174,174,174, 173,173,173,173,173,173,173,173, 172, 172, 172, 172, 171,171,171,171,171,171,171,171, 169,170,170,170,170,170,170,169 & ],shape(BCT_INTERACTIONSLIPSLIP)) - - + + if (len_trim(structure) /= 3) & call IO_error(137,ext_msg='lattice_interaction_SlipBySlip: '//trim(structure)) - + select case(structure(1:3)) case('fcc') interactionTypes = FCC_INTERACTIONSLIPSLIP @@ -1254,26 +1232,26 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul case default call IO_error(137,ext_msg='lattice_interaction_SlipBySlip: '//trim(structure)) end select - + interactionMatrix = buildInteraction(Nslip,Nslip,NslipMax,NslipMax,interactionValues,interactionTypes) - + end function lattice_interaction_SlipBySlip - - + + !-------------------------------------------------------------------------------------------------- !> @brief Twin-twin interaction matrix !> details only active twin systems are considered !-------------------------------------------------------------------------------------------------- function lattice_interaction_TwinByTwin(Ntwin,interactionValues,structure) result(interactionMatrix) - + integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family real(pReal), dimension(:), intent(in) :: interactionValues !< values for twin-twin interaction character(len=*), intent(in) :: structure !< lattice structure real(pReal), dimension(sum(Ntwin),sum(Ntwin)) :: interactionMatrix - + integer, dimension(:), allocatable :: NtwinMax integer, dimension(:,:), allocatable :: interactionTypes - + integer, dimension(LATTICE_FCC_NTWIN,LATTICE_FCC_NTWIN), parameter :: & FCC_INTERACTIONTWINTWIN = reshape( [& 1,1,1,2,2,2,2,2,2,2,2,2, & ! -----> acting @@ -1289,7 +1267,7 @@ function lattice_interaction_TwinByTwin(Ntwin,interactionValues,structure) resul 2,2,2,2,2,2,2,2,2,1,1,1, & 2,2,2,2,2,2,2,2,2,1,1,1 & ],shape(FCC_INTERACTIONTWINTWIN)) !< Twin-twin interaction types for fcc - + integer, dimension(LATTICE_BCC_NTWIN,LATTICE_BCC_NTWIN), parameter :: & BCC_INTERACTIONTWINTWIN = reshape( [& 1,3,3,3,3,3,3,2,3,3,2,3, & ! -----> acting @@ -1338,10 +1316,10 @@ function lattice_interaction_TwinByTwin(Ntwin,interactionValues,structure) resul 20,20,20,20,20,20, 19,19,19,19,19,19, 18,18,18,18,18,18, 17,17,17,17,16,17, & 20,20,20,20,20,20, 19,19,19,19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,16 & ],shape(HEX_INTERACTIONTWINTWIN)) !< Twin-twin interaction types for hex - + if (len_trim(structure) /= 3) & call IO_error(137,ext_msg='lattice_interaction_TwinByTwin: '//trim(structure)) - + select case(structure(1:3)) case('fcc') interactionTypes = FCC_INTERACTIONTWINTWIN @@ -1355,26 +1333,26 @@ function lattice_interaction_TwinByTwin(Ntwin,interactionValues,structure) resul case default call IO_error(137,ext_msg='lattice_interaction_TwinByTwin: '//trim(structure)) end select - + interactionMatrix = buildInteraction(Ntwin,Ntwin,NtwinMax,NtwinMax,interactionValues,interactionTypes) - + end function lattice_interaction_TwinByTwin - - + + !-------------------------------------------------------------------------------------------------- !> @brief Trans-trans interaction matrix !> details only active trans systems are considered !-------------------------------------------------------------------------------------------------- function lattice_interaction_TransByTrans(Ntrans,interactionValues,structure) result(interactionMatrix) - + integer, dimension(:), intent(in) :: Ntrans !< number of active trans systems per family real(pReal), dimension(:), intent(in) :: interactionValues !< values for trans-trans interaction character(len=*), intent(in) :: structure !< lattice structure (parent crystal) real(pReal), dimension(sum(Ntrans),sum(Ntrans)) :: interactionMatrix - + integer, dimension(:), allocatable :: NtransMax integer, dimension(:,:), allocatable :: interactionTypes - + integer, dimension(LATTICE_FCC_NTRANS,LATTICE_FCC_NTRANS), parameter :: & FCC_INTERACTIONTRANSTRANS = reshape( [& 1,1,1,2,2,2,2,2,2,2,2,2, & ! -----> acting @@ -1390,44 +1368,44 @@ function lattice_interaction_TransByTrans(Ntrans,interactionValues,structure) re 2,2,2,2,2,2,2,2,2,1,1,1, & 2,2,2,2,2,2,2,2,2,1,1,1 & ],shape(FCC_INTERACTIONTRANSTRANS)) !< Trans-trans interaction types for fcc - + if (len_trim(structure) /= 3) & call IO_error(137,ext_msg='lattice_interaction_TransByTrans: '//trim(structure)) - + if(structure(1:3) == 'fcc') then interactionTypes = FCC_INTERACTIONTRANSTRANS NtransMax = LATTICE_FCC_NTRANSSYSTEM else call IO_error(137,ext_msg='lattice_interaction_TransByTrans: '//trim(structure)) end if - + interactionMatrix = buildInteraction(Ntrans,Ntrans,NtransMax,NtransMax,interactionValues,interactionTypes) - + end function lattice_interaction_TransByTrans - - + + !-------------------------------------------------------------------------------------------------- !> @brief Slip-twin interaction matrix !> details only active slip and twin systems are considered !-------------------------------------------------------------------------------------------------- function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,structure) result(interactionMatrix) - + integer, dimension(:), intent(in) :: Nslip, & !< number of active slip systems per family Ntwin !< number of active twin systems per family real(pReal), dimension(:), intent(in) :: interactionValues !< values for slip-twin interaction character(len=*), intent(in) :: structure !< lattice structure real(pReal), dimension(sum(Nslip),sum(Ntwin)) :: interactionMatrix - + integer, dimension(:), allocatable :: NslipMax, & NtwinMax integer, dimension(:,:), allocatable :: interactionTypes - + integer, dimension(LATTICE_FCC_NTWIN,LATTICE_FCC_NSLIP), parameter :: & FCC_INTERACTIONSLIPTWIN = reshape( [& 1,1,1,3,3,3,2,2,2,3,3,3, & ! -----> twin (acting) 1,1,1,3,3,3,3,3,3,2,2,2, & ! | 1,1,1,2,2,2,3,3,3,3,3,3, & ! | - 3,3,3,1,1,1,3,3,3,2,2,2, & ! v + 3,3,3,1,1,1,3,3,3,2,2,2, & ! v 3,3,3,1,1,1,2,2,2,3,3,3, & ! slip (reacting) 2,2,2,1,1,1,3,3,3,3,3,3, & 2,2,2,3,3,3,1,1,1,3,3,3, & @@ -1436,7 +1414,7 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,structure) 3,3,3,2,2,2,3,3,3,1,1,1, & 2,2,2,3,3,3,3,3,3,1,1,1, & 3,3,3,3,3,3,2,2,2,1,1,1, & - + 4,4,4,4,4,4,4,4,4,4,4,4, & 4,4,4,4,4,4,4,4,4,4,4,4, & 4,4,4,4,4,4,4,4,4,4,4,4, & @@ -1520,10 +1498,10 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,structure) 21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24 & ! ],shape(HEX_INTERACTIONSLIPTWIN)) !< Slip-twin interaction types for hex - + if (len_trim(structure) /= 3) & call IO_error(137,ext_msg='lattice_interaction_SlipByTwin: '//trim(structure)) - + select case(structure(1:3)) case('fcc') interactionTypes = FCC_INTERACTIONSLIPTWIN @@ -1540,28 +1518,28 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,structure) case default call IO_error(137,ext_msg='lattice_interaction_SlipByTwin: '//trim(structure)) end select - + interactionMatrix = buildInteraction(Nslip,Ntwin,NslipMax,NtwinMax,interactionValues,interactionTypes) - + end function lattice_interaction_SlipByTwin - - + + !-------------------------------------------------------------------------------------------------- !> @brief Slip-trans interaction matrix !> details only active slip and trans systems are considered !-------------------------------------------------------------------------------------------------- function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,structure) result(interactionMatrix) - + integer, dimension(:), intent(in) :: Nslip, & !< number of active slip systems per family Ntrans !< number of active trans systems per family real(pReal), dimension(:), intent(in) :: interactionValues !< values for slip-trans interaction character(len=*), intent(in) :: structure !< lattice structure (parent crystal) real(pReal), dimension(sum(Nslip),sum(Ntrans)) :: interactionMatrix - + integer, dimension(:), allocatable :: NslipMax, & NtransMax integer, dimension(:,:), allocatable :: interactionTypes - + integer, dimension(LATTICE_FCC_NTRANS,LATTICE_FCC_NSLIP), parameter :: & FCC_INTERACTIONSLIPTRANS = reshape( [& 1,1,1,3,3,3,2,2,2,3,3,3, & ! -----> trans (acting) @@ -1576,7 +1554,7 @@ function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,structur 3,3,3,2,2,2,3,3,3,1,1,1, & 2,2,2,3,3,3,3,3,3,1,1,1, & 3,3,3,3,3,3,2,2,2,1,1,1, & - + 4,4,4,4,4,4,4,4,4,4,4,4, & 4,4,4,4,4,4,4,4,4,4,4,4, & 4,4,4,4,4,4,4,4,4,4,4,4, & @@ -1584,10 +1562,10 @@ function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,structur 4,4,4,4,4,4,4,4,4,4,4,4, & 4,4,4,4,4,4,4,4,4,4,4,4 & ],shape(FCC_INTERACTIONSLIPTRANS)) !< Slip-trans interaction types for fcc - + if (len_trim(structure) /= 3) & call IO_error(137,ext_msg='lattice_interaction_SlipByTrans: '//trim(structure)) - + select case(structure(1:3)) case('fcc') interactionTypes = FCC_INTERACTIONSLIPTRANS @@ -1596,34 +1574,34 @@ function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,structur case default call IO_error(137,ext_msg='lattice_interaction_SlipByTrans: '//trim(structure)) end select - + interactionMatrix = buildInteraction(Nslip,Ntrans,NslipMax,NtransMax,interactionValues,interactionTypes) - + end function lattice_interaction_SlipByTrans - - + + !-------------------------------------------------------------------------------------------------- !> @brief Twin-slip interaction matrix !> details only active twin and slip systems are considered !-------------------------------------------------------------------------------------------------- function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,structure) result(interactionMatrix) - + integer, dimension(:), intent(in) :: Ntwin, & !< number of active twin systems per family Nslip !< number of active slip systems per family real(pReal), dimension(:), intent(in) :: interactionValues !< values for twin-twin interaction character(len=*), intent(in) :: structure !< lattice structure real(pReal), dimension(sum(Ntwin),sum(Nslip)) :: interactionMatrix - + integer, dimension(:), allocatable :: NtwinMax, & NslipMax integer, dimension(:,:), allocatable :: interactionTypes - + integer, dimension(LATTICE_FCC_NSLIP,LATTICE_FCC_NTWIN), parameter :: & FCC_INTERACTIONTWINSLIP = 1 !< Twin-slip interaction types for fcc - + integer, dimension(LATTICE_BCC_NSLIP,LATTICE_BCC_NTWIN), parameter :: & BCC_INTERACTIONTWINSLIP = 1 !< Twin-slip interaction types for bcc - + integer, dimension(LATTICE_HEX_NSLIP,LATTICE_HEX_NTWIN), parameter :: & HEX_INTERACTIONTWINSLIP = reshape( [& 1, 1, 1, 5, 5, 5, 9, 9, 9, 13,13,13,13,13,13, 17,17,17,17,17,17,17,17,17,17,17,17, 21,21,21,21,21,21, & ! ----> slip (acting) @@ -1654,10 +1632,10 @@ function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,structure) 4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24, & 4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24 & ],shape(HEX_INTERACTIONTWINSLIP)) !< Twin-slip interaction types for hex - + if (len_trim(structure) /= 3) & call IO_error(137,ext_msg='lattice_interaction_TwinBySlip: '//trim(structure)) - + select case(structure(1:3)) case('fcc') interactionTypes = FCC_INTERACTIONTWINSLIP @@ -1674,31 +1652,31 @@ function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,structure) case default call IO_error(137,ext_msg='lattice_interaction_TwinBySlip: '//trim(structure)) end select - + interactionMatrix = buildInteraction(Ntwin,Nslip,NtwinMax,NslipMax,interactionValues,interactionTypes) - + end function lattice_interaction_TwinBySlip - - + + !-------------------------------------------------------------------------------------------------- !> @brief Schmid matrix for slip !> details only active slip systems are considered !-------------------------------------------------------------------------------------------------- function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) result(SchmidMatrix) - + integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family character(len=*), intent(in) :: structure !< lattice structure real(pReal), intent(in) :: cOverA real(pReal), dimension(3,3,sum(Nslip)) :: SchmidMatrix - + real(pReal), dimension(3,3,sum(Nslip)) :: coordinateSystem real(pReal), dimension(:,:), allocatable :: slipSystems integer, dimension(:), allocatable :: NslipMax integer :: i - + if (len_trim(structure) /= 3) & call IO_error(137,ext_msg='lattice_SchmidMatrix_slip: '//trim(structure)) - + select case(structure(1:3)) case('fcc') NslipMax = LATTICE_FCC_NSLIPSYSTEM @@ -1715,42 +1693,42 @@ function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) result(SchmidMatrix) case default call IO_error(137,ext_msg='lattice_SchmidMatrix_slip: '//trim(structure)) end select - + if (any(NslipMax(1:size(Nslip)) - Nslip < 0)) & call IO_error(145,ext_msg='Nslip '//trim(structure)) if (any(Nslip < 0)) & call IO_error(144,ext_msg='Nslip '//trim(structure)) - + coordinateSystem = buildCoordinateSystem(Nslip,NslipMax,slipSystems,structure,cOverA) - + do i = 1, sum(Nslip) SchmidMatrix(1:3,1:3,i) = math_outer(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i)) if (abs(math_trace33(SchmidMatrix(1:3,1:3,i))) > tol_math_check) & call IO_error(0,i,ext_msg = 'dilatational Schmid matrix for slip') enddo - + end function lattice_SchmidMatrix_slip - - + + !-------------------------------------------------------------------------------------------------- !> @brief Schmid matrix for twinning !> details only active twin systems are considered !-------------------------------------------------------------------------------------------------- function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) result(SchmidMatrix) - + integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family character(len=*), intent(in) :: structure !< lattice structure real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), dimension(3,3,sum(Ntwin)) :: SchmidMatrix - + real(pReal), dimension(3,3,sum(Ntwin)) :: coordinateSystem real(pReal), dimension(:,:), allocatable :: twinSystems integer, dimension(:), allocatable :: NtwinMax integer :: i - + if (len_trim(structure) /= 3) & call IO_error(137,ext_msg='lattice_SchmidMatrix_twin: '//trim(structure)) - + select case(structure(1:3)) case('fcc') NtwinMax = LATTICE_FCC_NTWINSYSTEM @@ -1764,37 +1742,37 @@ function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) result(SchmidMatrix) case default call IO_error(137,ext_msg='lattice_SchmidMatrix_twin: '//trim(structure)) end select - + if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0)) & call IO_error(145,ext_msg='Ntwin '//trim(structure)) if (any(Ntwin < 0)) & call IO_error(144,ext_msg='Ntwin '//trim(structure)) - + coordinateSystem = buildCoordinateSystem(Ntwin,NtwinMax,twinSystems,structure,cOverA) - + do i = 1, sum(Ntwin) SchmidMatrix(1:3,1:3,i) = math_outer(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i)) if (abs(math_trace33(SchmidMatrix(1:3,1:3,i))) > tol_math_check) & call IO_error(0,i,ext_msg = 'dilatational Schmid matrix for twin') enddo - + end function lattice_SchmidMatrix_twin - - + + !-------------------------------------------------------------------------------------------------- !> @brief Schmid matrix for twinning !> details only active twin systems are considered !-------------------------------------------------------------------------------------------------- function lattice_SchmidMatrix_trans(Ntrans,structure_target,cOverA,a_bcc,a_fcc) result(SchmidMatrix) - + integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family character(len=*), intent(in) :: structure_target !< lattice structure real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), dimension(3,3,sum(Ntrans)) :: SchmidMatrix - + real(pReal), dimension(3,3,sum(Ntrans)) :: devNull real(pReal) :: a_bcc, a_fcc - + if (len_trim(structure_target) /= 3) & call IO_error(137,ext_msg='lattice_SchmidMatrix_trans: '//trim(structure_target)) if (structure_target(1:3) /= 'bcc' .and. structure_target(1:3) /= 'hex') & @@ -1802,15 +1780,15 @@ function lattice_SchmidMatrix_trans(Ntrans,structure_target,cOverA,a_bcc,a_fcc) if (structure_target(1:3) == 'hex' .and. (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal)) & call IO_error(131,ext_msg='lattice_SchmidMatrix_trans: '//trim(structure_target)) - + if (structure_target(1:3) == 'bcc' .and. (a_bcc <= 0.0_pReal .or. a_fcc <= 0.0_pReal)) & call IO_error(134,ext_msg='lattice_SchmidMatrix_trans: '//trim(structure_target)) - + call buildTransformationSystem(devNull,SchmidMatrix,Ntrans,cOverA,a_fcc,a_bcc) - + end function lattice_SchmidMatrix_trans - - + + !-------------------------------------------------------------------------------------------------- !> @brief Schmid matrix for cleavage !> details only active cleavage systems are considered @@ -1821,15 +1799,15 @@ function lattice_SchmidMatrix_cleavage(Ncleavage,structure,cOverA) result(Schmid character(len=*), intent(in) :: structure !< lattice structure real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), dimension(3,3,3,sum(Ncleavage)) :: SchmidMatrix - + real(pReal), dimension(3,3,sum(Ncleavage)) :: coordinateSystem real(pReal), dimension(:,:), allocatable :: cleavageSystems integer, dimension(:), allocatable :: NcleavageMax integer :: i - + if (len_trim(structure) /= 3) & call IO_error(137,ext_msg='lattice_SchmidMatrix_cleavage: '//trim(structure)) - + select case(structure(1:3)) case('iso') NcleavageMax = LATTICE_ISO_NCLEAVAGESYSTEM @@ -1849,56 +1827,56 @@ function lattice_SchmidMatrix_cleavage(Ncleavage,structure,cOverA) result(Schmid case default call IO_error(137,ext_msg='lattice_SchmidMatrix_cleavage: '//trim(structure)) end select - + if (any(NcleavageMax(1:size(Ncleavage)) - Ncleavage < 0)) & call IO_error(145,ext_msg='Ncleavage '//trim(structure)) if (any(Ncleavage < 0)) & call IO_error(144,ext_msg='Ncleavage '//trim(structure)) - + coordinateSystem = buildCoordinateSystem(Ncleavage,NcleavageMax,cleavageSystems,structure,cOverA) - + do i = 1, sum(Ncleavage) SchmidMatrix(1:3,1:3,1,i) = math_outer(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i)) SchmidMatrix(1:3,1:3,2,i) = math_outer(coordinateSystem(1:3,3,i),coordinateSystem(1:3,2,i)) SchmidMatrix(1:3,1:3,3,i) = math_outer(coordinateSystem(1:3,2,i),coordinateSystem(1:3,2,i)) enddo - + end function lattice_SchmidMatrix_cleavage - - + + !-------------------------------------------------------------------------------------------------- !> @brief Slip direction of slip systems (|| b) !-------------------------------------------------------------------------------------------------- function lattice_slip_direction(Nslip,structure,cOverA) result(d) - + integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family character(len=*), intent(in) :: structure !< lattice structure real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), dimension(3,sum(Nslip)) :: d - + real(pReal), dimension(3,3,sum(Nslip)) :: coordinateSystem - + coordinateSystem = coordinateSystem_slip(Nslip,structure,cOverA) d = coordinateSystem(1:3,1,1:sum(Nslip)) - + end function lattice_slip_direction - - + + !-------------------------------------------------------------------------------------------------- !> @brief Normal direction of slip systems (|| n) !-------------------------------------------------------------------------------------------------- function lattice_slip_normal(Nslip,structure,cOverA) result(n) - + integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family character(len=*), intent(in) :: structure !< lattice structure real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), dimension(3,sum(Nslip)) :: n - + real(pReal), dimension(3,3,sum(Nslip)) :: coordinateSystem - + coordinateSystem = coordinateSystem_slip(Nslip,structure,cOverA) n = coordinateSystem(1:3,2,1:sum(Nslip)) - + end function lattice_slip_normal @@ -1906,41 +1884,41 @@ end function lattice_slip_normal !> @brief Transverse direction of slip systems ( || t = b x n) !-------------------------------------------------------------------------------------------------- function lattice_slip_transverse(Nslip,structure,cOverA) result(t) - + integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family character(len=*), intent(in) :: structure !< lattice structure real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), dimension(3,sum(Nslip)) :: t - + real(pReal), dimension(3,3,sum(Nslip)) :: coordinateSystem - + coordinateSystem = coordinateSystem_slip(Nslip,structure,cOverA) t = coordinateSystem(1:3,3,1:sum(Nslip)) - + end function lattice_slip_transverse - - + + !-------------------------------------------------------------------------------------------------- !> @brief Projection of the transverse direction onto the slip plane !> @details: This projection is used to calculate forest hardening for edge dislocations !-------------------------------------------------------------------------------------------------- function slipProjection_transverse(Nslip,structure,cOverA) result(projection) - + integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family character(len=*), intent(in) :: structure !< lattice structure real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), dimension(sum(Nslip),sum(Nslip)) :: projection - + real(pReal), dimension(3,sum(Nslip)) :: n, t integer :: i, j - + n = lattice_slip_normal (Nslip,structure,cOverA) t = lattice_slip_transverse(Nslip,structure,cOverA) - + do i=1, sum(Nslip); do j=1, sum(Nslip) projection(i,j) = abs(math_inner(n(:,i),t(:,j))) enddo; enddo - + end function slipProjection_transverse @@ -1952,15 +1930,15 @@ function lattice_labels_slip(Nslip,structure) result(labels) integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family character(len=*), intent(in) :: structure !< lattice structure - + character(len=:), dimension(:), allocatable :: labels real(pReal), dimension(:,:), allocatable :: slipSystems integer, dimension(:), allocatable :: NslipMax - + if (len_trim(structure) /= 3) & call IO_error(137,ext_msg='lattice_labels_slip: '//trim(structure)) - + select case(structure(1:3)) case('fcc') NslipMax = LATTICE_FCC_NSLIPSYSTEM @@ -1977,14 +1955,14 @@ function lattice_labels_slip(Nslip,structure) result(labels) case default call IO_error(137,ext_msg='lattice_labels_slip: '//trim(structure)) end select - + if (any(NslipMax(1:size(Nslip)) - Nslip < 0)) & call IO_error(145,ext_msg='Nslip '//trim(structure)) if (any(Nslip < 0)) & call IO_error(144,ext_msg='Nslip '//trim(structure)) - + labels = getLabels(Nslip,NslipMax,slipSystems) - + end function lattice_labels_slip @@ -1996,15 +1974,15 @@ function lattice_labels_twin(Ntwin,structure) result(labels) integer, dimension(:), intent(in) :: Ntwin !< number of active slip systems per family character(len=*), intent(in) :: structure !< lattice structure - + character(len=:), dimension(:), allocatable :: labels real(pReal), dimension(:,:), allocatable :: twinSystems integer, dimension(:), allocatable :: NtwinMax - + if (len_trim(structure) /= 3) & call IO_error(137,ext_msg='lattice_labels_twin: '//trim(structure)) - + select case(structure(1:3)) case('fcc') NtwinMax = LATTICE_FCC_NTWINSYSTEM @@ -2018,14 +1996,14 @@ function lattice_labels_twin(Ntwin,structure) result(labels) case default call IO_error(137,ext_msg='lattice_labels_twin: '//trim(structure)) end select - + if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0)) & call IO_error(145,ext_msg='Ntwin '//trim(structure)) if (any(Ntwin < 0)) & call IO_error(144,ext_msg='Ntwin '//trim(structure)) - + labels = getLabels(Ntwin,NtwinMax,twinSystems) - + end function lattice_labels_twin @@ -2034,25 +2012,25 @@ end function lattice_labels_twin !> @details: This projection is used to calculate forest hardening for screw dislocations !-------------------------------------------------------------------------------------------------- function slipProjection_direction(Nslip,structure,cOverA) result(projection) - + integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family character(len=*), intent(in) :: structure !< lattice structure real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), dimension(sum(Nslip),sum(Nslip)) :: projection - + real(pReal), dimension(3,sum(Nslip)) :: n, d integer :: i, j - + n = lattice_slip_normal (Nslip,structure,cOverA) d = lattice_slip_direction(Nslip,structure,cOverA) - + do i=1, sum(Nslip); do j=1, sum(Nslip) projection(i,j) = abs(math_inner(n(:,i),d(:,j))) enddo; enddo - + end function slipProjection_direction - - + + !-------------------------------------------------------------------------------------------------- !> @brief build a local coordinate system on slip systems !> @details Order: Direction, plane (normal), and common perpendicular @@ -2063,13 +2041,13 @@ function coordinateSystem_slip(Nslip,structure,cOverA) result(coordinateSystem) character(len=*), intent(in) :: structure !< lattice structure real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), dimension(3,3,sum(Nslip)) :: coordinateSystem - + real(pReal), dimension(:,:), allocatable :: slipSystems integer, dimension(:), allocatable :: NslipMax - + if (len_trim(structure) /= 3) & call IO_error(137,ext_msg='coordinateSystem_slip: '//trim(structure)) - + select case(structure(1:3)) case('fcc') NslipMax = LATTICE_FCC_NSLIPSYSTEM @@ -2086,22 +2064,22 @@ function coordinateSystem_slip(Nslip,structure,cOverA) result(coordinateSystem) case default call IO_error(137,ext_msg='coordinateSystem_slip: '//trim(structure)) end select - + if (any(NslipMax(1:size(Nslip)) - Nslip < 0)) & call IO_error(145,ext_msg='Nslip '//trim(structure)) if (any(Nslip < 0)) & call IO_error(144,ext_msg='Nslip '//trim(structure)) - + coordinateSystem = buildCoordinateSystem(Nslip,NslipMax,slipSystems,structure,cOverA) - + end function coordinateSystem_slip - - + + !-------------------------------------------------------------------------------------------------- !> @brief Populates reduced interaction matrix !-------------------------------------------------------------------------------------------------- function buildInteraction(reacting_used,acting_used,reacting_max,acting_max,values,matrix) - + integer, dimension(:), intent(in) :: & reacting_used, & !< # of reacting systems per family as specified in material.config acting_used, & !< # of acting systems per family as specified in material.config @@ -2110,42 +2088,42 @@ function buildInteraction(reacting_used,acting_used,reacting_max,acting_max,valu real(pReal), dimension(:), intent(in) :: values !< interaction values integer, dimension(:,:), intent(in) :: matrix !< interaction types real(pReal), dimension(sum(reacting_used),sum(acting_used)) :: buildInteraction - + integer :: & acting_family_index, acting_family, acting_system, & reacting_family_index, reacting_family, reacting_system, & i,j,k,l - + do acting_family = 1,size(acting_used,1) acting_family_index = sum(acting_used(1:acting_family-1)) do acting_system = 1,acting_used(acting_family) - + do reacting_family = 1,size(reacting_used,1) reacting_family_index = sum(reacting_used(1:reacting_family-1)) do reacting_system = 1,reacting_used(reacting_family) - + i = sum( acting_max(1: acting_family-1)) + acting_system j = sum(reacting_max(1:reacting_family-1)) + reacting_system - + k = acting_family_index + acting_system l = reacting_family_index + reacting_system - + if (matrix(i,j) > size(values)) call IO_error(138,ext_msg='buildInteraction') - + buildInteraction(l,k) = values(matrix(i,j)) - + enddo; enddo enddo; enddo - + end function buildInteraction - - + + !-------------------------------------------------------------------------------------------------- !> @brief build a local coordinate system on slip, twin, trans, cleavage systems !> @details Order: Direction, plane (normal), and common perpendicular !-------------------------------------------------------------------------------------------------- function buildCoordinateSystem(active,potential,system,structure,cOverA) - + integer, dimension(:), intent(in) :: & active, & !< # of active systems per family potential !< # of potential systems per family @@ -2157,7 +2135,7 @@ function buildCoordinateSystem(active,potential,system,structure,cOverA) cOverA real(pReal), dimension(3,3,sum(active)) :: & buildCoordinateSystem - + real(pReal), dimension(3) :: & direction, normal integer :: & @@ -2165,26 +2143,26 @@ function buildCoordinateSystem(active,potential,system,structure,cOverA) p, & !< index in potential system matrix f, & !< index of my family s !< index of my system in current family - + if (len_trim(structure) /= 3) & call IO_error(137,ext_msg='buildCoordinateSystem: '//trim(structure)) if (trim(structure(1:3)) == 'bct' .and. cOverA > 2.0_pReal) & call IO_error(131,ext_msg='buildCoordinateSystem:'//trim(structure)) if (trim(structure(1:3)) == 'hex' .and. (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal)) & call IO_error(131,ext_msg='buildCoordinateSystem:'//trim(structure)) - + a = 0 activeFamilies: do f = 1,size(active,1) activeSystems: do s = 1,active(f) a = a + 1 p = sum(potential(1:f-1))+s - + select case(trim(structure(1:3))) - + case ('fcc','bcc','iso','ort','bct') direction = system(1:3,p) normal = system(4:6,p) - + case ('hex') direction = [ system(1,p)*1.5_pReal, & (system(1,p)+2.0_pReal*system(2,p))*sqrt(0.75_pReal), & @@ -2192,23 +2170,23 @@ function buildCoordinateSystem(active,potential,system,structure,cOverA) normal = [ system(5,p), & (system(5,p)+2.0_pReal*system(6,p))/sqrt(3.0_pReal), & system(8,p)/cOverA ] ! plane (hkil)->(h (h+2k)/sqrt(3) l/(p/a)) - + case default call IO_error(137,ext_msg='buildCoordinateSystem: '//trim(structure)) - + end select - + buildCoordinateSystem(1:3,1,a) = direction/norm2(direction) buildCoordinateSystem(1:3,2,a) = normal /norm2(normal) buildCoordinateSystem(1:3,3,a) = math_cross(direction/norm2(direction),& normal /norm2(normal)) - + enddo activeSystems enddo activeFamilies - + end function buildCoordinateSystem - - + + !-------------------------------------------------------------------------------------------------- !> @brief Helper function to define transformation systems ! Needed to calculate Schmid matrix and rotated stiffness matrices. @@ -2216,7 +2194,7 @@ end function buildCoordinateSystem ! set a_Xcc = 0.0 for fcc -> hex transformation !-------------------------------------------------------------------------------------------------- subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc) - + integer, dimension(:), intent(in) :: & Ntrans real(pReal), dimension(3,3,sum(Ntrans)), intent(out) :: & @@ -2226,10 +2204,10 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc) cOverA, & !< c/a for target hex structure a_bcc, & !< lattice parameter a for target bcc structure a_fcc !< lattice parameter a for parent fcc structure - + type(rotation) :: & R, & !< Pitsch rotation - B !< Rotation of fcc to Bain coordinate system + B !< Rotation of fcc to Bain coordinate system real(pReal), dimension(3,3) :: & U, & !< Bain deformation ss, sd @@ -2267,7 +2245,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc) 0.0, 1.0, 0.0, 10.26, & 0.0,-1.0, 0.0, 10.26 & ],shape(LATTICE_FCCTOBCC_SYSTEMTRANS)) - + integer, dimension(9,LATTICE_fcc_Ntrans), parameter :: & LATTICE_FCCTOBCC_BAINVARIANT = reshape( [& 1, 0, 0, 0, 1, 0, 0, 0, 1, & ! Pitsch OR (Ma & Hartmaier 2014, Table 3) @@ -2283,7 +2261,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc) 0, 0, 1, 1, 0, 0, 0, 1, 0, & 0, 0, 1, 1, 0, 0, 0, 1, 0 & ],shape(LATTICE_FCCTOBCC_BAINVARIANT)) - + real(pReal), dimension(4,LATTICE_fcc_Ntrans), parameter :: & LATTICE_FCCTOBCC_BAINROT = reshape([& 1.0, 0.0, 0.0, 45.0, & ! Rotate fcc austensite to bain variant @@ -2299,7 +2277,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc) 0.0, 0.0, 1.0, 45.0, & 0.0, 0.0, 1.0, 45.0 & ],shape(LATTICE_FCCTOBCC_BAINROT)) - + if (a_bcc > 0.0_pReal .and. a_fcc > 0.0_pReal .and. dEq0(cOverA)) then ! fcc -> bcc transformation do i = 1,sum(Ntrans) call R%fromAxisAngle(LATTICE_FCCTOBCC_SYSTEMTRANS(:,i),degrees=.true.,P=1) @@ -2307,7 +2285,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc) x = real(LATTICE_FCCTOBCC_BAINVARIANT(1:3,i),pReal) y = real(LATTICE_FCCTOBCC_BAINVARIANT(4:6,i),pReal) z = real(LATTICE_FCCTOBCC_BAINVARIANT(7:9,i),pReal) - + U = (a_bcc/a_fcc)*math_outer(x,x) & + (a_bcc/a_fcc)*math_outer(y,y) * sqrt(2.0_pReal) & + (a_bcc/a_fcc)*math_outer(z,z) * sqrt(2.0_pReal) @@ -2319,7 +2297,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc) sd = MATH_I3 ss(1,3) = sqrt(2.0_pReal)/4.0_pReal sd(3,3) = cOverA/sqrt(8.0_pReal/3.0_pReal) - + do i = 1,sum(Ntrans) x = LATTICE_FCCTOHEX_SYSTEMTRANS(1:3,i)/norm2(LATTICE_FCCTOHEX_SYSTEMTRANS(1:3,i)) z = LATTICE_FCCTOHEX_SYSTEMTRANS(4:6,i)/norm2(LATTICE_FCCTOHEX_SYSTEMTRANS(4:6,i)) @@ -2340,7 +2318,7 @@ end subroutine buildTransformationSystem !> @brief select active systems as strings !-------------------------------------------------------------------------------------------------- function getlabels(active,potential,system) result(labels) - + integer, dimension(:), intent(in) :: & active, & !< # of active systems per family potential !< # of potential systems per family @@ -2350,7 +2328,7 @@ function getlabels(active,potential,system) result(labels) character(len=:), dimension(:), allocatable :: labels character(len=:), allocatable :: label - integer :: i,j + integer :: i,j integer :: & a, & !< index of active system p, & !< index in potential system matrix @@ -2359,13 +2337,13 @@ function getlabels(active,potential,system) result(labels) i = 2*size(system,1) + (size(system,1) - 2) + 4 ! 2 letters per index + spaces + brackets allocate(character(len=i) :: labels(sum(active)), label) - + a = 0 activeFamilies: do f = 1,size(active,1) activeSystems: do s = 1,active(f) a = a + 1 p = sum(potential(1:f-1))+s - + i = 1 label(i:i) = '[' direction: do j = 1, size(system,1)/2 @@ -2374,7 +2352,7 @@ function getlabels(active,potential,system) result(labels) i = i + 3 enddo direction label(i:i) = ']' - + i = i +1 label(i:i) = '(' normal: do j = size(system,1)/2+1, size(system,1) @@ -2383,12 +2361,12 @@ function getlabels(active,potential,system) result(labels) i = i + 3 enddo normal label(i:i) = ')' - + labels(s) = label enddo activeSystems enddo activeFamilies - + end function getlabels diff --git a/src/material.f90 b/src/material.f90 index 7e68049b8..9c5009484 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -257,7 +257,7 @@ subroutine material_init allocate(damage (material_Nhomogenization)) allocate(temperatureRate (material_Nhomogenization)) - + do m = 1,size(config_microstructure) if(minval(microstructure_phase(1:microstructure_Nconstituents(m),m)) < 1 .or. & maxval(microstructure_phase(1:microstructure_Nconstituents(m),m)) > size(config_phase)) & @@ -268,6 +268,7 @@ subroutine material_init if(microstructure_Nconstituents(m) < 1) & call IO_error(151,m) enddo + if(homogenization_maxNgrains > size(microstructure_phase,1)) call IO_error(148) debugOut: if (iand(myDebug,debug_levelExtensive) /= 0) then write(6,'(/,a,/)') ' MATERIAL configuration' @@ -290,6 +291,7 @@ subroutine material_init !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! new mappings + allocate(material_phaseAt(homogenization_maxNgrains,discretization_nElem), source=0) allocate(material_texture(homogenization_maxNgrains,discretization_nIP,discretization_nElem),source=0) !this is only needed by plasticity nonlocal allocate(material_orientation0(homogenization_maxNgrains,discretization_nIP,discretization_nElem)) @@ -298,15 +300,24 @@ subroutine material_init do i = 1, discretization_nIP myMicro = discretization_microstructureAt(e) do c = 1, homogenization_Ngrains(discretization_homogenizationAt(e)) - material_phaseAt(c,e) = microstructure_phase(c,myMicro) - material_texture(c,i,e) = microstructure_texture(c,myMicro) - material_orientation0(c,i,e) = texture_orientation(material_texture(c,i,e)) + if(microstructure_phase(c,myMicro) > 0) then + material_phaseAt(c,e) = microstructure_phase(c,myMicro) + else + call IO_error(150,ext_msg='phase') + endif + if(microstructure_texture(c,myMicro) > 0) then + material_texture(c,i,e) = microstructure_texture(c,myMicro) + material_orientation0(c,i,e) = texture_orientation(material_texture(c,i,e)) + else + call IO_error(150,ext_msg='texture') + endif enddo enddo enddo deallocate(microstructure_phase) deallocate(microstructure_texture) + deallocate(texture_orientation) allocate(material_homogenizationAt,source=discretization_homogenizationAt) @@ -464,7 +475,7 @@ subroutine material_parseMicrostructure real(pReal), dimension(:,:), allocatable :: & microstructure_fraction !< vol fraction of each constituent in microstructure integer :: & - microstructure_maxNconstituents !< max number of constituents in any phase + maxNconstituents !< max number of constituents in any phase allocate(microstructure_Nconstituents(size(config_microstructure)), source=0) @@ -475,10 +486,10 @@ subroutine material_parseMicrostructure microstructure_Nconstituents(m) = config_microstructure(m)%countKeys('(constituent)') enddo - microstructure_maxNconstituents = maxval(microstructure_Nconstituents) - allocate(microstructure_phase (microstructure_maxNconstituents,size(config_microstructure)),source=0) - allocate(microstructure_texture (microstructure_maxNconstituents,size(config_microstructure)),source=0) - allocate(microstructure_fraction(microstructure_maxNconstituents,size(config_microstructure)),source=0.0_pReal) + maxNconstituents = maxval(microstructure_Nconstituents) + allocate(microstructure_phase (maxNconstituents,size(config_microstructure)),source=0) + allocate(microstructure_texture (maxNconstituents,size(config_microstructure)),source=0) + allocate(microstructure_fraction(maxNconstituents,size(config_microstructure)),source=0.0_pReal) allocate(strings(1)) ! Intel 16.0 Bug do m=1, size(config_microstructure) diff --git a/src/results.f90 b/src/results.f90 index 8708f7856..cfd495a71 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -485,15 +485,15 @@ end subroutine results_writeScalarDataset_rotation !-------------------------------------------------------------------------------------------------- !> @brief adds the unique mapping from spatial position and constituent ID to results !-------------------------------------------------------------------------------------------------- -subroutine results_mapping_constituent(phaseAt,memberAt,label) +subroutine results_mapping_constituent(phaseAt,memberAtLocal,label) integer, dimension(:,:), intent(in) :: phaseAt !< phase section at (constituent,element) - integer, dimension(:,:,:), intent(in) :: memberAt !< phase member at (constituent,IP,element) - character(len=*), dimension(:), intent(in) :: label !< label of each phase section + integer, dimension(:,:,:), intent(in) :: memberAtLocal !< phase member at (constituent,IP,element) + character(len=pStringLen), dimension(:), intent(in) :: label !< label of each phase section - integer, dimension(size(memberAt,1),size(memberAt,2),size(memberAt,3)) :: & - phaseAt_perIP, & - memberAt_total + integer, dimension(size(memberAtLocal,1),size(memberAtLocal,2),size(memberAtLocal,3)) :: & + phaseAtMaterialpoint, & + memberAtGlobal integer, dimension(size(label),0:worldsize-1) :: memberOffset !< offset in member counting per process integer, dimension(0:worldsize-1) :: writeSize !< amount of data written per process integer(HSIZE_T), dimension(2) :: & @@ -543,10 +543,10 @@ subroutine results_mapping_constituent(phaseAt,memberAt,label) call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, ierr) memberOffset = 0 do i=1, size(label) - memberOffset(i,worldrank) = count(phaseAt == i)*size(memberAt,2) ! number of points/instance of this process + memberOffset(i,worldrank) = count(phaseAt == i)*size(memberAtLocal,2) ! number of points/instance of this process enddo writeSize = 0 - writeSize(worldrank) = size(memberAt(1,:,:)) ! total number of points by this process + writeSize(worldrank) = size(memberAtLocal(1,:,:)) ! total number of points by this process !-------------------------------------------------------------------------------------------------- ! MPI settings and communication @@ -578,14 +578,14 @@ subroutine results_mapping_constituent(phaseAt,memberAt,label) !--------------------------------------------------------------------------------------------------- ! expand phaseAt to consider IPs (is not stored per IP) - do i = 1, size(phaseAt_perIP,2) - phaseAt_perIP(:,i,:) = phaseAt + do i = 1, size(phaseAtMaterialpoint,2) + phaseAtMaterialpoint(:,i,:) = phaseAt enddo !--------------------------------------------------------------------------------------------------- ! renumber member from my process to all processes do i = 1, size(label) - where(phaseAt_perIP == i) memberAt_total = memberAt + sum(memberOffset(i,0:worldrank-1)) -1 ! convert to 0-based + where(phaseAtMaterialpoint == i) memberAtGlobal = memberAtLocal + sum(memberOffset(i,0:worldrank-1)) -1 ! convert to 0-based enddo !-------------------------------------------------------------------------------------------------- @@ -596,10 +596,10 @@ subroutine results_mapping_constituent(phaseAt,memberAt,label) call h5dcreate_f(loc_id, 'constituent', dtype_id, filespace_id, dset_id, ierr) if (ierr < 0) call IO_error(1,ext_msg='results_mapping_constituent: h5dcreate_f') - call h5dwrite_f(dset_id, name_id, reshape(label(pack(phaseAt_perIP,.true.)),myShape), & + call h5dwrite_f(dset_id, name_id, reshape(label(pack(phaseAtMaterialpoint,.true.)),myShape), & myShape, ierr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) if (ierr < 0) call IO_error(1,ext_msg='results_mapping_constituent: h5dwrite_f/name_id') - call h5dwrite_f(dset_id, position_id, reshape(pack(memberAt_total,.true.),myShape), & + call h5dwrite_f(dset_id, position_id, reshape(pack(memberAtGlobal,.true.),myShape), & myShape, ierr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) if (ierr < 0) call IO_error(1,ext_msg='results_mapping_constituent: h5dwrite_f/position_id') @@ -620,15 +620,15 @@ end subroutine results_mapping_constituent !-------------------------------------------------------------------------------------------------- !> @brief adds the unique mapping from spatial position and constituent ID to results !-------------------------------------------------------------------------------------------------- -subroutine results_mapping_materialpoint(homogenizationAt,memberAt,label) +subroutine results_mapping_materialpoint(homogenizationAt,memberAtLocal,label) integer, dimension(:), intent(in) :: homogenizationAt !< homogenization section at (element) - integer, dimension(:,:), intent(in) :: memberAt !< homogenization member at (IP,element) - character(len=*), dimension(:), intent(in) :: label !< label of each homogenization section + integer, dimension(:,:), intent(in) :: memberAtLocal !< homogenization member at (IP,element) + character(len=pStringLen), dimension(:), intent(in) :: label !< label of each homogenization section - integer, dimension(size(memberAt,1),size(memberAt,2)) :: & - homogenizationAt_perIP, & - memberAt_total + integer, dimension(size(memberAtLocal,1),size(memberAtLocal,2)) :: & + homogenizationAtMaterialpoint, & + memberAtGlobal integer, dimension(size(label),0:worldsize-1) :: memberOffset !< offset in member counting per process integer, dimension(0:worldsize-1) :: writeSize !< amount of data written per process integer(HSIZE_T), dimension(1) :: & @@ -678,10 +678,10 @@ subroutine results_mapping_materialpoint(homogenizationAt,memberAt,label) call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, ierr) memberOffset = 0 do i=1, size(label) - memberOffset(i,worldrank) = count(homogenizationAt == i)*size(memberAt,1) ! number of points/instance of this process + memberOffset(i,worldrank) = count(homogenizationAt == i)*size(memberAtLocal,1) ! number of points/instance of this process enddo writeSize = 0 - writeSize(worldrank) = size(memberAt) ! total number of points by this process + writeSize(worldrank) = size(memberAtLocal) ! total number of points by this process !-------------------------------------------------------------------------------------------------- ! MPI settings and communication @@ -713,14 +713,14 @@ subroutine results_mapping_materialpoint(homogenizationAt,memberAt,label) !--------------------------------------------------------------------------------------------------- ! expand phaseAt to consider IPs (is not stored per IP) - do i = 1, size(homogenizationAt_perIP,1) - homogenizationAt_perIP(i,:) = homogenizationAt + do i = 1, size(homogenizationAtMaterialpoint,1) + homogenizationAtMaterialpoint(i,:) = homogenizationAt enddo !--------------------------------------------------------------------------------------------------- ! renumber member from my process to all processes do i = 1, size(label) - where(homogenizationAt_perIP == i) memberAt_total = memberAt + sum(memberOffset(i,0:worldrank-1)) - 1 ! convert to 0-based + where(homogenizationAtMaterialpoint == i) memberAtGlobal = memberAtLocal + sum(memberOffset(i,0:worldrank-1)) - 1 ! convert to 0-based enddo !-------------------------------------------------------------------------------------------------- @@ -731,10 +731,10 @@ subroutine results_mapping_materialpoint(homogenizationAt,memberAt,label) call h5dcreate_f(loc_id, 'materialpoint', dtype_id, filespace_id, dset_id, ierr) if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5dcreate_f') - call h5dwrite_f(dset_id, name_id, reshape(label(pack(homogenizationAt_perIP,.true.)),myShape), & + call h5dwrite_f(dset_id, name_id, reshape(label(pack(homogenizationAtMaterialpoint,.true.)),myShape), & myShape, ierr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5dwrite_f/name_id') - call h5dwrite_f(dset_id, position_id, reshape(pack(memberAt_total,.true.),myShape), & + call h5dwrite_f(dset_id, position_id, reshape(pack(memberAtGlobal,.true.),myShape), & myShape, ierr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5dwrite_f/position_id')