diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index bb4a2e7a0..40662cb40 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -344,6 +344,7 @@ Phenopowerlaw_singleSlip: Pytest_grid: stage: grid script: + - module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel - cd pytest - pytest except: diff --git a/PRIVATE b/PRIVATE index 2550bf655..dc568df60 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 2550bf655aebe61ce459718323ba38b1ac205a7f +Subproject commit dc568df60e36b659d9a1f84ac93fd4abb1b8fe3c diff --git a/VERSION b/VERSION index a8cc04d22..425e852a2 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v3.0.0-alpha-321-g0f6495430 +v3.0.0-alpha-335-gdb8f6400f diff --git a/python/damask/__init__.py b/python/damask/__init__.py index 3f2ff6813..e21ebd03d 100644 --- a/python/damask/__init__.py +++ b/python/damask/__init__.py @@ -20,6 +20,12 @@ from ._result import Result # noqa from ._geom import Geom # noqa from ._material import Material # noqa from . import solver # noqa +from . import util # noqa +from . import seeds # noqa +from . import grid_filters # noqa +from . import mechanics # noqa + + # deprecated Environment = _ diff --git a/python/damask/_geom.py b/python/damask/_geom.py index 6af265021..50d4e2b33 100644 --- a/python/damask/_geom.py +++ b/python/damask/_geom.py @@ -212,7 +212,7 @@ class Geom: return np.argmin(np.sum((np.broadcast_to(point,(len(seeds),3))-seeds)**2,axis=1) - weights) @staticmethod - def from_Laguerre_tessellation(grid,size,seeds,weights,periodic=True): + def from_Laguerre_tessellation(grid,size,seeds,weights,material=None,periodic=True): """ Generate geometry from Laguerre tessellation. @@ -226,6 +226,9 @@ class Geom: Position of the seed points in meter. All points need to lay within the box. weights : numpy.ndarray of shape (seeds.shape[0]) Weights of the seeds. Setting all weights to 1.0 gives a standard Voronoi tessellation. + material : numpy.ndarray of shape (seeds.shape[0]), optional + Material ID of the seeds. Defaults to None, in which case materials are + consecutively numbered. periodic : Boolean, optional Perform a periodic tessellation. Defaults to True. @@ -245,22 +248,22 @@ class Geom: result = pool.map_async(partial(Geom._find_closest_seed,seeds_p,weights_p), [coord for coord in coords]) pool.close() pool.join() - material = np.array(result.get()) + material_ = np.array(result.get()) if periodic: - material = material.reshape(grid*3) - material = material[grid[0]:grid[0]*2,grid[1]:grid[1]*2,grid[2]:grid[2]*2]%seeds.shape[0] + material_ = material_.reshape(grid*3) + material_ = material_[grid[0]:grid[0]*2,grid[1]:grid[1]*2,grid[2]:grid[2]*2]%seeds.shape[0] else: - material = material.reshape(grid) + material_ = material_.reshape(grid) - return Geom(material = material+1, + return Geom(material = material_+1 if material is None else material[material_], size = size, comments = util.execution_stamp('Geom','from_Laguerre_tessellation'), ) @staticmethod - def from_Voronoi_tessellation(grid,size,seeds,periodic=True): + def from_Voronoi_tessellation(grid,size,seeds,material=None,periodic=True): """ Generate geometry from Voronoi tessellation. @@ -272,15 +275,18 @@ class Geom: Physical size of the geometry in meter. seeds : numpy.ndarray of shape (:,3) Position of the seed points in meter. All points need to lay within the box. + material : numpy.ndarray of shape (seeds.shape[0]), optional + Material ID of the seeds. Defaults to None, in which case materials are + consecutively numbered. periodic : Boolean, optional Perform a periodic tessellation. Defaults to True. """ coords = grid_filters.cell_coord0(grid,size).reshape(-1,3) KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds) - devNull,material = KDTree.query(coords) + devNull,material_ = KDTree.query(coords) - return Geom(material = material.reshape(grid)+1, + return Geom(material = (material_+1 if material is None else material[material_]).reshape(grid), size = size, comments = util.execution_stamp('Geom','from_Voronoi_tessellation'), ) diff --git a/python/damask/seeds.py b/python/damask/seeds.py new file mode 100644 index 000000000..ea40794da --- /dev/null +++ b/python/damask/seeds.py @@ -0,0 +1,113 @@ +from scipy import spatial as _spatial +import numpy as _np + +from . import util +from . import grid_filters + + +def from_random(size,N_seeds,grid=None,seed=None): + """ + Random seeding in space. + + Parameters + ---------- + size : numpy.ndarray of shape (3) + Physical size of the seeding domain. + N_seeds : int + Number of seeds. + grid : numpy.ndarray of shape (3), optional. + If given, ensures that all seeds initiate one grain if using a + standard Voronoi tessellation. + seed : {None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional + A seed to initialize the BitGenerator. Defaults to None. + If None, then fresh, unpredictable entropy will be pulled from the OS. + + """ + rng = _np.random.default_rng(seed) + if grid is None: + coords = rng.random((N_seeds,3)) * size + else: + grid_coords = grid_filters.cell_coord0(grid,size).reshape(-1,3,order='F') + coords = grid_coords[rng.choice(_np.prod(grid),N_seeds, replace=False)] \ + + _np.broadcast_to(size/grid,(N_seeds,3))*(rng.random((N_seeds,3))*.5-.25) # wobble without leaving grid + + return coords + + +def from_Poisson_disc(size,N_seeds,N_candidates,distance,periodic=True,seed=None): + """ + Seeding in space according to a Poisson disc distribution. + + Parameters + ---------- + size : numpy.ndarray of shape (3) + Physical size of the seeding domain. + N_seeds : int + Number of seeds. + N_candidates : int + Number of candidates to consider for finding best candidate. + distance : float + Minimum acceptable distance to other seeds. + periodic : boolean, optional + Calculate minimum distance for periodically repeated grid. + seed : {None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional + A seed to initialize the BitGenerator. Defaults to None. + If None, then fresh, unpredictable entropy will be pulled from the OS. + + """ + rng = _np.random.default_rng(seed) + coords = _np.empty((N_seeds,3)) + coords[0] = rng.random(3) * size + + i = 1 + progress = util._ProgressBar(N_seeds+1,'',50) + while i < N_seeds: + candidates = rng.random((N_candidates,3))*_np.broadcast_to(size,(N_candidates,3)) + tree = _spatial.cKDTree(coords[:i],boxsize=size) if periodic else \ + _spatial.cKDTree(coords[:i]) + distances, dev_null = tree.query(candidates) + best = distances.argmax() + if distances[best] > distance: # require minimum separation + coords[i] = candidates[best] # maximum separation to existing point cloud + i += 1 + progress.update(i) + + return coords + + +def from_geom(geom,selection=None,invert=False,average=False,periodic=True): + """ + Create seed from existing geometry description. + + Parameters + ---------- + geom : damask.Geom + Geometry, from which the material IDs are used as seeds. + selection : iterable of integers, optional + Material IDs to consider. + invert : boolean, false + Do not consider the material IDs given in selection. Defaults to False. + average : boolean, optional + Seed corresponds to center of gravity of material ID cloud. + periodic : boolean, optional + Center of gravity with periodic boundaries. + + """ + material = geom.material.reshape((-1,1),order='F') + mask = _np.full(geom.grid.prod(),True,dtype=bool) if selection is None else \ + _np.isin(material,selection,invert=invert) + coords = grid_filters.cell_coord0(geom.grid,geom.size).reshape(-1,3,order='F') + + if not average: + return (coords[mask],material[mask]) + else: + materials = _np.unique(material[mask]) + coords_ = _np.zeros((materials.size,3),dtype=float) + for i,mat in enumerate(materials): + pc = (2*_np.pi*coords[material[:,0]==mat,:]-geom.origin)/geom.size + coords_[i] = geom.origin + geom.size / 2 / _np.pi * (_np.pi + + _np.arctan2(-_np.average(_np.sin(pc),axis=0), + -_np.average(_np.cos(pc),axis=0))) \ + if periodic else \ + _np.average(coords[material[:,0]==mat,:],axis=0) + return (coords_,materials) diff --git a/python/tests/test_Geom.py b/python/tests/test_Geom.py index ac4a8eac3..9aff60186 100644 --- a/python/tests/test_Geom.py +++ b/python/tests/test_Geom.py @@ -83,6 +83,13 @@ class TestGeom: with pytest.raises(ValueError): Geom.load(tmpdir/'no_materialpoint.vtr') + def test_invalid_material(self): + with pytest.raises(TypeError): + Geom(np.zeros((3,3,3),dtype='complex'),np.ones(3)) + + def test_cast_to_int(self): + g = Geom(np.zeros((3,3,3)),np.ones(3)) + assert g.material.dtype in np.sctypes['int'] @pytest.mark.parametrize('compress',[True,False]) def test_compress(self,default,tmpdir,compress): @@ -324,8 +331,8 @@ class TestGeom: size = np.random.random(3) + 1.0 N_seeds= np.random.randint(10,30) seeds = np.random.rand(N_seeds,3) * np.broadcast_to(size,(N_seeds,3)) - Voronoi = Geom.from_Voronoi_tessellation( grid,size,seeds, periodic) - Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,np.ones(N_seeds),periodic) + Voronoi = Geom.from_Voronoi_tessellation( grid,size,seeds, np.arange(N_seeds)+5,periodic) + Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,np.ones(N_seeds),np.arange(N_seeds)+5,periodic) assert geom_equal(Laguerre,Voronoi) @@ -337,7 +344,7 @@ class TestGeom: weights= np.full((N_seeds),-np.inf) ms = np.random.randint(1, N_seeds+1) weights[ms-1] = np.random.random() - Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,weights,np.random.random()>0.5) + Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,weights,periodic=np.random.random()>0.5) assert np.all(Laguerre.material == ms) @@ -349,7 +356,7 @@ class TestGeom: material = np.ones(grid) material[:,grid[1]//2:,:] = 2 if approach == 'Laguerre': - geom = Geom.from_Laguerre_tessellation(grid,size,seeds,np.ones(2),np.random.random()>0.5) + geom = Geom.from_Laguerre_tessellation(grid,size,seeds,np.ones(2),periodic=np.random.random()>0.5) elif approach == 'Voronoi': - geom = Geom.from_Voronoi_tessellation(grid,size,seeds, np.random.random()>0.5) + geom = Geom.from_Voronoi_tessellation(grid,size,seeds, periodic=np.random.random()>0.5) assert np.all(geom.material == material) diff --git a/python/tests/test_seeds.py b/python/tests/test_seeds.py new file mode 100644 index 000000000..49d9a6bfd --- /dev/null +++ b/python/tests/test_seeds.py @@ -0,0 +1,36 @@ +import pytest +import numpy as np +from scipy.spatial import cKDTree + +from damask import seeds +from damask import Geom + +class TestSeeds: + + @pytest.mark.parametrize('grid',[None,np.ones(3,dtype='i')*10]) + def test_from_random(self,grid): + N_seeds = np.random.randint(30,300) + size = np.ones(3) + np.random.random(3) + coords = seeds.from_random(size,N_seeds,grid) + assert (0<=coords).all() and (coords=distance + + def test_from_geom(self): + grid = np.random.randint(10,20,3) + N_seeds = np.random.randint(30,300) + size = np.ones(3) + np.random.random(3) + coords = seeds.from_random(size,N_seeds,grid) + geom_1 = Geom.from_Voronoi_tessellation(grid,size,coords) + coords,material = seeds.from_geom(geom_1) + geom_2 = Geom.from_Voronoi_tessellation(grid,size,coords,material) + assert (geom_2.material==geom_1.material).all() diff --git a/src/constitutive_plastic_disloTungsten.f90 b/src/constitutive_plastic_disloTungsten.f90 index 98a9ce2f6..54c01b912 100644 --- a/src/constitutive_plastic_disloTungsten.f90 +++ b/src/constitutive_plastic_disloTungsten.f90 @@ -17,18 +17,18 @@ submodule(constitutive:constitutive_plastic) plastic_disloTungsten D_0 = 1.0_pReal, & !< prefactor for self-diffusion coefficient Q_cl = 1.0_pReal !< activation energy for dislocation climb real(pReal), allocatable, dimension(:) :: & - b_sl, & !< magnitude of burgers vector [m] + b_sl, & !< magnitude of Burgers vector [m] D_a, & i_sl, & !< Adj. parameter for distance between 2 forest dislocations - atomicVolume, & !< factor to calculate atomic volume - tau_0, & !< Peierls stress + f_at, & !< factor to calculate atomic volume + tau_Peierls, & !< Peierls stress !* mobility law parameters - delta_F, & !< activation energy for glide [J] - v0, & !< dislocation velocity prefactor [m/s] + Q_s, & !< activation energy for glide [J] + v_0, & !< dislocation velocity prefactor [m/s] p, & !< p-exponent in glide velocity q, & !< q-exponent in glide velocity B, & !< friction coefficient - kink_height, & !< height of the kink pair + h, & !< height of the kink pair w, & !< width of the kink pair omega !< attempt frequency for kink pair nucleation real(pReal), allocatable, dimension(:,:) :: & @@ -142,7 +142,7 @@ module function plastic_disloTungsten_init() result(myPlasticity) phase%get_asFloat('c/a',defaultVal=0.0_pReal)) if(trim(phase%get_asString('lattice')) == 'bcc') then - a = pl%get_asFloats('nonSchmid_coefficients',defaultVal = emptyRealArray) + a = pl%get_asFloats('a_nonSchmid',defaultVal = emptyRealArray) prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,a,+1) prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,a,-1) else @@ -158,17 +158,17 @@ module function plastic_disloTungsten_init() result(myPlasticity) rho_mob_0 = pl%get_asFloats('rho_mob_0', requiredSize=size(N_sl)) rho_dip_0 = pl%get_asFloats('rho_dip_0', requiredSize=size(N_sl)) - prm%v0 = pl%get_asFloats('v_0', requiredSize=size(N_sl)) + prm%v_0 = pl%get_asFloats('v_0', requiredSize=size(N_sl)) prm%b_sl = pl%get_asFloats('b_sl', requiredSize=size(N_sl)) - prm%delta_F = pl%get_asFloats('Q_s', requiredSize=size(N_sl)) + prm%Q_s = pl%get_asFloats('Q_s', requiredSize=size(N_sl)) prm%i_sl = pl%get_asFloats('i_sl', requiredSize=size(N_sl)) - prm%tau_0 = pl%get_asFloats('tau_peierls', requiredSize=size(N_sl)) + prm%tau_Peierls = pl%get_asFloats('tau_Peierls', requiredSize=size(N_sl)) prm%p = pl%get_asFloats('p_sl', requiredSize=size(N_sl), & defaultVal=[(1.0_pReal,i=1,size(N_sl))]) prm%q = pl%get_asFloats('q_sl', requiredSize=size(N_sl), & defaultVal=[(1.0_pReal,i=1,size(N_sl))]) - prm%kink_height = pl%get_asFloats('h', requiredSize=size(N_sl)) + prm%h = pl%get_asFloats('h', requiredSize=size(N_sl)) prm%w = pl%get_asFloats('w', requiredSize=size(N_sl)) prm%omega = pl%get_asFloats('omega', requiredSize=size(N_sl)) prm%B = pl%get_asFloats('B', requiredSize=size(N_sl)) @@ -176,7 +176,7 @@ module function plastic_disloTungsten_init() result(myPlasticity) prm%D = pl%get_asFloat('D') prm%D_0 = pl%get_asFloat('D_0') prm%Q_cl = pl%get_asFloat('Q_cl') - prm%atomicVolume = pl%get_asFloat('f_at') * prm%b_sl**3.0_pReal + prm%f_at = pl%get_asFloat('f_at') * prm%b_sl**3.0_pReal prm%D_a = pl%get_asFloat('D_a') * prm%b_sl prm%dipoleformation = pl%get_asBool('dipole_formation_factor', defaultVal = .true.) @@ -186,16 +186,16 @@ module function plastic_disloTungsten_init() result(myPlasticity) rho_dip_0 = math_expand(rho_dip_0, N_sl) prm%q = math_expand(prm%q, N_sl) prm%p = math_expand(prm%p, N_sl) - prm%delta_F = math_expand(prm%delta_F, N_sl) + prm%Q_s = math_expand(prm%Q_s, N_sl) prm%b_sl = math_expand(prm%b_sl, N_sl) - prm%kink_height = math_expand(prm%kink_height, N_sl) + prm%h = math_expand(prm%h, N_sl) prm%w = math_expand(prm%w, N_sl) prm%omega = math_expand(prm%omega, N_sl) - prm%tau_0 = math_expand(prm%tau_0, N_sl) - prm%v0 = math_expand(prm%v0, N_sl) + prm%tau_Peierls = math_expand(prm%tau_Peierls, N_sl) + prm%v_0 = math_expand(prm%v_0, N_sl) prm%B = math_expand(prm%B, N_sl) prm%i_sl = math_expand(prm%i_sl, N_sl) - prm%atomicVolume = math_expand(prm%atomicVolume, N_sl) + prm%f_at = math_expand(prm%f_at, N_sl) prm%D_a = math_expand(prm%D_a, N_sl) ! sanity checks @@ -203,17 +203,17 @@ module function plastic_disloTungsten_init() result(myPlasticity) if ( prm%Q_cl <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl' if (any(rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_mob_0' if (any(rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_dip_0' - if (any(prm%v0 < 0.0_pReal)) extmsg = trim(extmsg)//' v_0' + if (any(prm%v_0 < 0.0_pReal)) extmsg = trim(extmsg)//' v_0' if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' b_sl' - if (any(prm%delta_F <= 0.0_pReal)) extmsg = trim(extmsg)//' Q_s' - if (any(prm%tau_0 < 0.0_pReal)) extmsg = trim(extmsg)//' tau_peierls' + if (any(prm%Q_s <= 0.0_pReal)) extmsg = trim(extmsg)//' Q_s' + if (any(prm%tau_Peierls < 0.0_pReal)) extmsg = trim(extmsg)//' tau_Peierls' if (any(prm%D_a <= 0.0_pReal)) extmsg = trim(extmsg)//' D_a or b_sl' - if (any(prm%atomicVolume <= 0.0_pReal)) extmsg = trim(extmsg)//' f_at or b_sl' + if (any(prm%f_at <= 0.0_pReal)) extmsg = trim(extmsg)//' f_at or b_sl' else slipActive rho_mob_0= emptyRealArray; rho_dip_0 = emptyRealArray - allocate(prm%b_sl,prm%D_a,prm%i_sl,prm%atomicVolume,prm%tau_0, & - prm%delta_F,prm%v0,prm%p,prm%q,prm%B,prm%kink_height,prm%w,prm%omega, & + allocate(prm%b_sl,prm%D_a,prm%i_sl,prm%f_at,prm%tau_Peierls, & + prm%Q_s,prm%v_0,prm%p,prm%q,prm%B,prm%h,prm%w,prm%omega, & source = emptyRealArray) allocate(prm%forestProjection(0,0)) allocate(prm%h_sl_sl (0,0)) @@ -354,7 +354,7 @@ module subroutine plastic_disloTungsten_dotState(Mp,T,instance,of) dot_rho_dip_formation = merge(2.0_pReal*dip_distance* stt%rho_mob(:,of)*abs(dot%gamma_sl(:,of))/prm%b_sl, & ! ToDo: ignore region of spontaneous annihilation 0.0_pReal, & prm%dipoleformation) - v_cl = (3.0_pReal*prm%mu*VacancyDiffusion*prm%atomicVolume/(2.0_pReal*pi*kB*T)) & + v_cl = (3.0_pReal*prm%mu*VacancyDiffusion*prm%f_at/(2.0_pReal*pi*kB*T)) & * (1.0_pReal/(dip_distance+prm%D_a)) dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,of))/(dip_distance-prm%D_a) ! ToDo: Discuss with Franz: Stress dependency? end where @@ -477,12 +477,12 @@ pure subroutine kinetics(Mp,T,instance,of, & if (present(tau_pos_out)) tau_pos_out = tau_pos if (present(tau_neg_out)) tau_neg_out = tau_neg - associate(BoltzmannRatio => prm%delta_F/(kB*T), & - dot_gamma_0 => stt%rho_mob(:,of)*prm%b_sl*prm%v0, & + associate(BoltzmannRatio => prm%Q_s/(kB*T), & + dot_gamma_0 => stt%rho_mob(:,of)*prm%b_sl*prm%v_0, & effectiveLength => dst%Lambda_sl(:,of) - prm%w) significantPositiveTau: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check) - StressRatio = (abs(tau_pos)-dst%threshold_stress(:,of))/prm%tau_0 + StressRatio = (abs(tau_pos)-dst%threshold_stress(:,of))/prm%tau_Peierls StressRatio_p = StressRatio** prm%p StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) @@ -490,7 +490,7 @@ pure subroutine kinetics(Mp,T,instance,of, & t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength) t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos) - vel = prm%kink_height/(t_n + t_k) + vel = prm%h/(t_n + t_k) dot_gamma_pos = dot_gamma_0 * sign(vel,tau_pos) * 0.5_pReal else where significantPositiveTau @@ -500,10 +500,10 @@ pure subroutine kinetics(Mp,T,instance,of, & if (present(ddot_gamma_dtau_pos)) then significantPositiveTau2: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check) dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & - * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_0 + * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls dtk = -1.0_pReal * t_k / tau_pos - dvel = -1.0_pReal * prm%kink_height * (dtk + dtn) / (t_n + t_k)**2.0_pReal + dvel = -1.0_pReal * prm%h * (dtk + dtn) / (t_n + t_k)**2.0_pReal ddot_gamma_dtau_pos = dot_gamma_0 * dvel* 0.5_pReal else where significantPositiveTau2 @@ -512,7 +512,7 @@ pure subroutine kinetics(Mp,T,instance,of, & endif significantNegativeTau: where(abs(tau_neg)-dst%threshold_stress(:,of) > tol_math_check) - StressRatio = (abs(tau_neg)-dst%threshold_stress(:,of))/prm%tau_0 + StressRatio = (abs(tau_neg)-dst%threshold_stress(:,of))/prm%tau_Peierls StressRatio_p = StressRatio** prm%p StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) @@ -520,7 +520,7 @@ pure subroutine kinetics(Mp,T,instance,of, & t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength) t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos) - vel = prm%kink_height/(t_n + t_k) + vel = prm%h/(t_n + t_k) dot_gamma_neg = dot_gamma_0 * sign(vel,tau_neg) * 0.5_pReal else where significantNegativeTau @@ -530,10 +530,10 @@ pure subroutine kinetics(Mp,T,instance,of, & if (present(ddot_gamma_dtau_neg)) then significantNegativeTau2: where(abs(tau_neg)-dst%threshold_stress(:,of) > tol_math_check) dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & - * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_0 + * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls dtk = -1.0_pReal * t_k / tau_neg - dvel = -1.0_pReal * prm%kink_height * (dtk + dtn) / (t_n + t_k)**2.0_pReal + dvel = -1.0_pReal * prm%h * (dtk + dtn) / (t_n + t_k)**2.0_pReal ddot_gamma_dtau_neg = dot_gamma_0 * dvel * 0.5_pReal else where significantNegativeTau2 diff --git a/src/constitutive_plastic_dislotwin.f90 b/src/constitutive_plastic_dislotwin.f90 index 4890316b8..62dcdd83e 100644 --- a/src/constitutive_plastic_dislotwin.f90 +++ b/src/constitutive_plastic_dislotwin.f90 @@ -16,38 +16,38 @@ submodule(constitutive:constitutive_plastic) plastic_dislotwin real(pReal) :: & mu = 1.0_pReal, & !< equivalent shear modulus nu = 1.0_pReal, & !< equivalent shear Poisson's ratio - D0 = 1.0_pReal, & !< prefactor for self-diffusion coefficient - Qsd = 1.0_pReal, & !< activation energy for dislocation climb + D_0 = 1.0_pReal, & !< prefactor for self-diffusion coefficient + Q_cl = 1.0_pReal, & !< activation energy for dislocation climb omega = 1.0_pReal, & !< frequency factor for dislocation climb D = 1.0_pReal, & !< grain size p_sb = 1.0_pReal, & !< p-exponent in shear band velocity q_sb = 1.0_pReal, & !< q-exponent in shear band velocity - CEdgeDipMinDistance = 1.0_pReal, & !< adjustment parameter to calculate minimum dipole distance + D_a = 1.0_pReal, & !< adjustment parameter to calculate minimum dipole distance i_tw = 1.0_pReal, & !< adjustment parameter to calculate MFP for twinning tau_0 = 1.0_pReal, & !< strength due to elements in solid solution L_tw = 1.0_pReal, & !< Length of twin nuclei in Burgers vectors L_tr = 1.0_pReal, & !< Length of trans nuclei in Burgers vectors - xc_twin = 1.0_pReal, & !< critical distance for formation of twin nucleus - xc_trans = 1.0_pReal, & !< critical distance for formation of trans nucleus + x_c_tw = 1.0_pReal, & !< critical distance for formation of twin nucleus + x_c_tr = 1.0_pReal, & !< critical distance for formation of trans nucleus V_cs = 1.0_pReal, & !< cross slip volume - sbResistance = 1.0_pReal, & !< value for shearband resistance - sbVelocity = 1.0_pReal, & !< value for shearband velocity_0 + xi_sb = 1.0_pReal, & !< value for shearband resistance + v_sb = 1.0_pReal, & !< value for shearband velocity_0 E_sb = 1.0_pReal, & !< activation energy for shear bands - SFE_0K = 1.0_pReal, & !< stacking fault energy at zero K - dSFE_dT = 1.0_pReal, & !< temperature dependence of stacking fault energy - gamma_fcc_hex = 1.0_pReal, & !< Free energy difference between austensite and martensite + Gamma_sf_0K = 1.0_pReal, & !< stacking fault energy at zero K + dGamma_sf_dT = 1.0_pReal, & !< temperature dependence of stacking fault energy + delta_G = 1.0_pReal, & !< Free energy difference between austensite and martensite i_tr = 1.0_pReal, & !< adjustment parameter to calculate MFP for transformation h = 1.0_pReal !< Stack height of hex nucleus real(pReal), allocatable, dimension(:) :: & - b_sl, & !< absolute length of burgers vector [m] for each slip system - b_tw, & !< absolute length of burgers vector [m] for each twin system - b_tr, & !< absolute length of burgers vector [m] for each transformation system - Delta_F,& !< activation energy for glide [J] for each slip system - v0, & !< dislocation velocity prefactor [m/s] for each slip system + b_sl, & !< absolute length of Burgers vector [m] for each slip system + b_tw, & !< absolute length of Burgers vector [m] for each twin system + b_tr, & !< absolute length of Burgers vector [m] for each transformation system + Q_s,& !< activation energy for glide [J] for each slip system + v_0, & !< dislocation velocity prefactor [m/s] for each slip system dot_N_0_tw, & !< twin nucleation rate [1/m³s] for each twin system dot_N_0_tr, & !< trans nucleation rate [1/m³s] for each trans system t_tw, & !< twin thickness [m] for each twin system - CLambdaSlip, & !< Adj. parameter for distance between 2 forest dislocations for each slip system + i_sl, & !< Adj. parameter for distance between 2 forest dislocations for each slip system t_tr, & !< martensite lamellar thickness [m] for each trans system and instance p, & !< p-exponent in glide velocity q, & !< q-exponent in glide velocity @@ -209,23 +209,23 @@ module function plastic_dislotwin_init() result(myPlasticity) rho_mob_0 = pl%get_asFloats('rho_mob_0', requiredSize=size(N_sl)) rho_dip_0 = pl%get_asFloats('rho_dip_0', requiredSize=size(N_sl)) - prm%v0 = pl%get_asFloats('v_0', requiredSize=size(N_sl)) - prm%b_sl = pl%get_asFloats('b_sl', requiredSize=size(N_sl)) - prm%Delta_F = pl%get_asFloats('Q_s', requiredSize=size(N_sl)) - prm%CLambdaSlip = pl%get_asFloats('i_sl', requiredSize=size(N_sl)) - prm%p = pl%get_asFloats('p_sl', requiredSize=size(N_sl)) - prm%q = pl%get_asFloats('q_sl', requiredSize=size(N_sl)) - prm%B = pl%get_asFloats('B', requiredSize=size(N_sl), & + prm%v_0 = pl%get_asFloats('v_0', requiredSize=size(N_sl)) + prm%b_sl = pl%get_asFloats('b_sl', requiredSize=size(N_sl)) + prm%Q_s = pl%get_asFloats('Q_s', requiredSize=size(N_sl)) + prm%i_sl = pl%get_asFloats('i_sl', requiredSize=size(N_sl)) + prm%p = pl%get_asFloats('p_sl', requiredSize=size(N_sl)) + prm%q = pl%get_asFloats('q_sl', requiredSize=size(N_sl)) + prm%B = pl%get_asFloats('B', requiredSize=size(N_sl), & defaultVal=[(0.0_pReal, i=1,size(N_sl))]) prm%tau_0 = pl%get_asFloat('tau_0') - prm%CEdgeDipMinDistance = pl%get_asFloat('D_a') - prm%D0 = pl%get_asFloat('D_0') - prm%Qsd = pl%get_asFloat('Q_cl') + prm%D_a = pl%get_asFloat('D_a') + prm%D_0 = pl%get_asFloat('D_0') + prm%Q_cl = pl%get_asFloat('Q_cl') prm%ExtendedDislocations = pl%get_asBool('extend_dislocations',defaultVal = .false.) if (prm%ExtendedDislocations) then - prm%SFE_0K = pl%get_asFloat('Gamma_sf_0K') - prm%dSFE_dT = pl%get_asFloat('dGamma_sf_dT') + prm%Gamma_sf_0K = pl%get_asFloat('Gamma_sf_0K') + prm%dGamma_sf_dT = pl%get_asFloat('dGamma_sf_dT') endif prm%dipoleformation = .not. pl%get_asBool('no_dipole_formation',defaultVal = .false.) @@ -238,29 +238,29 @@ module function plastic_dislotwin_init() result(myPlasticity) ! expand: family => system rho_mob_0 = math_expand(rho_mob_0, N_sl) rho_dip_0 = math_expand(rho_dip_0, N_sl) - prm%v0 = math_expand(prm%v0, N_sl) + prm%v_0 = math_expand(prm%v_0, N_sl) prm%b_sl = math_expand(prm%b_sl, N_sl) - prm%Delta_F = math_expand(prm%Delta_F, N_sl) - prm%CLambdaSlip = math_expand(prm%CLambdaSlip, N_sl) + prm%Q_s = math_expand(prm%Q_s, N_sl) + prm%i_sl = math_expand(prm%i_sl, N_sl) prm%p = math_expand(prm%p, N_sl) prm%q = math_expand(prm%q, N_sl) prm%B = math_expand(prm%B, N_sl) ! sanity checks - if ( prm%D0 <= 0.0_pReal) extmsg = trim(extmsg)//' D_0' - if ( prm%Qsd <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl' - if (any(rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_mob_0' - if (any(rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_dip_0' - if (any(prm%v0 < 0.0_pReal)) extmsg = trim(extmsg)//' v_0' - if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' b_sl' - if (any(prm%Delta_F <= 0.0_pReal)) extmsg = trim(extmsg)//' Q_s' - if (any(prm%CLambdaSlip <= 0.0_pReal)) extmsg = trim(extmsg)//' i_sl' - if (any(prm%B < 0.0_pReal)) extmsg = trim(extmsg)//' B' - if (any(prm%p<=0.0_pReal .or. prm%p>1.0_pReal)) extmsg = trim(extmsg)//' p_sl' - if (any(prm%q< 1.0_pReal .or. prm%q>2.0_pReal)) extmsg = trim(extmsg)//' q_sl' + if ( prm%D_0 <= 0.0_pReal) extmsg = trim(extmsg)//' D_0' + if ( prm%Q_cl <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl' + if (any(rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_mob_0' + if (any(rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_dip_0' + if (any(prm%v_0 < 0.0_pReal)) extmsg = trim(extmsg)//' v_0' + if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' b_sl' + if (any(prm%Q_s <= 0.0_pReal)) extmsg = trim(extmsg)//' Q_s' + if (any(prm%i_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' i_sl' + if (any(prm%B < 0.0_pReal)) extmsg = trim(extmsg)//' B' + if (any(prm%p<=0.0_pReal .or. prm%p>1.0_pReal)) extmsg = trim(extmsg)//' p_sl' + if (any(prm%q< 1.0_pReal .or. prm%q>2.0_pReal)) extmsg = trim(extmsg)//' q_sl' else slipActive rho_mob_0 = emptyRealArray; rho_dip_0 = emptyRealArray - allocate(prm%b_sl,prm%Delta_F,prm%v0,prm%CLambdaSlip,prm%p,prm%q,prm%B,source=emptyRealArray) + allocate(prm%b_sl,prm%Q_s,prm%v_0,prm%i_sl,prm%p,prm%q,prm%B,source=emptyRealArray) allocate(prm%forestProjection(0,0),prm%h_sl_sl(0,0)) endif slipActive @@ -279,7 +279,7 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%t_tw = pl%get_asFloats('t_tw', requiredSize=size(N_tw)) prm%r = pl%get_asFloats('p_tw', requiredSize=size(N_tw)) - prm%xc_twin = pl%get_asFloat('x_c_tw') + prm%x_c_tw = pl%get_asFloat('x_c_tw') prm%L_tw = pl%get_asFloat('L_tw') prm%i_tw = pl%get_asFloat('i_tw') @@ -300,7 +300,7 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%r = math_expand(prm%r,N_tw) ! sanity checks - if ( prm%xc_twin < 0.0_pReal) extmsg = trim(extmsg)//' x_c_twin' + if ( prm%x_c_tw < 0.0_pReal) extmsg = trim(extmsg)//' x_c_twin' if ( prm%L_tw < 0.0_pReal) extmsg = trim(extmsg)//' L_tw' if ( prm%i_tw < 0.0_pReal) extmsg = trim(extmsg)//' i_tw' if (any(prm%b_tw < 0.0_pReal)) extmsg = trim(extmsg)//' b_tw' @@ -324,8 +324,8 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%h = pl%get_asFloat('h', defaultVal=0.0_pReal) ! ToDo: How to handle that??? prm%i_tr = pl%get_asFloat('i_tr', defaultVal=0.0_pReal) ! ToDo: How to handle that??? - prm%gamma_fcc_hex = pl%get_asFloat('delta_G') - prm%xc_trans = pl%get_asFloat('x_c_tr', defaultVal=0.0_pReal) ! ToDo: How to handle that??? + prm%delta_G = pl%get_asFloat('delta_G') + prm%x_c_tr = pl%get_asFloat('x_c_tr', defaultVal=0.0_pReal) ! ToDo: How to handle that??? prm%L_tr = pl%get_asFloat('L_tr') prm%h_tr_tr = lattice_interaction_TransByTrans(N_tr,pl%get_asFloats('h_tr_tr'), & @@ -351,7 +351,7 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%s = math_expand(prm%s,N_tr) ! sanity checks - if ( prm%xc_trans < 0.0_pReal) extmsg = trim(extmsg)//' x_c_trans' + if ( prm%x_c_tr < 0.0_pReal) extmsg = trim(extmsg)//' x_c_trans' if ( prm%L_tr < 0.0_pReal) extmsg = trim(extmsg)//' L_tr' if ( prm%i_tr < 0.0_pReal) extmsg = trim(extmsg)//' i_tr' if (any(prm%t_tr < 0.0_pReal)) extmsg = trim(extmsg)//' t_tr' @@ -366,15 +366,15 @@ module function plastic_dislotwin_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! shearband related parameters - prm%sbVelocity = pl%get_asFloat('v_sb',defaultVal=0.0_pReal) - if (prm%sbVelocity > 0.0_pReal) then - prm%sbResistance = pl%get_asFloat('xi_sb') + prm%v_sb = pl%get_asFloat('v_sb',defaultVal=0.0_pReal) + if (prm%v_sb > 0.0_pReal) then + prm%xi_sb = pl%get_asFloat('xi_sb') prm%E_sb = pl%get_asFloat('Q_sb') prm%p_sb = pl%get_asFloat('p_sb') prm%q_sb = pl%get_asFloat('q_sb') ! sanity checks - if (prm%sbResistance < 0.0_pReal) extmsg = trim(extmsg)//' xi_sb' + if (prm%xi_sb < 0.0_pReal) extmsg = trim(extmsg)//' xi_sb' if (prm%E_sb < 0.0_pReal) extmsg = trim(extmsg)//' Q_sb' if (prm%p_sb <= 0.0_pReal) extmsg = trim(extmsg)//' p_sb' if (prm%q_sb <= 0.0_pReal) extmsg = trim(extmsg)//' q_sb' @@ -386,8 +386,8 @@ module function plastic_dislotwin_init() result(myPlasticity) prm%D = pl%get_asFloat('D') twinOrSlipActive: if (prm%sum_N_tw + prm%sum_N_tr > 0) then - prm%SFE_0K = pl%get_asFloat('Gamma_sf_0K') - prm%dSFE_dT = pl%get_asFloat('dGamma_sf_dT') + prm%Gamma_sf_0K = pl%get_asFloat('Gamma_sf_0K') + prm%dGamma_sf_dT = pl%get_asFloat('dGamma_sf_dT') prm%V_cs = pl%get_asFloat('V_cs') endif twinOrSlipActive @@ -602,7 +602,7 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) Lp = Lp * f_unrotated dLp_dMp = dLp_dMp * f_unrotated - shearBandingContribution: if(dNeq0(prm%sbVelocity)) then + shearBandingContribution: if(dNeq0(prm%v_sb)) then BoltzmannRatio = prm%E_sb/(kB*T) call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design? @@ -613,10 +613,10 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) tau = math_tensordot(Mp,P_sb) significantShearBandStress: if (abs(tau) > tol_math_check) then - StressRatio_p = (abs(tau)/prm%sbResistance)**prm%p_sb - dot_gamma_sb = sign(prm%sbVelocity*exp(-BoltzmannRatio*(1-StressRatio_p)**prm%q_sb), tau) - ddot_gamma_dtau = abs(dot_gamma_sb)*BoltzmannRatio* prm%p_sb*prm%q_sb/ prm%sbResistance & - * (abs(tau)/prm%sbResistance)**(prm%p_sb-1.0_pReal) & + StressRatio_p = (abs(tau)/prm%xi_sb)**prm%p_sb + dot_gamma_sb = sign(prm%v_sb*exp(-BoltzmannRatio*(1-StressRatio_p)**prm%q_sb), tau) + ddot_gamma_dtau = abs(dot_gamma_sb)*BoltzmannRatio* prm%p_sb*prm%q_sb/ prm%xi_sb & + * (abs(tau)/prm%xi_sb)**(prm%p_sb-1.0_pReal) & * (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal) Lp = Lp + dot_gamma_sb * P_sb @@ -654,7 +654,7 @@ module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) Gamma, & !< stacking fault energy tau, & sigma_cl, & !< climb stress - b_d !< ratio of burgers vector to stacking fault width + b_d !< ratio of Burgers vector to stacking fault width real(pReal), dimension(param(instance)%sum_N_sl) :: & dot_rho_dip_formation, & dot_rho_dip_climb, & @@ -675,7 +675,7 @@ module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) call kinetics_slip(Mp,T,instance,of,dot_gamma_sl) dot%gamma_sl(:,of) = abs(dot_gamma_sl) - rho_dip_distance_min = prm%CEdgeDipMinDistance*prm%b_sl + rho_dip_distance_min = prm%D_a*prm%b_sl slipState: do i = 1, prm%sum_N_sl tau = math_tensordot(Mp,prm%P_sl(1:3,1:3,i)) @@ -701,12 +701,12 @@ module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) !@details: Refer: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 sigma_cl = dot_product(prm%n0_sl(1:3,i),matmul(Mp,prm%n0_sl(1:3,i))) if (prm%ExtendedDislocations) then - Gamma = prm%SFE_0K + prm%dSFE_dT * T + Gamma = prm%Gamma_sf_0K + prm%dGamma_sf_dT * T b_d = 24.0_pReal*PI*(1.0_pReal - prm%nu)/(2.0_pReal + prm%nu)* Gamma/(prm%mu*prm%b_sl(i)) else b_d = 1.0_pReal endif - v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Qsd/(kB*T)) & + v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Q_cl/(kB*T)) & * (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal) dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,of) & @@ -768,7 +768,7 @@ module subroutine plastic_dislotwin_dependentState(T,instance,of) sumf_twin = sum(stt%f_tw(1:prm%sum_N_tw,of)) sumf_trans = sum(stt%f_tr(1:prm%sum_N_tr,of)) - Gamma = prm%SFE_0K + prm%dSFE_dT * T + Gamma = prm%Gamma_sf_0K + prm%dGamma_sf_dT * T !* rescaled volume fraction for topology f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,of)/prm%t_tw ! this is per system ... @@ -776,7 +776,7 @@ module subroutine plastic_dislotwin_dependentState(T,instance,of) ! ToDo ...Physically correct, but naming could be adjusted inv_lambda_sl_sl = sqrt(matmul(prm%forestProjection, & - stt%rho_mob(:,of)+stt%rho_dip(:,of)))/prm%CLambdaSlip + stt%rho_mob(:,of)+stt%rho_dip(:,of)))/prm%i_sl if (prm%sum_N_tw > 0 .and. prm%sum_N_sl > 0) & inv_lambda_sl_tw = matmul(prm%h_sl_tw,f_over_t_tw)/(1.0_pReal-sumf_twin) @@ -805,21 +805,21 @@ module subroutine plastic_dislotwin_dependentState(T,instance,of) !* threshold stress for growing twin/martensite if(prm%sum_N_tw == prm%sum_N_sl) & dst%tau_hat_tw(:,of) = Gamma/(3.0_pReal*prm%b_tw) & - + 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_sl) ! slip burgers here correct? + + 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_sl) ! slip Burgers here correct? if(prm%sum_N_tr == prm%sum_N_sl) & dst%tau_hat_tr(:,of) = Gamma/(3.0_pReal*prm%b_tr) & - + 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_sl) & ! slip burgers here correct? - + prm%h*prm%gamma_fcc_hex/ (3.0_pReal*prm%b_tr) + + 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_sl) & ! slip Burgers here correct? + + prm%h*prm%delta_G/ (3.0_pReal*prm%b_tr) dst%V_tw(:,of) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,of)**2.0_pReal dst%V_tr(:,of) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,of)**2.0_pReal - x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the burgers vector for slip and is the same for twin and trans - dst%tau_r_tw(:,of) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%xc_twin)+cos(pi/3.0_pReal)/x0) + x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip and is the same for twin and trans + dst%tau_r_tw(:,of) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0) - x0 = prm%mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the burgers vector for slip - dst%tau_r_tr(:,of) = prm%mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%xc_trans)+cos(pi/3.0_pReal)/x0) + x0 = prm%mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip + dst%tau_r_tr(:,of) = prm%mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(pi/3.0_pReal)/x0) end associate @@ -927,8 +927,8 @@ pure subroutine kinetics_slip(Mp,T,instance,of, & significantStress: where(tau_eff > tol_math_check) stressRatio = tau_eff/prm%tau_0 StressRatio_p = stressRatio** prm%p - BoltzmannRatio = prm%Delta_F/(kB*T) - v_wait_inverse = prm%v0**(-1.0_pReal) * exp(BoltzmannRatio*(1.0_pReal-StressRatio_p)** prm%q) + BoltzmannRatio = prm%Q_s/(kB*T) + v_wait_inverse = prm%v_0**(-1.0_pReal) * exp(BoltzmannRatio*(1.0_pReal-StressRatio_p)** prm%q) v_run_inverse = prm%B/(tau_eff*prm%b_sl) dot_gamma_sl = sign(stt%rho_mob(:,of)*prm%b_sl/(v_wait_inverse+v_run_inverse),tau) diff --git a/src/constitutive_plastic_isotropic.f90 b/src/constitutive_plastic_isotropic.f90 index 477938145..2c9028671 100644 --- a/src/constitutive_plastic_isotropic.f90 +++ b/src/constitutive_plastic_isotropic.f90 @@ -14,7 +14,7 @@ submodule(constitutive:constitutive_plastic) plastic_isotropic M, & !< Taylor factor dot_gamma_0, & !< reference strain rate n, & !< stress exponent - h0, & + h_0, & h_ln, & xi_inf, & !< maximum critical stress a, & @@ -109,7 +109,7 @@ module function plastic_isotropic_init() result(myPlasticity) prm%xi_inf = pl%get_asFloat('xi_inf') prm%dot_gamma_0 = pl%get_asFloat('dot_gamma_0') prm%n = pl%get_asFloat('n') - prm%h0 = pl%get_asFloat('h_0') + prm%h_0 = pl%get_asFloat('h_0') prm%M = pl%get_asFloat('M') prm%h_ln = pl%get_asFloat('h_ln', defaultVal=0.0_pReal) prm%c_1 = pl%get_asFloat('c_1', defaultVal=0.0_pReal) @@ -310,7 +310,7 @@ module subroutine plastic_isotropic_dotState(Mp,instance,of) / prm%c_4 * (dot_gamma / prm%dot_gamma_0)**(1.0_pReal / prm%n) endif dot%xi(of) = dot_gamma & - * ( prm%h0 + prm%h_ln * log(dot_gamma) ) & + * ( prm%h_0 + prm%h_ln * log(dot_gamma) ) & * abs( 1.0_pReal - stt%xi(of)/xi_inf_star )**prm%a & * sign(1.0_pReal, 1.0_pReal - stt%xi(of)/xi_inf_star) else diff --git a/src/constitutive_plastic_kinehardening.f90 b/src/constitutive_plastic_kinehardening.f90 index 55e482db6..777691242 100644 --- a/src/constitutive_plastic_kinehardening.f90 +++ b/src/constitutive_plastic_kinehardening.f90 @@ -9,15 +9,15 @@ submodule(constitutive:constitutive_plastic) plastic_kinehardening type :: tParameters real(pReal) :: & - gdot0 = 1.0_pReal, & !< reference shear strain rate for slip - n = 1.0_pReal !< stress exponent for slip + n = 1.0_pReal, & !< stress exponent for slip + dot_gamma_0 = 1.0_pReal !< reference shear strain rate for slip real(pReal), allocatable, dimension(:) :: & - theta0, & !< initial hardening rate of forward stress for each slip - theta1, & !< asymptotic hardening rate of forward stress for each slip - theta0_b, & !< initial hardening rate of back stress for each slip - theta1_b, & !< asymptotic hardening rate of back stress for each slip - tau1, & - tau1_b + h_0_f, & !< initial hardening rate of forward stress for each slip + h_inf_f, & !< asymptotic hardening rate of forward stress for each slip + h_0_b, & !< initial hardening rate of back stress for each slip + h_inf_b, & !< asymptotic hardening rate of back stress for each slip + xi_inf_f, & + xi_inf_b real(pReal), allocatable, dimension(:,:) :: & interaction_slipslip !< slip resistance from slip activity real(pReal), allocatable, dimension(:,:,:) :: & @@ -125,7 +125,7 @@ module function plastic_kinehardening_init() result(myPlasticity) phase%get_asFloat('c/a',defaultVal=0.0_pReal)) if(trim(phase%get_asString('lattice')) == 'bcc') then - a = pl%get_asFloats('nonSchmid_coefficients',defaultVal = emptyRealArray) + a = pl%get_asFloats('a_nonSchmid',defaultVal = emptyRealArray) if(size(a) > 0) prm%nonSchmidActive = .true. prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,a,+1) prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,a,-1) @@ -137,38 +137,38 @@ module function plastic_kinehardening_init() result(myPlasticity) pl%get_asFloats('h_sl_sl'), & phase%get_asString('lattice')) - xi_0 = pl%get_asFloats('xi_0', requiredSize=size(N_sl)) - prm%tau1 = pl%get_asFloats('xi_inf_f', requiredSize=size(N_sl)) - prm%tau1_b = pl%get_asFloats('xi_inf_b', requiredSize=size(N_sl)) - prm%theta0 = pl%get_asFloats('h_0_f', requiredSize=size(N_sl)) - prm%theta1 = pl%get_asFloats('h_inf_f', requiredSize=size(N_sl)) - prm%theta0_b = pl%get_asFloats('h_0_b', requiredSize=size(N_sl)) - prm%theta1_b = pl%get_asFloats('h_inf_b', requiredSize=size(N_sl)) + xi_0 = pl%get_asFloats('xi_0', requiredSize=size(N_sl)) + prm%xi_inf_f = pl%get_asFloats('xi_inf_f', requiredSize=size(N_sl)) + prm%xi_inf_b = pl%get_asFloats('xi_inf_b', requiredSize=size(N_sl)) + prm%h_0_f = pl%get_asFloats('h_0_f', requiredSize=size(N_sl)) + prm%h_inf_f = pl%get_asFloats('h_inf_f', requiredSize=size(N_sl)) + prm%h_0_b = pl%get_asFloats('h_0_b', requiredSize=size(N_sl)) + prm%h_inf_b = pl%get_asFloats('h_inf_b', requiredSize=size(N_sl)) - prm%gdot0 = pl%get_asFloat('dot_gamma_0') - prm%n = pl%get_asFloat('n') + prm%dot_gamma_0 = pl%get_asFloat('dot_gamma_0') + prm%n = pl%get_asFloat('n') ! expand: family => system - xi_0 = math_expand(xi_0, N_sl) - prm%tau1 = math_expand(prm%tau1, N_sl) - prm%tau1_b = math_expand(prm%tau1_b, N_sl) - prm%theta0 = math_expand(prm%theta0, N_sl) - prm%theta1 = math_expand(prm%theta1, N_sl) - prm%theta0_b = math_expand(prm%theta0_b,N_sl) - prm%theta1_b = math_expand(prm%theta1_b,N_sl) + xi_0 = math_expand(xi_0, N_sl) + prm%xi_inf_f = math_expand(prm%xi_inf_f, N_sl) + prm%xi_inf_b = math_expand(prm%xi_inf_b, N_sl) + prm%h_0_f = math_expand(prm%h_0_f, N_sl) + prm%h_inf_f = math_expand(prm%h_inf_f, N_sl) + prm%h_0_b = math_expand(prm%h_0_b, N_sl) + prm%h_inf_b = math_expand(prm%h_inf_b, N_sl) !-------------------------------------------------------------------------------------------------- ! sanity checks - if ( prm%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' dot_gamma_0' - if ( prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' n' - if (any(xi_0 <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_0' - if (any(prm%tau1 <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_inf_f' - if (any(prm%tau1_b <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_inf_b' + if ( prm%dot_gamma_0 <= 0.0_pReal) extmsg = trim(extmsg)//' dot_gamma_0' + if ( prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' n' + if (any(xi_0 <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_0' + if (any(prm%xi_inf_f <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_inf_f' + if (any(prm%xi_inf_b <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_inf_b' !ToDo: Any sensible checks for theta? else slipActive xi_0 = emptyRealArray - allocate(prm%tau1,prm%tau1_b,prm%theta0,prm%theta1,prm%theta0_b,prm%theta1_b,source=emptyRealArray) + allocate(prm%xi_inf_f,prm%xi_inf_b,prm%h_0_f,prm%h_inf_f,prm%h_0_b,prm%h_inf_b,source=emptyRealArray) allocate(prm%interaction_SlipSlip(0,0)) endif slipActive @@ -303,16 +303,16 @@ module subroutine plastic_kinehardening_dotState(Mp,instance,of) dot%crss(:,of) = matmul(prm%interaction_SlipSlip,dot%accshear(:,of)) & - * ( prm%theta1 & - + (prm%theta0 - prm%theta1 + prm%theta0*prm%theta1*sumGamma/prm%tau1) & - * exp(-sumGamma*prm%theta0/prm%tau1) & + * ( prm%h_inf_f & + + (prm%h_0_f - prm%h_inf_f + prm%h_0_f*prm%h_inf_f*sumGamma/prm%xi_inf_f) & + * exp(-sumGamma*prm%h_0_f/prm%xi_inf_f) & ) dot%crss_back(:,of) = stt%sense(:,of)*dot%accshear(:,of) * & - ( prm%theta1_b + & - (prm%theta0_b - prm%theta1_b & - + prm%theta0_b*prm%theta1_b/(prm%tau1_b+stt%chi0(:,of))*(stt%accshear(:,of)-stt%gamma0(:,of))& - ) *exp(-(stt%accshear(:,of)-stt%gamma0(:,of)) *prm%theta0_b/(prm%tau1_b+stt%chi0(:,of))) & + ( prm%h_inf_b + & + (prm%h_0_b - prm%h_inf_b & + + prm%h_0_b*prm%h_inf_b/(prm%xi_inf_b+stt%chi0(:,of))*(stt%accshear(:,of)-stt%gamma0(:,of))& + ) *exp(-(stt%accshear(:,of)-stt%gamma0(:,of)) *prm%h_0_b/(prm%xi_inf_b+stt%chi0(:,of))) & ) end associate @@ -442,14 +442,14 @@ pure subroutine kinetics(Mp,instance,of, & enddo where(dNeq0(tau_pos)) - gdot_pos = prm%gdot0 * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active + gdot_pos = prm%dot_gamma_0 * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active * sign(abs(tau_pos/stt%crss(:,of))**prm%n, tau_pos) else where gdot_pos = 0.0_pReal end where where(dNeq0(tau_neg)) - gdot_neg = prm%gdot0 * 0.5_pReal & ! only used if non-Schmid active, always 1/2 + gdot_neg = prm%dot_gamma_0 * 0.5_pReal & ! only used if non-Schmid active, always 1/2 * sign(abs(tau_neg/stt%crss(:,of))**prm%n, tau_neg) else where gdot_neg = 0.0_pReal diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index 21523161b..8238f17c9 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -46,58 +46,58 @@ submodule(constitutive:constitutive_plastic) plastic_nonlocal type :: tInitialParameters !< container type for internal constitutive parameters real(pReal) :: & - rhoSglScatter, & !< standard deviation of scatter in initial dislocation density - rhoSglRandom, & - rhoSglRandomBinning + sigma_rho_u, & !< standard deviation of scatter in initial dislocation density + random_rho_u, & + random_rho_u_binning real(pReal), dimension(:), allocatable :: & - rhoSglEdgePos0, & !< initial edge_pos dislocation density - rhoSglEdgeNeg0, & !< initial edge_neg dislocation density - rhoSglScrewPos0, & !< initial screw_pos dislocation density - rhoSglScrewNeg0, & !< initial screw_neg dislocation density - rhoDipEdge0, & !< initial edge dipole dislocation density - rhoDipScrew0 !< initial screw dipole dislocation density + rho_u_ed_pos_0, & !< initial edge_pos dislocation density + rho_u_ed_neg_0, & !< initial edge_neg dislocation density + rho_u_sc_pos_0, & !< initial screw_pos dislocation density + rho_u_sc_neg_0, & !< initial screw_neg dislocation density + rho_d_ed_0, & !< initial edge dipole dislocation density + rho_d_sc_0 !< initial screw dipole dislocation density integer, dimension(:), allocatable :: & N_sl end type tInitialParameters type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & - atomicVolume, & !< atomic volume - Dsd0, & !< prefactor for self-diffusion coefficient - selfDiffusionEnergy, & !< activation enthalpy for diffusion + V_at, & !< atomic volume + D_0, & !< prefactor for self-diffusion coefficient + Q_cl, & !< activation enthalpy for diffusion atol_rho, & !< absolute tolerance for dislocation density in state integration - significantRho, & !< density considered significant - significantN, & !< number of dislocations considered significant - doublekinkwidth, & !< width of a doubkle kink in multiples of the burgers vector length b - solidSolutionEnergy, & !< activation energy for solid solution in J - solidSolutionSize, & !< solid solution obstacle size in multiples of the burgers vector length - solidSolutionConcentration, & !< concentration of solid solution in atomic parts + rho_significant, & !< density considered significant + rho_min, & !< number of dislocations considered significant + w, & !< width of a doubkle kink in multiples of the Burgers vector length b + Q_sol, & !< activation energy for solid solution in J + f_sol, & !< solid solution obstacle size in multiples of the Burgers vector length + c_sol, & !< concentration of solid solution in atomic parts p, & !< parameter for kinetic law (Kocks,Argon,Ashby) q, & !< parameter for kinetic law (Kocks,Argon,Ashby) - viscosity, & !< viscosity for dislocation glide in Pa s - fattack, & !< attack frequency in Hz - surfaceTransmissivity, & !< transmissivity at free surface - grainboundaryTransmissivity, & !< transmissivity at grain boundary (identified by different texture) - CFLfactor, & !< safety factor for CFL flux condition - fEdgeMultiplication, & !< factor that determines how much edge dislocations contribute to multiplication (0...1) - linetensionEffect, & - edgeJogFactor, & + eta, & !< viscosity for dislocation glide in Pa s + nu_a, & !< attack frequency in Hz + chi_surface, & !< transmissivity at free surface + chi_GB, & !< transmissivity at grain boundary (identified by different texture) + f_c, & !< safety factor for CFL flux condition + f_ed_mult, & !< factor that determines how much edge dislocations contribute to multiplication (0...1) + f_F, & + f_ed, & mu, & nu real(pReal), dimension(:), allocatable :: & - minDipoleHeight_edge, & !< minimum stable edge dipole height - minDipoleHeight_screw, & !< minimum stable screw dipole height - peierlsstress_edge, & - peierlsstress_screw, & - lambda0, & !< mean free path prefactor for each - burgers !< absolute length of burgers vector [m] + d_ed, & !< minimum stable edge dipole height + d_sc, & !< minimum stable screw dipole height + tau_Peierls_ed, & + tau_Peierls_sc, & + i_sl, & !< mean free path prefactor for each + b_sl !< absolute length of Burgers vector [m] real(pReal), dimension(:,:), allocatable :: & slip_normal, & slip_direction, & slip_transverse, & minDipoleHeight, & ! edge and screw peierlsstress, & ! edge and screw - interactionSlipSlip ,& !< coefficients for slip-slip interaction + h_sl_sl ,& !< coefficients for slip-slip interaction forestProjection_Edge, & !< matrix of forest projections of edge dislocations forestProjection_Screw !< matrix of forest projections of screw dislocations real(pReal), dimension(:,:,:), allocatable :: & @@ -244,7 +244,7 @@ module function plastic_nonlocal_init() result(myPlasticity) phase%get_asFloat('c/a',defaultVal=0.0_pReal)) if(trim(phase%get_asString('lattice')) == 'bcc') then - a = pl%get_asFloats('nonSchmid_coefficients',defaultVal = emptyRealArray) + a = pl%get_asFloats('a_nonSchmid',defaultVal = emptyRealArray) if(size(a) > 0) prm%nonSchmidActive = .true. prm%nonSchmid_pos = lattice_nonSchmidMatrix(ini%N_sl,a,+1) prm%nonSchmid_neg = lattice_nonSchmidMatrix(ini%N_sl,a,-1) @@ -253,9 +253,9 @@ module function plastic_nonlocal_init() result(myPlasticity) prm%nonSchmid_neg = prm%Schmid endif - prm%interactionSlipSlip = lattice_interaction_SlipBySlip(ini%N_sl, & - pl%get_asFloats('h_sl_sl'), & - phase%get_asString('lattice')) + prm%h_sl_sl = lattice_interaction_SlipBySlip(ini%N_sl, & + pl%get_asFloats('h_sl_sl'), & + phase%get_asString('lattice')) prm%forestProjection_edge = lattice_forestProjection_edge (ini%N_sl,phase%get_asString('lattice'),& phase%get_asFloat('c/a',defaultVal=0.0_pReal)) @@ -279,113 +279,113 @@ module function plastic_nonlocal_init() result(myPlasticity) enddo enddo - ini%rhoSglEdgePos0 = pl%get_asFloats('rho_u_ed_pos_0', requiredSize=size(ini%N_sl)) - ini%rhoSglEdgeNeg0 = pl%get_asFloats('rho_u_ed_neg_0', requiredSize=size(ini%N_sl)) - ini%rhoSglScrewPos0 = pl%get_asFloats('rho_u_sc_pos_0', requiredSize=size(ini%N_sl)) - ini%rhoSglScrewNeg0 = pl%get_asFloats('rho_u_sc_neg_0', requiredSize=size(ini%N_sl)) - ini%rhoDipEdge0 = pl%get_asFloats('rho_d_ed_0', requiredSize=size(ini%N_sl)) - ini%rhoDipScrew0 = pl%get_asFloats('rho_d_sc_0', requiredSize=size(ini%N_sl)) + ini%rho_u_ed_pos_0 = pl%get_asFloats('rho_u_ed_pos_0', requiredSize=size(ini%N_sl)) + ini%rho_u_ed_neg_0 = pl%get_asFloats('rho_u_ed_neg_0', requiredSize=size(ini%N_sl)) + ini%rho_u_sc_pos_0 = pl%get_asFloats('rho_u_sc_pos_0', requiredSize=size(ini%N_sl)) + ini%rho_u_sc_neg_0 = pl%get_asFloats('rho_u_sc_neg_0', requiredSize=size(ini%N_sl)) + ini%rho_d_ed_0 = pl%get_asFloats('rho_d_ed_0', requiredSize=size(ini%N_sl)) + ini%rho_d_sc_0 = pl%get_asFloats('rho_d_sc_0', requiredSize=size(ini%N_sl)) - prm%lambda0 = pl%get_asFloats('i_sl', requiredSize=size(ini%N_sl)) - prm%burgers = pl%get_asFloats('b_sl', requiredSize=size(ini%N_sl)) + prm%i_sl = pl%get_asFloats('i_sl', requiredSize=size(ini%N_sl)) + prm%b_sl = pl%get_asFloats('b_sl', requiredSize=size(ini%N_sl)) - prm%lambda0 = math_expand(prm%lambda0,ini%N_sl) - prm%burgers = math_expand(prm%burgers,ini%N_sl) + prm%i_sl = math_expand(prm%i_sl,ini%N_sl) + prm%b_sl = math_expand(prm%b_sl,ini%N_sl) - prm%minDipoleHeight_edge = pl%get_asFloats('d_ed', requiredSize=size(ini%N_sl)) - prm%minDipoleHeight_screw = pl%get_asFloats('d_sc', requiredSize=size(ini%N_sl)) - prm%minDipoleHeight_edge = math_expand(prm%minDipoleHeight_edge, ini%N_sl) - prm%minDipoleHeight_screw = math_expand(prm%minDipoleHeight_screw,ini%N_sl) + prm%d_ed = pl%get_asFloats('d_ed', requiredSize=size(ini%N_sl)) + prm%d_sc = pl%get_asFloats('d_sc', requiredSize=size(ini%N_sl)) + prm%d_ed = math_expand(prm%d_ed,ini%N_sl) + prm%d_sc = math_expand(prm%d_sc,ini%N_sl) allocate(prm%minDipoleHeight(prm%sum_N_sl,2)) - prm%minDipoleHeight(:,1) = prm%minDipoleHeight_edge - prm%minDipoleHeight(:,2) = prm%minDipoleHeight_screw + prm%minDipoleHeight(:,1) = prm%d_ed + prm%minDipoleHeight(:,2) = prm%d_sc - prm%peierlsstress_edge = pl%get_asFloats('tau_peierls_ed', requiredSize=size(ini%N_sl)) - prm%peierlsstress_screw = pl%get_asFloats('tau_peierls_sc', requiredSize=size(ini%N_sl)) - prm%peierlsstress_edge = math_expand(prm%peierlsstress_edge, ini%N_sl) - prm%peierlsstress_screw = math_expand(prm%peierlsstress_screw,ini%N_sl) + prm%tau_Peierls_ed = pl%get_asFloats('tau_Peierls_ed', requiredSize=size(ini%N_sl)) + prm%tau_Peierls_sc = pl%get_asFloats('tau_Peierls_sc', requiredSize=size(ini%N_sl)) + prm%tau_Peierls_ed = math_expand(prm%tau_Peierls_ed,ini%N_sl) + prm%tau_Peierls_sc = math_expand(prm%tau_Peierls_sc,ini%N_sl) allocate(prm%peierlsstress(prm%sum_N_sl,2)) - prm%peierlsstress(:,1) = prm%peierlsstress_edge - prm%peierlsstress(:,2) = prm%peierlsstress_screw + prm%peierlsstress(:,1) = prm%tau_Peierls_ed + prm%peierlsstress(:,2) = prm%tau_Peierls_sc - prm%significantRho = pl%get_asFloat('rho_significant') - prm%significantN = pl%get_asFloat('rho_num_significant', 0.0_pReal) - prm%CFLfactor = pl%get_asFloat('f_c',defaultVal=2.0_pReal) + prm%rho_significant = pl%get_asFloat('rho_significant') + prm%rho_min = pl%get_asFloat('rho_min', 0.0_pReal) + prm%f_c = pl%get_asFloat('f_c',defaultVal=2.0_pReal) - prm%atomicVolume = pl%get_asFloat('V_at') - prm%Dsd0 = pl%get_asFloat('D_0') !,'dsd0' - prm%selfDiffusionEnergy = pl%get_asFloat('Q_cl') !,'qsd' - prm%linetensionEffect = pl%get_asFloat('f_F') - prm%edgeJogFactor = pl%get_asFloat('f_ed') !,'edgejogs' - prm%doublekinkwidth = pl%get_asFloat('w') - prm%solidSolutionEnergy = pl%get_asFloat('Q_sol') - prm%solidSolutionSize = pl%get_asFloat('f_sol') - prm%solidSolutionConcentration = pl%get_asFloat('c_sol') + prm%V_at = pl%get_asFloat('V_at') + prm%D_0 = pl%get_asFloat('D_0') + prm%Q_cl = pl%get_asFloat('Q_cl') + prm%f_F = pl%get_asFloat('f_F') + prm%f_ed = pl%get_asFloat('f_ed') !,'edgejogs' + prm%w = pl%get_asFloat('w') + prm%Q_sol = pl%get_asFloat('Q_sol') + prm%f_sol = pl%get_asFloat('f_sol') + prm%c_sol = pl%get_asFloat('c_sol') - prm%p = pl%get_asFloat('p_sl') - prm%q = pl%get_asFloat('q_sl') - prm%viscosity = pl%get_asFloat('eta') - prm%fattack = pl%get_asFloat('nu_a') + prm%p = pl%get_asFloat('p_sl') + prm%q = pl%get_asFloat('q_sl') + prm%eta = pl%get_asFloat('eta') + prm%nu_a = pl%get_asFloat('nu_a') ! ToDo: discuss logic - ini%rhoSglScatter = pl%get_asFloat('sigma_rho_u') - ini%rhoSglRandom = pl%get_asFloat('random_rho_u',defaultVal= 0.0_pReal) + ini%sigma_rho_u = pl%get_asFloat('sigma_rho_u') + ini%random_rho_u = pl%get_asFloat('random_rho_u',defaultVal= 0.0_pReal) if (pl%contains('random_rho_u')) & - ini%rhoSglRandomBinning = pl%get_asFloat('random_rho_u_binning',defaultVal=0.0_pReal) !ToDo: useful default? + ini%random_rho_u_binning = pl%get_asFloat('random_rho_u_binning',defaultVal=0.0_pReal) !ToDo: useful default? ! if (rhoSglRandom(instance) < 0.0_pReal) & ! if (rhoSglRandomBinning(instance) <= 0.0_pReal) & - prm%surfaceTransmissivity = pl%get_asFloat('chi_surface',defaultVal=1.0_pReal) - prm%grainboundaryTransmissivity = pl%get_asFloat('chi_GB', defaultVal=-1.0_pReal) - prm%fEdgeMultiplication = pl%get_asFloat('f_ed_mult') + prm%chi_surface = pl%get_asFloat('chi_surface',defaultVal=1.0_pReal) + prm%chi_GB = pl%get_asFloat('chi_GB', defaultVal=-1.0_pReal) + prm%f_ed_mult = pl%get_asFloat('f_ed_mult') prm%shortRangeStressCorrection = pl%get_asBool('short_range_stress_correction', defaultVal = .false.) !-------------------------------------------------------------------------------------------------- ! sanity checks - if (any(prm%burgers < 0.0_pReal)) extmsg = trim(extmsg)//' b_sl' - if (any(prm%lambda0 <= 0.0_pReal)) extmsg = trim(extmsg)//' i_sl' + if (any(prm%b_sl < 0.0_pReal)) extmsg = trim(extmsg)//' b_sl' + if (any(prm%i_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' i_sl' - if (any(ini%rhoSglEdgePos0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_ed_pos_0' - if (any(ini%rhoSglEdgeNeg0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_ed_neg_0' - if (any(ini%rhoSglScrewPos0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_sc_pos_0' - if (any(ini%rhoSglScrewNeg0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_sc_neg_0' - if (any(ini%rhoDipEdge0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_d_ed_0' - if (any(ini%rhoDipScrew0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_d_sc_0' + if (any(ini%rho_u_ed_pos_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_ed_pos_0' + if (any(ini%rho_u_ed_neg_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_ed_neg_0' + if (any(ini%rho_u_sc_pos_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_sc_pos_0' + if (any(ini%rho_u_sc_neg_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_sc_neg_0' + if (any(ini%rho_d_ed_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_d_ed_0' + if (any(ini%rho_d_sc_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_d_sc_0' if (any(prm%peierlsstress < 0.0_pReal)) extmsg = trim(extmsg)//' tau_peierls' if (any(prm%minDipoleHeight < 0.0_pReal)) extmsg = trim(extmsg)//' d_ed or d_sc' - if (prm%viscosity <= 0.0_pReal) extmsg = trim(extmsg)//' eta' - if (prm%selfDiffusionEnergy <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl' - if (prm%fattack <= 0.0_pReal) extmsg = trim(extmsg)//' nu_a' - if (prm%doublekinkwidth <= 0.0_pReal) extmsg = trim(extmsg)//' w' - if (prm%Dsd0 < 0.0_pReal) extmsg = trim(extmsg)//' D_0' - if (prm%atomicVolume <= 0.0_pReal) extmsg = trim(extmsg)//' V_at' ! ToDo: in disloTungsten, the atomic volume is given as a factor + if (prm%eta <= 0.0_pReal) extmsg = trim(extmsg)//' eta' + if (prm%Q_cl <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl' + if (prm%nu_a <= 0.0_pReal) extmsg = trim(extmsg)//' nu_a' + if (prm%w <= 0.0_pReal) extmsg = trim(extmsg)//' w' + if (prm%D_0 < 0.0_pReal) extmsg = trim(extmsg)//' D_0' + if (prm%V_at <= 0.0_pReal) extmsg = trim(extmsg)//' V_at' ! ToDo: in disloTungsten, the atomic volume is given as a factor - if (prm%significantN < 0.0_pReal) extmsg = trim(extmsg)//' rho_num_significant' - if (prm%significantrho < 0.0_pReal) extmsg = trim(extmsg)//' rho_significant' - if (prm%atol_rho < 0.0_pReal) extmsg = trim(extmsg)//' atol_rho' - if (prm%CFLfactor < 0.0_pReal) extmsg = trim(extmsg)//' f_c' + if (prm%rho_min < 0.0_pReal) extmsg = trim(extmsg)//' rho_min' + if (prm%rho_significant < 0.0_pReal) extmsg = trim(extmsg)//' rho_significant' + if (prm%atol_rho < 0.0_pReal) extmsg = trim(extmsg)//' atol_rho' + if (prm%f_c < 0.0_pReal) extmsg = trim(extmsg)//' f_c' - if (prm%p <= 0.0_pReal .or. prm%p > 1.0_pReal) extmsg = trim(extmsg)//' p_sl' - if (prm%q < 1.0_pReal .or. prm%q > 2.0_pReal) extmsg = trim(extmsg)//' q_sl' + if (prm%p <= 0.0_pReal .or. prm%p > 1.0_pReal) extmsg = trim(extmsg)//' p_sl' + if (prm%q < 1.0_pReal .or. prm%q > 2.0_pReal) extmsg = trim(extmsg)//' q_sl' - if (prm%linetensionEffect < 0.0_pReal .or. prm%linetensionEffect > 1.0_pReal) & + if (prm%f_F < 0.0_pReal .or. prm%f_F > 1.0_pReal) & extmsg = trim(extmsg)//' f_F' - if (prm%edgeJogFactor < 0.0_pReal .or. prm%edgeJogFactor > 1.0_pReal) & + if (prm%f_ed < 0.0_pReal .or. prm%f_ed > 1.0_pReal) & extmsg = trim(extmsg)//' f_ed' - if (prm%solidSolutionEnergy <= 0.0_pReal) extmsg = trim(extmsg)//' Q_sol' - if (prm%solidSolutionSize <= 0.0_pReal) extmsg = trim(extmsg)//' f_sol' - if (prm%solidSolutionConcentration <= 0.0_pReal) extmsg = trim(extmsg)//' c_sol' + if (prm%Q_sol <= 0.0_pReal) extmsg = trim(extmsg)//' Q_sol' + if (prm%f_sol <= 0.0_pReal) extmsg = trim(extmsg)//' f_sol' + if (prm%c_sol <= 0.0_pReal) extmsg = trim(extmsg)//' c_sol' - if (prm%grainboundaryTransmissivity > 1.0_pReal) extmsg = trim(extmsg)//' chi_GB' - if (prm%surfaceTransmissivity < 0.0_pReal .or. prm%surfaceTransmissivity > 1.0_pReal) & - extmsg = trim(extmsg)//' chi_surface' + if (prm%chi_GB > 1.0_pReal) extmsg = trim(extmsg)//' chi_GB' + if (prm%chi_surface < 0.0_pReal .or. prm%chi_surface > 1.0_pReal) & + extmsg = trim(extmsg)//' chi_surface' - if (prm%fEdgeMultiplication < 0.0_pReal .or. prm%fEdgeMultiplication > 1.0_pReal) & - extmsg = trim(extmsg)//' f_ed_mult' + if (prm%f_ed_mult < 0.0_pReal .or. prm%f_ed_mult > 1.0_pReal) & + extmsg = trim(extmsg)//' f_ed_mult' endif slipActive @@ -412,7 +412,7 @@ module function plastic_nonlocal_init() result(myPlasticity) call IO_error(212,ext_msg='IPneighborhood does not exist') - plasticState(p)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention + plasticState(p)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention st0%rho => plasticState(p)%state0 (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:) stt%rho => plasticState(p)%state (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:) @@ -620,16 +620,16 @@ module subroutine plastic_nonlocal_dependentState(F, Fp, instance, of, ip, el) ! coefficients are corrected for the line tension effect ! (see Kubin,Devincre,Hoc; 2008; Modeling dislocation storage rates and mean free paths in face-centered cubic crystals) if (any(lattice_structure(material_phaseAt(1,el)) == [LATTICE_bcc_ID,LATTICE_fcc_ID])) then - myInteractionMatrix = prm%interactionSlipSlip & - * spread(( 1.0_pReal - prm%linetensionEffect & - + prm%linetensionEffect & - * log(0.35_pReal * prm%burgers * sqrt(max(stt%rho_forest(:,of),prm%significantRho))) & - / log(0.35_pReal * prm%burgers * 1e6_pReal))** 2.0_pReal,2,prm%sum_N_sl) + myInteractionMatrix = prm%h_sl_sl & + * spread(( 1.0_pReal - prm%f_F & + + prm%f_F & + * log(0.35_pReal * prm%b_sl * sqrt(max(stt%rho_forest(:,of),prm%rho_significant))) & + / log(0.35_pReal * prm%b_sl * 1e6_pReal))** 2.0_pReal,2,prm%sum_N_sl) else - myInteractionMatrix = prm%interactionSlipSlip + myInteractionMatrix = prm%h_sl_sl endif - dst%tau_pass(:,of) = prm%mu * prm%burgers & + dst%tau_pass(:,of) = prm%mu * prm%b_sl & * sqrt(matmul(myInteractionMatrix,sum(abs(rho),2))) !*** calculate the dislocation stress of the neighboring excess dislocation densities @@ -731,7 +731,7 @@ module subroutine plastic_nonlocal_dependentState(F, Fp, instance, of, ip, el) where(rhoTotal > 0.0_pReal) rhoExcessGradient_over_rho = rhoExcessGradient / rhoTotal ! ... gives the local stress correction when multiplied with a factor - dst%tau_back(s,of) = - prm%mu * prm%burgers(s) / (2.0_pReal * PI) & + dst%tau_back(s,of) = - prm%mu * prm%b_sl(s) / (2.0_pReal * PI) & * ( rhoExcessGradient_over_rho(1) / (1.0_pReal - prm%nu) & + rhoExcessGradient_over_rho(2)) enddo @@ -841,7 +841,7 @@ module subroutine plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp, & forall (s = 1:ns, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) & rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t)) - gdotTotal = sum(rhoSgl(:,1:4) * v, 2) * prm%burgers + gdotTotal = sum(rhoSgl(:,1:4) * v, 2) * prm%b_sl Lp = 0.0_pReal dLp_dMp = 0.0_pReal @@ -850,10 +850,10 @@ module subroutine plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp, & forall (i=1:3,j=1:3,k=1:3,l=1:3) & dLp_dMp(i,j,k,l) = dLp_dMp(i,j,k,l) & + prm%Schmid(i,j,s) * prm%Schmid(k,l,s) & - * sum(rhoSgl(s,1:4) * dv_dtau(s,1:4)) * prm%burgers(s) & + * sum(rhoSgl(s,1:4) * dv_dtau(s,1:4)) * prm%b_sl(s) & + prm%Schmid(i,j,s) & * ( prm%nonSchmid_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) & - - prm%nonSchmid_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4)) * prm%burgers(s) + - prm%nonSchmid_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4)) * prm%b_sl(s) enddo end associate @@ -929,8 +929,8 @@ module subroutine plastic_nonlocal_deltaState(Mp,instance,of,ip,el) if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal enddo - dUpper(:,1) = prm%mu * prm%burgers/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau)) - dUpper(:,2) = prm%mu * prm%burgers/(4.0_pReal * PI * abs(tau)) + dUpper(:,1) = prm%mu * prm%b_sl/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau)) + dUpper(:,2) = prm%mu * prm%b_sl/(4.0_pReal * PI * abs(tau)) where(dNeq0(sqrt(sum(abs(rho(:,edg)),2)))) & dUpper(:,1) = min(1.0_pReal/sqrt(sum(abs(rho(:,edg)),2)),dUpper(:,1)) @@ -1041,7 +1041,7 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & my_rhoSgl0 = rho0(:,sgl) forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),of) - gdot = rhoSgl(:,1:4) * v * spread(prm%burgers,2,4) + gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) #ifdef DEBUG if (debugConstitutive%basic & @@ -1060,8 +1060,8 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & enddo dLower = prm%minDipoleHeight - dUpper(:,1) = prm%mu * prm%burgers/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau)) - dUpper(:,2) = prm%mu * prm%burgers/(4.0_pReal * PI * abs(tau)) + dUpper(:,1) = prm%mu * prm%b_sl/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau)) + dUpper(:,2) = prm%mu * prm%b_sl/(4.0_pReal * PI * abs(tau)) where(dNeq0(sqrt(sum(abs(rho(:,edg)),2)))) & dUpper(:,1) = min(1.0_pReal/sqrt(sum(abs(rho(:,edg)),2)),dUpper(:,1)) @@ -1075,18 +1075,18 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & rhoDotMultiplication = 0.0_pReal isBCC: if (lattice_structure(ph) == LATTICE_bcc_ID) then forall (s = 1:ns, sum(abs(v(s,1:4))) > 0.0_pReal) - rhoDotMultiplication(s,1:2) = sum(abs(gdot(s,3:4))) / prm%burgers(s) & ! assuming double-cross-slip of screws to be decisive for multiplication - * sqrt(stt%rho_forest(s,of)) / prm%lambda0(s) ! & ! mean free path + rhoDotMultiplication(s,1:2) = sum(abs(gdot(s,3:4))) / prm%b_sl(s) & ! assuming double-cross-slip of screws to be decisive for multiplication + * sqrt(stt%rho_forest(s,of)) / prm%i_sl(s) ! & ! mean free path ! * 2.0_pReal * sum(abs(v(s,3:4))) / sum(abs(v(s,1:4))) ! ratio of screw to overall velocity determines edge generation - rhoDotMultiplication(s,3:4) = sum(abs(gdot(s,3:4))) /prm%burgers(s) & ! assuming double-cross-slip of screws to be decisive for multiplication - * sqrt(stt%rho_forest(s,of)) / prm%lambda0(s) ! & ! mean free path + rhoDotMultiplication(s,3:4) = sum(abs(gdot(s,3:4))) /prm%b_sl(s) & ! assuming double-cross-slip of screws to be decisive for multiplication + * sqrt(stt%rho_forest(s,of)) / prm%i_sl(s) ! & ! mean free path ! * 2.0_pReal * sum(abs(v(s,1:2))) / sum(abs(v(s,1:4))) ! ratio of edge to overall velocity determines screw generation endforall else isBCC rhoDotMultiplication(:,1:4) = spread( & - (sum(abs(gdot(:,1:2)),2) * prm%fEdgeMultiplication + sum(abs(gdot(:,3:4)),2)) & - * sqrt(stt%rho_forest(:,of)) / prm%lambda0 / prm%burgers, 2, 4) + (sum(abs(gdot(:,1:2)),2) * prm%f_ed_mult + sum(abs(gdot(:,3:4)),2)) & + * sqrt(stt%rho_forest(:,of)) / prm%i_sl / prm%b_sl, 2, 4) endif isBCC forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,instance),of) @@ -1097,20 +1097,20 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & !*** formation by glide do c = 1,2 - rhoDotSingle2DipoleGlide(:,2*c-1) = -2.0_pReal * dUpper(:,c) / prm%burgers & + rhoDotSingle2DipoleGlide(:,2*c-1) = -2.0_pReal * dUpper(:,c) / prm%b_sl & * ( rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) & ! negative mobile --> positive mobile + rhoSgl(:,2*c) * abs(gdot(:,2*c-1)) & ! positive mobile --> negative mobile + abs(rhoSgl(:,2*c+4)) * abs(gdot(:,2*c-1))) ! positive mobile --> negative immobile - rhoDotSingle2DipoleGlide(:,2*c) = -2.0_pReal * dUpper(:,c) / prm%burgers & + rhoDotSingle2DipoleGlide(:,2*c) = -2.0_pReal * dUpper(:,c) / prm%b_sl & * ( rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) & ! negative mobile --> positive mobile + rhoSgl(:,2*c) * abs(gdot(:,2*c-1)) & ! positive mobile --> negative mobile + abs(rhoSgl(:,2*c+3)) * abs(gdot(:,2*c))) ! negative mobile --> positive immobile - rhoDotSingle2DipoleGlide(:,2*c+3) = -2.0_pReal * dUpper(:,c) / prm%burgers & + rhoDotSingle2DipoleGlide(:,2*c+3) = -2.0_pReal * dUpper(:,c) / prm%b_sl & * rhoSgl(:,2*c+3) * abs(gdot(:,2*c)) ! negative mobile --> positive immobile - rhoDotSingle2DipoleGlide(:,2*c+4) = -2.0_pReal * dUpper(:,c) / prm%burgers & + rhoDotSingle2DipoleGlide(:,2*c+4) = -2.0_pReal * dUpper(:,c) / prm%b_sl & * rhoSgl(:,2*c+4) * abs(gdot(:,2*c-1)) ! positive mobile --> negative immobile rhoDotSingle2DipoleGlide(:,c+8) = abs(rhoDotSingle2DipoleGlide(:,2*c+3)) & @@ -1123,27 +1123,27 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & !*** athermal annihilation rhoDotAthermalAnnihilation = 0.0_pReal forall (c=1:2) & - rhoDotAthermalAnnihilation(:,c+8) = -2.0_pReal * dLower(:,c) / prm%burgers & - * ( 2.0_pReal * (rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) + rhoSgl(:,2*c) * abs(gdot(:,2*c-1))) & ! was single hitting single + rhoDotAthermalAnnihilation(:,c+8) = -2.0_pReal * dLower(:,c) / prm%b_sl & + * ( 2.0_pReal * (rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) + rhoSgl(:,2*c) * abs(gdot(:,2*c-1))) & ! was single hitting single + 2.0_pReal * (abs(rhoSgl(:,2*c+3)) * abs(gdot(:,2*c)) + abs(rhoSgl(:,2*c+4)) * abs(gdot(:,2*c-1))) & ! was single hitting immobile single or was immobile single hit by single - + rhoDip(:,c) * (abs(gdot(:,2*c-1)) + abs(gdot(:,2*c)))) ! single knocks dipole constituent + + rhoDip(:,c) * (abs(gdot(:,2*c-1)) + abs(gdot(:,2*c)))) ! single knocks dipole constituent ! annihilated screw dipoles leave edge jogs behind on the colinear system if (lattice_structure(ph) == LATTICE_fcc_ID) & forall (s = 1:ns, prm%colinearSystem(s) > 0) & rhoDotAthermalAnnihilation(prm%colinearSystem(s),1:2) = - rhoDotAthermalAnnihilation(s,10) & - * 0.25_pReal * sqrt(stt%rho_forest(s,of)) * (dUpper(s,2) + dLower(s,2)) * prm%edgeJogFactor + * 0.25_pReal * sqrt(stt%rho_forest(s,of)) * (dUpper(s,2) + dLower(s,2)) * prm%f_ed !*** thermally activated annihilation of edge dipoles by climb rhoDotThermalAnnihilation = 0.0_pReal - selfDiffusion = prm%Dsd0 * exp(-prm%selfDiffusionEnergy / (kB * Temperature)) - vClimb = prm%atomicVolume * selfDiffusion * prm%mu & + selfDiffusion = prm%D_0 * exp(-prm%Q_cl / (kB * Temperature)) + vClimb = prm%V_at * selfDiffusion * prm%mu & / ( kB * Temperature * PI * (1.0_pReal-prm%nu) * (dUpper(:,1) + dLower(:,1))) forall (s = 1:ns, dUpper(s,1) > dLower(s,1)) & rhoDotThermalAnnihilation(s,9) = max(- 4.0_pReal * rhoDip(s,1) * vClimb(s) / (dUpper(s,1) - dLower(s,1)), & - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & - - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have + - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have rhoDot = rhoDotFlux(F,Fp,timestep, instance,of,ip,el) & + rhoDotMultiplication & @@ -1252,7 +1252,7 @@ function rhoDotFlux(F,Fp,timestep, instance,of,ip,el) my_rhoSgl0 = rho0(:,sgl) forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),of) !ToDo: MD: I think we should use state0 here - gdot = rhoSgl(:,1:4) * v * spread(prm%burgers,2,4) + gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,instance),of) @@ -1264,14 +1264,14 @@ function rhoDotFlux(F,Fp,timestep, instance,of,ip,el) !*** check CFL (Courant-Friedrichs-Lewy) condition for flux if (any( abs(gdot) > 0.0_pReal & ! any active slip system ... - .and. prm%CFLfactor * abs(v0) * timestep & + .and. prm%f_c * abs(v0) * timestep & > IPvolume(ip,el) / maxval(IParea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here) #ifdef DEBUG if (debugConstitutive%extensive) then print'(a,i5,a,i2)', '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip print'(a,e10.3,a,e10.3)', '<< CONST >> velocity is at ', & maxval(abs(v0), abs(gdot) > 0.0_pReal & - .and. prm%CFLfactor * abs(v0) * timestep & + .and. prm%f_c * abs(v0) * timestep & > IPvolume(ip,el) / maxval(IParea(:,ip,el))), & ' at a timestep of ',timestep print*, '<< CONST >> enforcing cutback !!!' @@ -1334,8 +1334,8 @@ function rhoDotFlux(F,Fp,timestep, instance,of,ip,el) neighbor_rhoSgl0(s,t) = max(plasticState(np)%state0(iRhoU(s,t,neighbor_instance),no),0.0_pReal) endforall - where (neighbor_rhoSgl0 * IPvolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%significantN & - .or. neighbor_rhoSgl0 < prm%significantRho) & + where (neighbor_rhoSgl0 * IPvolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%rho_min & + .or. neighbor_rhoSgl0 < prm%rho_significant) & neighbor_rhoSgl0 = 0.0_pReal normal_neighbor2me_defConf = math_det33(Favg) * matmul(math_inv33(transpose(Favg)), & IPareaNormal(1:3,neighbor_n,neighbor_ip,neighbor_el)) ! normal of the interface in (average) deformed configuration (pointing neighbor => me) @@ -1460,7 +1460,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e) if (neighbor_e <= 0 .or. neighbor_i <= 0) then !* FREE SURFACE !* Set surface transmissivity to the value specified in the material.config - forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%surfaceTransmissivity) + forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_surface) elseif (neighbor_phase /= ph) then !* PHASE BOUNDARY !* If we encounter a different nonlocal phase at the neighbor, @@ -1469,13 +1469,13 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e) !* we do not consider this to be a phase boundary, so completely compatible. if (.not. phase_localPlasticity(neighbor_phase) .and. .not. phase_localPlasticity(ph)) & forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = 0.0_pReal - elseif (prm%grainboundaryTransmissivity >= 0.0_pReal) then + elseif (prm%chi_GB >= 0.0_pReal) then !* GRAIN BOUNDARY ! !* fixed transmissivity for adjacent ips with different texture (only if explicitly given in material.config) if (any(dNeq(material_orientation0(1,i,e)%asQuaternion(), & material_orientation0(1,neighbor_i,neighbor_e)%asQuaternion())) .and. & (.not. phase_localPlasticity(neighbor_phase))) & - forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%grainboundaryTransmissivity) + forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB) else !* GRAIN BOUNDARY ? !* Compatibility defined by relative orientation of slip systems: @@ -1631,7 +1631,7 @@ subroutine stateInit(ini,phase,NipcMyPhase,instance) associate(stt => state(instance)) - if (ini%rhoSglRandom > 0.0_pReal) then ! randomly distribute dislocation segments on random slip system and of random type in the volume + if (ini%random_rho_u > 0.0_pReal) then ! randomly distribute dislocation segments on random slip system and of random type in the volume do e = 1,discretization_nElem do i = 1,discretization_nIP if (material_phaseAt(1,e) == phase) volume(material_phasememberAt(1,i,e)) = IPvolume(i,e) @@ -1639,11 +1639,11 @@ subroutine stateInit(ini,phase,NipcMyPhase,instance) enddo totalVolume = sum(volume) minimumIPVolume = minval(volume) - densityBinning = ini%rhoSglRandomBinning / minimumIpVolume ** (2.0_pReal / 3.0_pReal) + densityBinning = ini%random_rho_u_binning / minimumIpVolume ** (2.0_pReal / 3.0_pReal) ! fill random material points with dislocation segments until the desired overall density is reached meanDensity = 0.0_pReal - do while(meanDensity < ini%rhoSglRandom) + do while(meanDensity < ini%random_rho_u) call random_number(rnd) phasemember = nint(rnd(1)*real(NipcMyPhase,pReal) + 0.5_pReal) s = nint(rnd(2)*real(sum(ini%N_sl),pReal)*4.0_pReal + 0.5_pReal) @@ -1656,15 +1656,15 @@ subroutine stateInit(ini,phase,NipcMyPhase,instance) from = 1 + sum(ini%N_sl(1:f-1)) upto = sum(ini%N_sl(1:f)) do s = from,upto - noise = [math_sampleGaussVar(0.0_pReal, ini%rhoSglScatter), & - math_sampleGaussVar(0.0_pReal, ini%rhoSglScatter)] - stt%rho_sgl_mob_edg_pos(s,e) = ini%rhoSglEdgePos0(f) + noise(1) - stt%rho_sgl_mob_edg_neg(s,e) = ini%rhoSglEdgeNeg0(f) + noise(1) - stt%rho_sgl_mob_scr_pos(s,e) = ini%rhoSglScrewPos0(f) + noise(2) - stt%rho_sgl_mob_scr_neg(s,e) = ini%rhoSglScrewNeg0(f) + noise(2) + noise = [math_sampleGaussVar(0.0_pReal, ini%sigma_rho_u), & + math_sampleGaussVar(0.0_pReal, ini%sigma_rho_u)] + stt%rho_sgl_mob_edg_pos(s,e) = ini%rho_u_ed_pos_0(f) + noise(1) + stt%rho_sgl_mob_edg_neg(s,e) = ini%rho_u_ed_neg_0(f) + noise(1) + stt%rho_sgl_mob_scr_pos(s,e) = ini%rho_u_sc_pos_0(f) + noise(2) + stt%rho_sgl_mob_scr_neg(s,e) = ini%rho_u_sc_neg_0(f) + noise(2) enddo - stt%rho_dip_edg(from:upto,e) = ini%rhoDipEdge0(f) - stt%rho_dip_scr(from:upto,e) = ini%rhoDipScrew0(f) + stt%rho_dip_edg(from:upto,e) = ini%rho_d_ed_0(f) + stt%rho_dip_scr(from:upto,e) = ini%rho_d_sc_0(f) enddo enddo endif @@ -1732,14 +1732,14 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Tem !* Effective stress includes non Schmid constributions !* The derivative only gives absolute values; the correct sign is taken care of in the formula for the derivative of the velocity tauEff = max(0.0_pReal, abs(tauNS(s)) - tauThreshold(s)) ! ensure that the effective stress is positive - meanfreepath_P = prm%burgers(s) - jumpWidth_P = prm%burgers(s) - activationLength_P = prm%doublekinkwidth *prm%burgers(s) - activationVolume_P = activationLength_P * jumpWidth_P * prm%burgers(s) + meanfreepath_P = prm%b_sl(s) + jumpWidth_P = prm%b_sl(s) + activationLength_P = prm%w *prm%b_sl(s) + activationVolume_P = activationLength_P * jumpWidth_P * prm%b_sl(s) criticalStress_P = prm%peierlsStress(s,c) activationEnergy_P = criticalStress_P * activationVolume_P tauRel_P = min(1.0_pReal, tauEff / criticalStress_P) ! ensure that the activation probability cannot become greater than one - tPeierls = 1.0_pReal / prm%fattack & + tPeierls = 1.0_pReal / prm%nu_a & * exp(activationEnergy_P / (kB * Temperature) & * (1.0_pReal - tauRel_P**prm%p)**prm%q) if (tauEff < criticalStress_P) then @@ -1752,14 +1752,14 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Tem !* Contribution from solid solution strengthening !* The derivative only gives absolute values; the correct sign is taken care of in the formula for the derivative of the velocity tauEff = abs(tau(s)) - tauThreshold(s) - meanfreepath_S = prm%burgers(s) / sqrt(prm%solidSolutionConcentration) - jumpWidth_S = prm%solidSolutionSize * prm%burgers(s) - activationLength_S = prm%burgers(s) / sqrt(prm%solidSolutionConcentration) - activationVolume_S = activationLength_S * jumpWidth_S * prm%burgers(s) - activationEnergy_S = prm%solidSolutionEnergy + meanfreepath_S = prm%b_sl(s) / sqrt(prm%c_sol) + jumpWidth_S = prm%f_sol * prm%b_sl(s) + activationLength_S = prm%b_sl(s) / sqrt(prm%c_sol) + activationVolume_S = activationLength_S * jumpWidth_S * prm%b_sl(s) + activationEnergy_S = prm%Q_sol criticalStress_S = activationEnergy_S / activationVolume_S tauRel_S = min(1.0_pReal, tauEff / criticalStress_S) ! ensure that the activation probability cannot become greater than one - tSolidSolution = 1.0_pReal / prm%fattack & + tSolidSolution = 1.0_pReal / prm%nu_a & * exp(activationEnergy_S / (kB * Temperature)* (1.0_pReal - tauRel_S**prm%p)**prm%q) if (tauEff < criticalStress_S) then dtSolidSolution_dtau = tSolidSolution * prm%p * prm%q * activationVolume_S / (kB * Temperature) & @@ -1770,7 +1770,7 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Tem !* viscous glide velocity tauEff = abs(tau(s)) - tauThreshold(s) - mobility = prm%burgers(s) / prm%viscosity + mobility = prm%b_sl(s) / prm%eta vViscous = mobility * tauEff !* Mean velocity results from waiting time at peierls barriers and solid solution obstacles with respective meanfreepath of @@ -1805,7 +1805,7 @@ pure function getRho(instance,of,ip,el) getRho(:,mob) = max(getRho(:,mob),0.0_pReal) getRho(:,dip) = max(getRho(:,dip),0.0_pReal) - where(abs(getRho) < max(prm%significantN/IPvolume(ip,el)**(2.0_pReal/3.0_pReal),prm%significantRho)) & + where(abs(getRho) < max(prm%rho_min/IPvolume(ip,el)**(2.0_pReal/3.0_pReal),prm%rho_significant)) & getRho = 0.0_pReal end associate @@ -1830,7 +1830,7 @@ pure function getRho0(instance,of,ip,el) getRho0(:,mob) = max(getRho0(:,mob),0.0_pReal) getRho0(:,dip) = max(getRho0(:,dip),0.0_pReal) - where(abs(getRho0) < max(prm%significantN/IPvolume(ip,el)**(2.0_pReal/3.0_pReal),prm%significantRho)) & + where(abs(getRho0) < max(prm%rho_min/IPvolume(ip,el)**(2.0_pReal/3.0_pReal),prm%rho_significant)) & getRho0 = 0.0_pReal end associate diff --git a/src/constitutive_plastic_phenopowerlaw.f90 b/src/constitutive_plastic_phenopowerlaw.f90 index 8cf63a54d..fd1b5fbbb 100644 --- a/src/constitutive_plastic_phenopowerlaw.f90 +++ b/src/constitutive_plastic_phenopowerlaw.f90 @@ -8,28 +8,28 @@ submodule(constitutive:constitutive_plastic) plastic_phenopowerlaw type :: tParameters real(pReal) :: & - gdot0_slip = 1.0_pReal, & !< reference shear strain rate for slip - gdot0_twin = 1.0_pReal, & !< reference shear strain rate for twin - n_slip = 1.0_pReal, & !< stress exponent for slip - n_twin = 1.0_pReal, & !< stress exponent for twin - spr = 1.0_pReal, & !< push-up factor for slip saturation due to twinning - c_1 = 1.0_pReal, & - c_2 = 1.0_pReal, & - c_3 = 1.0_pReal, & - c_4 = 1.0_pReal, & - h0_SlipSlip = 1.0_pReal, & !< reference hardening slip - slip - h0_TwinSlip = 1.0_pReal, & !< reference hardening twin - slip - h0_TwinTwin = 1.0_pReal, & !< reference hardening twin - twin - a_slip = 1.0_pReal + dot_gamma_0_sl = 1.0_pReal, & !< reference shear strain rate for slip + dot_gamma_0_tw = 1.0_pReal, & !< reference shear strain rate for twin + n_sl = 1.0_pReal, & !< stress exponent for slip + n_tw = 1.0_pReal, & !< stress exponent for twin + f_sl_sat_tw = 1.0_pReal, & !< push-up factor for slip saturation due to twinning + c_1 = 1.0_pReal, & + c_2 = 1.0_pReal, & + c_3 = 1.0_pReal, & + c_4 = 1.0_pReal, & + h_0_sl_sl = 1.0_pReal, & !< reference hardening slip - slip + h_0_tw_sl = 1.0_pReal, & !< reference hardening twin - slip + h_0_tw_tw = 1.0_pReal, & !< reference hardening twin - twin + a_sl = 1.0_pReal real(pReal), allocatable, dimension(:) :: & - xi_slip_sat, & !< maximum critical shear stress for slip - H_int, & !< per family hardening activity (optional) - gamma_twin_char !< characteristic shear for twins + xi_inf_sl, & !< maximum critical shear stress for slip + h_int, & !< per family hardening activity (optional) + gamma_tw_char !< characteristic shear for twins real(pReal), allocatable, dimension(:,:) :: & - interaction_SlipSlip, & !< slip resistance from slip activity - interaction_SlipTwin, & !< slip resistance from twin activity - interaction_TwinSlip, & !< twin resistance from slip activity - interaction_TwinTwin !< twin resistance from twin activity + h_sl_sl, & !< slip resistance from slip activity + h_sl_tw, & !< slip resistance from twin activity + h_tw_sl, & !< twin resistance from slip activity + h_tw_tw !< twin resistance from twin activity real(pReal), allocatable, dimension(:,:,:) :: & P_sl, & P_tw, & @@ -78,8 +78,8 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) integer, dimension(:), allocatable :: & N_sl, N_tw real(pReal), dimension(:), allocatable :: & - xi_slip_0, & !< initial critical shear stress for slip - xi_twin_0, & !< initial critical shear stress for twin + xi_0_sl, & !< initial critical shear stress for slip + xi_0_tw, & !< initial critical shear stress for twin a !< non-Schmid coefficients character(len=pStringLen) :: & extmsg = '' @@ -120,7 +120,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) phase%get_asFloat('c/a',defaultVal=0.0_pReal)) if(phase%get_asString('lattice') == 'bcc') then - a = pl%get_asFloats('nonSchmid_coefficients',defaultVal=emptyRealArray) + a = pl%get_asFloats('a_nonSchmid',defaultVal=emptyRealArray) if(size(a) > 0) prm%nonSchmidActive = .true. prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,a,+1) prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,a,-1) @@ -128,36 +128,36 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) prm%nonSchmid_pos = prm%P_sl prm%nonSchmid_neg = prm%P_sl endif - prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(N_sl, & - pl%get_asFloats('h_sl_sl'), & - phase%get_asString('lattice')) + prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl, & + pl%get_asFloats('h_sl_sl'), & + phase%get_asString('lattice')) - xi_slip_0 = pl%get_asFloats('xi_0_sl', requiredSize=size(N_sl)) - prm%xi_slip_sat = pl%get_asFloats('xi_inf_sl', requiredSize=size(N_sl)) - prm%H_int = pl%get_asFloats('h_int', requiredSize=size(N_sl), & - defaultVal=[(0.0_pReal,i=1,size(N_sl))]) + xi_0_sl = pl%get_asFloats('xi_0_sl', requiredSize=size(N_sl)) + prm%xi_inf_sl = pl%get_asFloats('xi_inf_sl', requiredSize=size(N_sl)) + prm%h_int = pl%get_asFloats('h_int', requiredSize=size(N_sl), & + defaultVal=[(0.0_pReal,i=1,size(N_sl))]) - prm%gdot0_slip = pl%get_asFloat('dot_gamma_0_sl') - prm%n_slip = pl%get_asFloat('n_sl') - prm%a_slip = pl%get_asFloat('a_sl') - prm%h0_SlipSlip = pl%get_asFloat('h_0_sl_sl') + prm%dot_gamma_0_sl = pl%get_asFloat('dot_gamma_0_sl') + prm%n_sl = pl%get_asFloat('n_sl') + prm%a_sl = pl%get_asFloat('a_sl') + prm%h_0_sl_sl = pl%get_asFloat('h_0_sl_sl') ! expand: family => system - xi_slip_0 = math_expand(xi_slip_0, N_sl) - prm%xi_slip_sat = math_expand(prm%xi_slip_sat,N_sl) - prm%H_int = math_expand(prm%H_int, N_sl) + xi_0_sl = math_expand(xi_0_sl, N_sl) + prm%xi_inf_sl = math_expand(prm%xi_inf_sl,N_sl) + prm%h_int = math_expand(prm%h_int, N_sl) ! sanity checks - if ( prm%gdot0_slip <= 0.0_pReal) extmsg = trim(extmsg)//' dot_gamma_0_sl' - if ( prm%a_slip <= 0.0_pReal) extmsg = trim(extmsg)//' a_sl' - if ( prm%n_slip <= 0.0_pReal) extmsg = trim(extmsg)//' n_sl' - if (any(xi_slip_0 <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_0_sl' - if (any(prm%xi_slip_sat <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_inf_sl' + if ( prm%dot_gamma_0_sl <= 0.0_pReal) extmsg = trim(extmsg)//' dot_gamma_0_sl' + if ( prm%a_sl <= 0.0_pReal) extmsg = trim(extmsg)//' a_sl' + if ( prm%n_sl <= 0.0_pReal) extmsg = trim(extmsg)//' n_sl' + if (any(xi_0_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_0_sl' + if (any(prm%xi_inf_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_inf_sl' else slipActive - xi_slip_0 = emptyRealArray - allocate(prm%xi_slip_sat,prm%H_int,source=emptyRealArray) - allocate(prm%interaction_SlipSlip(0,0)) + xi_0_sl = emptyRealArray + allocate(prm%xi_inf_sl,prm%h_int,source=emptyRealArray) + allocate(prm%h_sl_sl(0,0)) endif slipActive !-------------------------------------------------------------------------------------------------- @@ -165,52 +165,52 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) N_tw = pl%get_asInts('N_tw', defaultVal=emptyIntArray) prm%sum_N_tw = sum(abs(N_tw)) twinActive: if (prm%sum_N_tw > 0) then - prm%P_tw = lattice_SchmidMatrix_twin(N_tw,phase%get_asString('lattice'),& - phase%get_asFloat('c/a',defaultVal=0.0_pReal)) - prm%interaction_TwinTwin = lattice_interaction_TwinByTwin(N_tw,& - pl%get_asFloats('h_tw_tw'), & - phase%get_asString('lattice')) - prm%gamma_twin_char = lattice_characteristicShear_twin(N_tw,phase%get_asString('lattice'),& - phase%get_asFloat('c/a',defaultVal=0.0_pReal)) + prm%P_tw = lattice_SchmidMatrix_twin(N_tw,phase%get_asString('lattice'),& + phase%get_asFloat('c/a',defaultVal=0.0_pReal)) + prm%h_tw_tw = lattice_interaction_TwinByTwin(N_tw,& + pl%get_asFloats('h_tw_tw'), & + phase%get_asString('lattice')) + prm%gamma_tw_char = lattice_characteristicShear_twin(N_tw,phase%get_asString('lattice'),& + phase%get_asFloat('c/a',defaultVal=0.0_pReal)) - xi_twin_0 = pl%get_asFloats('xi_0_tw',requiredSize=size(N_tw)) + xi_0_tw = pl%get_asFloats('xi_0_tw',requiredSize=size(N_tw)) - prm%c_1 = pl%get_asFloat('c_1',defaultVal=0.0_pReal) - prm%c_2 = pl%get_asFloat('c_2',defaultVal=1.0_pReal) - prm%c_3 = pl%get_asFloat('c_3',defaultVal=0.0_pReal) - prm%c_4 = pl%get_asFloat('c_4',defaultVal=0.0_pReal) - prm%gdot0_twin = pl%get_asFloat('dot_gamma_0_tw') - prm%n_twin = pl%get_asFloat('n_tw') - prm%spr = pl%get_asFloat('f_sl_sat_tw') - prm%h0_TwinTwin = pl%get_asFloat('h_0_tw_tw') + prm%c_1 = pl%get_asFloat('c_1',defaultVal=0.0_pReal) + prm%c_2 = pl%get_asFloat('c_2',defaultVal=1.0_pReal) + prm%c_3 = pl%get_asFloat('c_3',defaultVal=0.0_pReal) + prm%c_4 = pl%get_asFloat('c_4',defaultVal=0.0_pReal) + prm%dot_gamma_0_tw = pl%get_asFloat('dot_gamma_0_tw') + prm%n_tw = pl%get_asFloat('n_tw') + prm%f_sl_sat_tw = pl%get_asFloat('f_sl_sat_tw') + prm%h_0_tw_tw = pl%get_asFloat('h_0_tw_tw') ! expand: family => system - xi_twin_0 = math_expand(xi_twin_0,N_tw) + xi_0_tw = math_expand(xi_0_tw,N_tw) ! sanity checks - if (prm%gdot0_twin <= 0.0_pReal) extmsg = trim(extmsg)//' dot_gamma_0_tw' - if (prm%n_twin <= 0.0_pReal) extmsg = trim(extmsg)//' n_tw' + if (prm%dot_gamma_0_tw <= 0.0_pReal) extmsg = trim(extmsg)//' dot_gamma_0_tw' + if (prm%n_tw <= 0.0_pReal) extmsg = trim(extmsg)//' n_tw' else twinActive - xi_twin_0 = emptyRealArray - allocate(prm%gamma_twin_char,source=emptyRealArray) - allocate(prm%interaction_TwinTwin(0,0)) + xi_0_tw = emptyRealArray + allocate(prm%gamma_tw_char,source=emptyRealArray) + allocate(prm%h_tw_tw(0,0)) endif twinActive !-------------------------------------------------------------------------------------------------- ! slip-twin related parameters slipAndTwinActive: if (prm%sum_N_sl > 0 .and. prm%sum_N_tw > 0) then - prm%h0_TwinSlip = pl%get_asFloat('h_0_tw_sl') - prm%interaction_SlipTwin = lattice_interaction_SlipByTwin(N_sl,N_tw,& - pl%get_asFloats('h_sl_tw'), & - phase%get_asString('lattice')) - prm%interaction_TwinSlip = lattice_interaction_TwinBySlip(N_tw,N_sl,& - pl%get_asFloats('h_tw_sl'), & - phase%get_asString('lattice')) + prm%h_0_tw_sl = pl%get_asFloat('h_0_tw_sl') + prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,N_tw,& + pl%get_asFloats('h_sl_tw'), & + phase%get_asString('lattice')) + prm%h_tw_sl = lattice_interaction_TwinBySlip(N_tw,N_sl,& + pl%get_asFloats('h_tw_sl'), & + phase%get_asString('lattice')) else slipAndTwinActive - allocate(prm%interaction_SlipTwin(prm%sum_N_sl,prm%sum_N_tw)) ! at least one dimension is 0 - allocate(prm%interaction_TwinSlip(prm%sum_N_tw,prm%sum_N_sl)) ! at least one dimension is 0 - prm%h0_TwinSlip = 0.0_pReal + allocate(prm%h_sl_tw(prm%sum_N_sl,prm%sum_N_tw)) ! at least one dimension is 0 + allocate(prm%h_tw_sl(prm%sum_N_tw,prm%sum_N_sl)) ! at least one dimension is 0 + prm%h_0_tw_sl = 0.0_pReal endif slipAndTwinActive !-------------------------------------------------------------------------------------------------- @@ -237,7 +237,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) startIndex = 1 endIndex = prm%sum_N_sl stt%xi_slip => plasticState(p)%state (startIndex:endIndex,:) - stt%xi_slip = spread(xi_slip_0, 2, NipcMyPhase) + stt%xi_slip = spread(xi_0_sl, 2, NipcMyPhase) dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' @@ -245,7 +245,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tw stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:) - stt%xi_twin = spread(xi_twin_0, 2, NipcMyPhase) + stt%xi_twin = spread(xi_0_tw, 2, NipcMyPhase) dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' @@ -354,20 +354,20 @@ module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) sumGamma = sum(stt%gamma_slip(:,of)) - sumF = sum(stt%gamma_twin(:,of)/prm%gamma_twin_char) + sumF = sum(stt%gamma_twin(:,of)/prm%gamma_tw_char) !-------------------------------------------------------------------------------------------------- ! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices - c_SlipSlip = prm%h0_slipslip * (1.0_pReal + prm%c_1*sumF** prm%c_2) - c_TwinSlip = prm%h0_TwinSlip * sumGamma**prm%c_3 - c_TwinTwin = prm%h0_TwinTwin * sumF**prm%c_4 + c_SlipSlip = prm%h_0_sl_sl * (1.0_pReal + prm%c_1*sumF** prm%c_2) + c_TwinSlip = prm%h_0_tw_sl * sumGamma**prm%c_3 + c_TwinTwin = prm%h_0_tw_tw * sumF**prm%c_4 !-------------------------------------------------------------------------------------------------- ! calculate left and right vectors - left_SlipSlip = 1.0_pReal + prm%H_int - xi_slip_sat_offset = prm%spr*sqrt(sumF) - right_SlipSlip = abs(1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset)) **prm%a_slip & - * sign(1.0_pReal,1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset)) + left_SlipSlip = 1.0_pReal + prm%h_int + xi_slip_sat_offset = prm%f_sl_sat_tw*sqrt(sumF) + right_SlipSlip = abs(1.0_pReal-stt%xi_slip(:,of) / (prm%xi_inf_sl+xi_slip_sat_offset)) **prm%a_sl & + * sign(1.0_pReal,1.0_pReal-stt%xi_slip(:,of) / (prm%xi_inf_sl+xi_slip_sat_offset)) !-------------------------------------------------------------------------------------------------- ! shear rates @@ -378,11 +378,11 @@ module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) !-------------------------------------------------------------------------------------------------- ! hardening dot%xi_slip(:,of) = c_SlipSlip * left_SlipSlip * & - matmul(prm%interaction_SlipSlip,dot%gamma_slip(:,of)*right_SlipSlip) & - + matmul(prm%interaction_SlipTwin,dot%gamma_twin(:,of)) + matmul(prm%h_sl_sl,dot%gamma_slip(:,of)*right_SlipSlip) & + + matmul(prm%h_sl_tw,dot%gamma_twin(:,of)) - dot%xi_twin(:,of) = c_TwinSlip * matmul(prm%interaction_TwinSlip,dot%gamma_slip(:,of)) & - + c_TwinTwin * matmul(prm%interaction_TwinTwin,dot%gamma_twin(:,of)) + dot%xi_twin(:,of) = c_TwinSlip * matmul(prm%h_tw_sl,dot%gamma_slip(:,of)) & + + c_TwinTwin * matmul(prm%h_tw_tw,dot%gamma_twin(:,of)) end associate end subroutine plastic_phenopowerlaw_dotState @@ -460,29 +460,29 @@ pure subroutine kinetics_slip(Mp,instance,of, & enddo where(dNeq0(tau_slip_pos)) - gdot_slip_pos = prm%gdot0_slip * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active - * sign(abs(tau_slip_pos/stt%xi_slip(:,of))**prm%n_slip, tau_slip_pos) + gdot_slip_pos = prm%dot_gamma_0_sl * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active + * sign(abs(tau_slip_pos/stt%xi_slip(:,of))**prm%n_sl, tau_slip_pos) else where gdot_slip_pos = 0.0_pReal end where where(dNeq0(tau_slip_neg)) - gdot_slip_neg = prm%gdot0_slip * 0.5_pReal & ! only used if non-Schmid active, always 1/2 - * sign(abs(tau_slip_neg/stt%xi_slip(:,of))**prm%n_slip, tau_slip_neg) + gdot_slip_neg = prm%dot_gamma_0_sl * 0.5_pReal & ! only used if non-Schmid active, always 1/2 + * sign(abs(tau_slip_neg/stt%xi_slip(:,of))**prm%n_sl, tau_slip_neg) else where gdot_slip_neg = 0.0_pReal end where if (present(dgdot_dtau_slip_pos)) then where(dNeq0(gdot_slip_pos)) - dgdot_dtau_slip_pos = gdot_slip_pos*prm%n_slip/tau_slip_pos + dgdot_dtau_slip_pos = gdot_slip_pos*prm%n_sl/tau_slip_pos else where dgdot_dtau_slip_pos = 0.0_pReal end where endif if (present(dgdot_dtau_slip_neg)) then where(dNeq0(gdot_slip_neg)) - dgdot_dtau_slip_neg = gdot_slip_neg*prm%n_slip/tau_slip_neg + dgdot_dtau_slip_neg = gdot_slip_neg*prm%n_sl/tau_slip_neg else where dgdot_dtau_slip_neg = 0.0_pReal end where @@ -524,15 +524,15 @@ pure subroutine kinetics_twin(Mp,instance,of,& enddo where(tau_twin > 0.0_pReal) - gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)) & ! only twin in untwinned volume fraction - * prm%gdot0_twin*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_twin + gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,of)/prm%gamma_tw_char)) & ! only twin in untwinned volume fraction + * prm%dot_gamma_0_tw*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_tw else where gdot_twin = 0.0_pReal end where if (present(dgdot_dtau_twin)) then where(dNeq0(gdot_twin)) - dgdot_dtau_twin = gdot_twin*prm%n_twin/tau_twin + dgdot_dtau_twin = gdot_twin*prm%n_tw/tau_twin else where dgdot_dtau_twin = 0.0_pReal end where diff --git a/src/damage_nonlocal.f90 b/src/damage_nonlocal.f90 index a4ea54493..f31b10242 100644 --- a/src/damage_nonlocal.f90 +++ b/src/damage_nonlocal.f90 @@ -134,7 +134,7 @@ function damage_nonlocal_getDiffusion(ip,el) damage_nonlocal_getDiffusion = 0.0_pReal do grain = 1, homogenization_Ngrains(homog) damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + & - crystallite_push33ToRef(grain,ip,el,lattice_DamageDiffusion(1:3,1:3,material_phaseAt(grain,el))) + crystallite_push33ToRef(grain,ip,el,lattice_D(1:3,1:3,material_phaseAt(grain,el))) enddo damage_nonlocal_getDiffusion = & @@ -157,7 +157,7 @@ real(pReal) function damage_nonlocal_getMobility(ip,el) damage_nonlocal_getMobility = 0.0_pReal do ipc = 1, homogenization_Ngrains(material_homogenizationAt(el)) - damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_DamageMobility(material_phaseAt(ipc,el)) + damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_M(material_phaseAt(ipc,el)) enddo damage_nonlocal_getMobility = damage_nonlocal_getMobility/& diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90 index b226c3d0e..27df0acd5 100644 --- a/src/grid/discretization_grid.f90 +++ b/src/grid/discretization_grid.f90 @@ -107,9 +107,9 @@ subroutine discretization_grid_init(restart) call discretization_init(materialAt, & IPcoordinates0(myGrid,mySize,grid3Offset), & Nodes0(myGrid,mySize,grid3Offset),& - merge((grid(1)+1) * (grid(2)+1) * (grid3+1),& ! write bottom layer - (grid(1)+1) * (grid(2)+1) * grid3,& ! do not write bottom layer (is top of rank-1) - worldrank<1)) + merge((grid(1)+1) * (grid(2)+1) * (grid3+1),& ! write top layer... + (grid(1)+1) * (grid(2)+1) * grid3,& ! ...unless not last process + worldrank+1==worldsize)) FEsolving_execElem = [1,product(myGrid)] ! parallel loop bounds set to comprise all elements FEsolving_execIP = [1,1] ! parallel loop bounds set to comprise the only IP diff --git a/src/homogenization_mech_RGC.f90 b/src/homogenization_mech_RGC.f90 index f0485b244..6c015eea8 100644 --- a/src/homogenization_mech_RGC.f90 +++ b/src/homogenization_mech_RGC.f90 @@ -4,20 +4,20 @@ !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @brief Relaxed grain cluster (RGC) homogenization scheme -!> Nconstituents is defined as p x q x r (cluster) +!> N_constituents is defined as p x q x r (cluster) !-------------------------------------------------------------------------------------------------- submodule(homogenization) homogenization_mech_RGC use rotations type :: tParameters integer, dimension(:), allocatable :: & - Nconstituents + N_constituents real(pReal) :: & - xiAlpha, & - ciAlpha + xi_alpha, & + c_Alpha real(pReal), dimension(:), allocatable :: & - dAlpha, & - angles + D_alpha, & + a_g integer :: & of_debug = 0 character(len=pStringLen), allocatable, dimension(:) :: & @@ -163,20 +163,20 @@ module subroutine mech_RGC_init(num_homogMech) prm%output = homogMech%get_asStrings('output',defaultVal=emptyStringArray) #endif - prm%Nconstituents = homogMech%get_asInts('cluster_size',requiredSize=3) - if (homogenization_Ngrains(h) /= product(prm%Nconstituents)) & + prm%N_constituents = homogMech%get_asInts('cluster_size',requiredSize=3) + if (homogenization_Ngrains(h) /= product(prm%N_constituents)) & call IO_error(211,ext_msg='clustersize (mech_rgc)') - prm%xiAlpha = homogMech%get_asFloat('xi_alpha') - prm%ciAlpha = homogMech%get_asFloat('c_alpha') + prm%xi_alpha = homogMech%get_asFloat('xi_alpha') + prm%c_alpha = homogMech%get_asFloat('c_alpha') - prm%dAlpha = homogMech%get_asFloats('D_alpha', requiredSize=3) - prm%angles = homogMech%get_asFloats('a_g', requiredSize=3) + prm%D_alpha = homogMech%get_asFloats('D_alpha', requiredSize=3) + prm%a_g = homogMech%get_asFloats('a_g', requiredSize=3) NofMyHomog = count(material_homogenizationAt == h) - nIntFaceTot = 3*( (prm%Nconstituents(1)-1)*prm%Nconstituents(2)*prm%Nconstituents(3) & - + prm%Nconstituents(1)*(prm%Nconstituents(2)-1)*prm%Nconstituents(3) & - + prm%Nconstituents(1)*prm%Nconstituents(2)*(prm%Nconstituents(3)-1)) + nIntFaceTot = 3*( (prm%N_constituents(1)-1)*prm%N_constituents(2)*prm%N_constituents(3) & + + prm%N_constituents(1)*(prm%N_constituents(2)-1)*prm%N_constituents(3) & + + prm%N_constituents(1)*prm%N_constituents(2)*(prm%N_constituents(3)-1)) sizeState = nIntFaceTot & + size(['avg constitutive work ','average penalty energy']) @@ -197,8 +197,8 @@ module subroutine mech_RGC_init(num_homogMech) !-------------------------------------------------------------------------------------------------- ! assigning cluster orientations - dependentState(homogenization_typeInstance(h))%orientation = spread(eu2om(prm%angles*inRad),3,NofMyHomog) - !dst%orientation = spread(eu2om(prm%angles*inRad),3,NofMyHomog) ifort version 18.0.1 crashes (for whatever reason) + dependentState(homogenization_typeInstance(h))%orientation = spread(eu2om(prm%a_g*inRad),3,NofMyHomog) + !dst%orientation = spread(eu2om(prm%a_g*inRad),3,NofMyHomog) ifort version 18.0.1 crashes (for whatever reason) end associate @@ -229,8 +229,8 @@ module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of) !-------------------------------------------------------------------------------------------------- ! compute the deformation gradient of individual grains due to relaxations F = 0.0_pReal - do iGrain = 1,product(prm%Nconstituents) - iGrain3 = grain1to3(iGrain,prm%Nconstituents) + do iGrain = 1,product(prm%N_constituents) + iGrain3 = grain1to3(iGrain,prm%N_constituents) do iFace = 1,6 intFace = getInterface(iFace,iGrain3) ! identifying 6 interfaces of each grain aVect = relaxationVector(intFace,instance,of) ! get the relaxation vectors for each interface from global relaxation vector array @@ -290,7 +290,7 @@ module procedure mech_RGC_updateState !-------------------------------------------------------------------------------------------------- ! get the dimension of the cluster (grains and interfaces) - nGDim = prm%Nconstituents + nGDim = prm%N_constituents nGrain = product(nGDim) nIntFaceTot = (nGDim(1)-1)*nGDim(2)*nGDim(3) & + nGDim(1)*(nGDim(2)-1)*nGDim(3) & @@ -324,12 +324,12 @@ module procedure mech_RGC_updateState !------------------------------------------------------------------------------------------------ ! computing the residual stress from the balance of traction at all (interior) interfaces do iNum = 1,nIntFaceTot - faceID = interface1to4(iNum,param(instance)%Nconstituents) ! identifying the interface ID in local coordinate system (4-dimensional index) + faceID = interface1to4(iNum,param(instance)%N_constituents) ! identifying the interface ID in local coordinate system (4-dimensional index) !-------------------------------------------------------------------------------------------------- ! identify the left/bottom/back grain (-|N) iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate system (3-dimensional index) - iGrN = grain3to1(iGr3N,param(instance)%Nconstituents) ! translate the local grain ID into global coordinate system (1-dimensional index) + iGrN = grain3to1(iGr3N,param(instance)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) intFaceN = getInterface(2*faceID(1),iGr3N) normN = interfaceNormal(intFaceN,instance,of) @@ -337,7 +337,7 @@ module procedure mech_RGC_updateState ! identify the right/up/front grain (+|P) iGr3P = iGr3N iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate system (3-dimensional index) - iGrP = grain3to1(iGr3P,param(instance)%Nconstituents) ! translate the local grain ID into global coordinate system (1-dimensional index) + iGrP = grain3to1(iGr3P,param(instance)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) intFaceP = getInterface(2*faceID(1)-1,iGr3P) normP = interfaceNormal(intFaceP,instance,of) @@ -393,7 +393,7 @@ module procedure mech_RGC_updateState !-------------------------------------------------------------------------------------------------- ! compute/update the state for postResult, i.e., all energy densities computed by time-integration - do iGrain = 1,product(prm%Nconstituents) + do iGrain = 1,product(prm%N_constituents) do i = 1,3;do j = 1,3 stt%work(of) = stt%work(of) & + P(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/real(nGrain,pReal) @@ -450,18 +450,18 @@ module procedure mech_RGC_updateState ! ... of the constitutive stress tangent, assembled from dPdF or material constitutive model "smatrix" allocate(smatrix(3*nIntFaceTot,3*nIntFaceTot), source=0.0_pReal) do iNum = 1,nIntFaceTot - faceID = interface1to4(iNum,param(instance)%Nconstituents) ! assembling of local dPdF into global Jacobian matrix + faceID = interface1to4(iNum,param(instance)%N_constituents) ! assembling of local dPdF into global Jacobian matrix !-------------------------------------------------------------------------------------------------- ! identify the left/bottom/back grain (-|N) iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate sytem - iGrN = grain3to1(iGr3N,param(instance)%Nconstituents) ! translate into global grain ID + iGrN = grain3to1(iGr3N,param(instance)%N_constituents) ! translate into global grain ID intFaceN = getInterface(2*faceID(1),iGr3N) ! identifying the connecting interface in local coordinate system normN = interfaceNormal(intFaceN,instance,of) do iFace = 1,6 intFaceN = getInterface(iFace,iGr3N) ! identifying all interfaces that influence relaxation of the above interface mornN = interfaceNormal(intFaceN,instance,of) - iMun = interface4to1(intFaceN,param(instance)%Nconstituents) ! translate the interfaces ID into local 4-dimensional index + iMun = interface4to1(intFaceN,param(instance)%N_constituents) ! translate the interfaces ID into local 4-dimensional index if (iMun > 0) then ! get the corresponding tangent do i=1,3; do j=1,3; do k=1,3; do l=1,3 smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) & @@ -476,13 +476,13 @@ module procedure mech_RGC_updateState ! identify the right/up/front grain (+|P) iGr3P = iGr3N iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate sytem - iGrP = grain3to1(iGr3P,param(instance)%Nconstituents) ! translate into global grain ID + iGrP = grain3to1(iGr3P,param(instance)%N_constituents) ! translate into global grain ID intFaceP = getInterface(2*faceID(1)-1,iGr3P) ! identifying the connecting interface in local coordinate system normP = interfaceNormal(intFaceP,instance,of) do iFace = 1,6 intFaceP = getInterface(iFace,iGr3P) ! identifying all interfaces that influence relaxation of the above interface mornP = interfaceNormal(intFaceP,instance,of) - iMun = interface4to1(intFaceP,param(instance)%Nconstituents) ! translate the interfaces ID into local 4-dimensional index + iMun = interface4to1(intFaceP,param(instance)%N_constituents) ! translate the interfaces ID into local 4-dimensional index if (iMun > 0) then ! get the corresponding tangent do i=1,3; do j=1,3; do k=1,3; do l=1,3 smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) & @@ -522,12 +522,12 @@ module procedure mech_RGC_updateState ! computing the global stress residual array from the perturbed state p_resid = 0.0_pReal do iNum = 1,nIntFaceTot - faceID = interface1to4(iNum,param(instance)%Nconstituents) ! identifying the interface ID in local coordinate system (4-dimensional index) + faceID = interface1to4(iNum,param(instance)%N_constituents) ! identifying the interface ID in local coordinate system (4-dimensional index) !-------------------------------------------------------------------------------------------------- ! identify the left/bottom/back grain (-|N) iGr3N = faceID(2:4) ! identify the grain ID in local coordinate system (3-dimensional index) - iGrN = grain3to1(iGr3N,param(instance)%Nconstituents) ! translate the local grain ID into global coordinate system (1-dimensional index) + iGrN = grain3to1(iGr3N,param(instance)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) intFaceN = getInterface(2*faceID(1),iGr3N) ! identify the interface ID of the grain normN = interfaceNormal(intFaceN,instance,of) @@ -535,7 +535,7 @@ module procedure mech_RGC_updateState ! identify the right/up/front grain (+|P) iGr3P = iGr3N iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identify the grain ID in local coordinate system (3-dimensional index) - iGrP = grain3to1(iGr3P,param(instance)%Nconstituents) ! translate the local grain ID into global coordinate system (1-dimensional index) + iGrP = grain3to1(iGr3P,param(instance)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) intFaceP = getInterface(2*faceID(1)-1,iGr3P) ! identify the interface ID of the grain normP = interfaceNormal(intFaceP,instance,of) @@ -664,7 +664,7 @@ module procedure mech_RGC_updateState real(pReal) :: muGrain,muGNghb,nDefNorm,bgGrain,bgGNghb real(pReal), parameter :: nDefToler = 1.0e-10_pReal - nGDim = param(instance)%Nconstituents + nGDim = param(instance)%N_constituents rPen = 0.0_pReal nMis = 0.0_pReal @@ -685,11 +685,11 @@ module procedure mech_RGC_updateState !----------------------------------------------------------------------------------------------- ! computing the mismatch and penalty stress tensor of all grains - grainLoop: do iGrain = 1,product(prm%Nconstituents) + grainLoop: do iGrain = 1,product(prm%N_constituents) Gmoduli = equivalentModuli(iGrain,ip,el) muGrain = Gmoduli(1) ! collecting the equivalent shear modulus of grain bgGrain = Gmoduli(2) ! and the lengthh of Burgers vector - iGrain3 = grain1to3(iGrain,prm%Nconstituents) ! get the grain ID in local 3-dimensional index (x,y,z)-position + iGrain3 = grain1to3(iGrain,prm%N_constituents) ! get the grain ID in local 3-dimensional index (x,y,z)-position interfaceLoop: do iFace = 1,6 intFace = getInterface(iFace,iGrain3) ! get the 4-dimensional index of the interface in local numbering system of the grain @@ -699,7 +699,7 @@ module procedure mech_RGC_updateState + int(real(intFace(1),pReal)/real(abs(intFace(1)),pReal)) where(iGNghb3 < 1) iGNghb3 = nGDim where(iGNghb3 >nGDim) iGNghb3 = 1 - iGNghb = grain3to1(iGNghb3,prm%Nconstituents) ! get the ID of the neighboring grain + iGNghb = grain3to1(iGNghb3,prm%N_constituents) ! get the ID of the neighboring grain Gmoduli = equivalentModuli(iGNghb,ip,el) ! collect the shear modulus and Burgers vector of the neighbor muGNghb = Gmoduli(1) bgGNghb = Gmoduli(2) @@ -728,9 +728,9 @@ module procedure mech_RGC_updateState !------------------------------------------------------------------------------------------- ! compute the stress penalty of all interfaces do i = 1,3; do j = 1,3; do k = 1,3; do l = 1,3 - rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*bgGrain + muGNghb*bgGNghb)*prm%xiAlpha & - *surfCorr(abs(intFace(1)))/prm%dAlpha(abs(intFace(1))) & - *cosh(prm%ciAlpha*nDefNorm) & + rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*bgGrain + muGNghb*bgGNghb)*prm%xi_alpha & + *surfCorr(abs(intFace(1)))/prm%D_alpha(abs(intFace(1))) & + *cosh(prm%c_alpha*nDefNorm) & *0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_LeviCivita(k,l,j) & *tanh(nDefNorm/num%xSmoo) enddo; enddo;enddo; enddo @@ -885,8 +885,8 @@ module procedure mech_RGC_updateState associate(prm => param(instance)) F = 0.0_pReal - do iGrain = 1,product(prm%Nconstituents) - iGrain3 = grain1to3(iGrain,prm%Nconstituents) + do iGrain = 1,product(prm%N_constituents) + iGrain3 = grain1to3(iGrain,prm%N_constituents) do iFace = 1,6 intFace = getInterface(iFace,iGrain3) aVect = relaxationVector(intFace,instance,of) @@ -916,8 +916,8 @@ module subroutine mech_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ins real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses integer, intent(in) :: instance - avgP = sum(P,3) /real(product(param(instance)%Nconstituents),pReal) - dAvgPdAvgF = sum(dPdF,5)/real(product(param(instance)%Nconstituents),pReal) + avgP = sum(P,3) /real(product(param(instance)%N_constituents),pReal) + dAvgPdAvgF = sum(dPdF,5)/real(product(param(instance)%N_constituents),pReal) end subroutine mech_RGC_averageStressAndItsTangent @@ -975,7 +975,7 @@ pure function relaxationVector(intFace,instance,of) !-------------------------------------------------------------------------------------------------- ! collect the interface relaxation vector from the global state array - iNum = interface4to1(intFace,param(instance)%Nconstituents) ! identify the position of the interface in global state array + iNum = interface4to1(intFace,param(instance)%N_constituents) ! identify the position of the interface in global state array if (iNum > 0) then relaxationVector = state(instance)%relaxationVector((3*iNum-2):(3*iNum),of) else diff --git a/src/homogenization_mech_isostrain.f90 b/src/homogenization_mech_isostrain.f90 index 5138afa73..6883f1c37 100644 --- a/src/homogenization_mech_isostrain.f90 +++ b/src/homogenization_mech_isostrain.f90 @@ -13,7 +13,7 @@ submodule(homogenization) homogenization_mech_isostrain type :: tParameters !< container type for internal constitutive parameters integer :: & - Nconstituents + N_constituents integer(kind(average_ID)) :: & mapping end type @@ -51,7 +51,7 @@ module subroutine mech_isostrain_init homogMech => homog%get('mech') associate(prm => param(homogenization_typeInstance(h))) - prm%Nconstituents = homogMech%get_asInt('N_constituents') + prm%N_constituents = homogMech%get_asInt('N_constituents') select case(homogMech%get_asString('mapping',defaultVal = 'sum')) case ('sum') prm%mapping = parallel_ID @@ -107,8 +107,8 @@ module subroutine mech_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dP avgP = sum(P,3) dAvgPdAvgF = sum(dPdF,5) case (average_ID) - avgP = sum(P,3) /real(prm%Nconstituents,pReal) - dAvgPdAvgF = sum(dPdF,5)/real(prm%Nconstituents,pReal) + avgP = sum(P,3) /real(prm%N_constituents,pReal) + dAvgPdAvgF = sum(dPdF,5)/real(prm%N_constituents,pReal) end select end associate diff --git a/src/kinematics_cleavage_opening.f90 b/src/kinematics_cleavage_opening.f90 index 23f348831..d52fdbc1c 100644 --- a/src/kinematics_cleavage_opening.f90 +++ b/src/kinematics_cleavage_opening.f90 @@ -12,10 +12,10 @@ submodule(constitutive:constitutive_damage) kinematics_cleavage_opening integer :: & sum_N_cl !< total number of cleavage planes real(pReal) :: & - sdot0, & !< opening rate of cleavage planes - n !< damage rate sensitivity + dot_o, & !< opening rate of cleavage planes + q !< damage rate sensitivity real(pReal), dimension(:), allocatable :: & - critLoad + g_crit real(pReal), dimension(:,:,:,:), allocatable :: & cleavage_systems end type tParameters @@ -70,21 +70,21 @@ module function kinematics_cleavage_opening_init(kinematics_length) result(myKin N_cl = kinematic_type%get_asInts('N_cl') prm%sum_N_cl = sum(abs(N_cl)) - prm%n = kinematic_type%get_asFloat('q') - prm%sdot0 = kinematic_type%get_asFloat('dot_o') + prm%q = kinematic_type%get_asFloat('q') + prm%dot_o = kinematic_type%get_asFloat('dot_o') - prm%critLoad = kinematic_type%get_asFloats('g_crit',requiredSize=size(N_cl)) + prm%g_crit = kinematic_type%get_asFloats('g_crit',requiredSize=size(N_cl)) prm%cleavage_systems = lattice_SchmidMatrix_cleavage(N_cl,phase%get_asString('lattice'),& phase%get_asFloat('c/a',defaultVal=0.0_pReal)) ! expand: family => system - prm%critLoad = math_expand(prm%critLoad,N_cl) + prm%g_crit = math_expand(prm%g_crit,N_cl) ! sanity checks - if (prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' q' - if (prm%sdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' dot_o' - if (any(prm%critLoad < 0.0_pReal)) extmsg = trim(extmsg)//' g_crit' + if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' + if (prm%dot_o <= 0.0_pReal) extmsg = trim(extmsg)//' dot_o' + if (any(prm%g_crit < 0.0_pReal)) extmsg = trim(extmsg)//' g_crit' !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range @@ -128,13 +128,13 @@ module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, dLd_dTstar = 0.0_pReal associate(prm => param(kinematics_cleavage_opening_instance(material_phaseAt(ipc,el)))) do i = 1,prm%sum_N_cl - traction_crit = prm%critLoad(i)* damage(homog)%p(damageOffset)**2.0_pReal + traction_crit = prm%g_crit(i)* damage(homog)%p(damageOffset)**2.0_pReal traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i)) if (abs(traction_d) > traction_crit + tol_math_check) then - udotd = sign(1.0_pReal,traction_d)* prm%sdot0 * ((abs(traction_d) - traction_crit)/traction_crit)**prm%n + udotd = sign(1.0_pReal,traction_d)* prm%dot_o * ((abs(traction_d) - traction_crit)/traction_crit)**prm%q Ld = Ld + udotd*prm%cleavage_systems(1:3,1:3,1,i) - dudotd_dt = sign(1.0_pReal,traction_d)*udotd*prm%n / (abs(traction_d) - traction_crit) + dudotd_dt = sign(1.0_pReal,traction_d)*udotd*prm%q / (abs(traction_d) - traction_crit) forall (k=1:3,l=1:3,m=1:3,n=1:3) & dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) & + dudotd_dt*prm%cleavage_systems(k,l,1,i) * prm%cleavage_systems(m,n,1,i) @@ -142,9 +142,9 @@ module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i)) if (abs(traction_t) > traction_crit + tol_math_check) then - udott = sign(1.0_pReal,traction_t)* prm%sdot0 * ((abs(traction_t) - traction_crit)/traction_crit)**prm%n + udott = sign(1.0_pReal,traction_t)* prm%dot_o * ((abs(traction_t) - traction_crit)/traction_crit)**prm%q Ld = Ld + udott*prm%cleavage_systems(1:3,1:3,2,i) - dudott_dt = sign(1.0_pReal,traction_t)*udott*prm%n / (abs(traction_t) - traction_crit) + dudott_dt = sign(1.0_pReal,traction_t)*udott*prm%q / (abs(traction_t) - traction_crit) forall (k=1:3,l=1:3,m=1:3,n=1:3) & dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) & + dudott_dt*prm%cleavage_systems(k,l,2,i) * prm%cleavage_systems(m,n,2,i) @@ -152,9 +152,9 @@ module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i)) if (abs(traction_n) > traction_crit + tol_math_check) then - udotn = sign(1.0_pReal,traction_n)* prm%sdot0 * ((abs(traction_n) - traction_crit)/traction_crit)**prm%n + udotn = sign(1.0_pReal,traction_n)* prm%dot_o * ((abs(traction_n) - traction_crit)/traction_crit)**prm%q Ld = Ld + udotn*prm%cleavage_systems(1:3,1:3,3,i) - dudotn_dt = sign(1.0_pReal,traction_n)*udotn*prm%n / (abs(traction_n) - traction_crit) + dudotn_dt = sign(1.0_pReal,traction_n)*udotn*prm%q / (abs(traction_n) - traction_crit) forall (k=1:3,l=1:3,m=1:3,n=1:3) & dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) & + dudotn_dt*prm%cleavage_systems(k,l,3,i) * prm%cleavage_systems(m,n,3,i) diff --git a/src/kinematics_slipplane_opening.f90 b/src/kinematics_slipplane_opening.f90 index 660483b90..e0de5e181 100644 --- a/src/kinematics_slipplane_opening.f90 +++ b/src/kinematics_slipplane_opening.f90 @@ -12,10 +12,10 @@ submodule(constitutive:constitutive_damage) kinematics_slipplane_opening integer :: & sum_N_sl !< total number of cleavage planes real(pReal) :: & - sdot0, & !< opening rate of cleavage planes - n !< damage rate sensitivity + dot_o, & !< opening rate of cleavage planes + q !< damage rate sensitivity real(pReal), dimension(:), allocatable :: & - critLoad + g_crit real(pReal), dimension(:,:,:), allocatable :: & P_d, & P_t, & @@ -70,8 +70,8 @@ module function kinematics_slipplane_opening_init(kinematics_length) result(myKi associate(prm => param(kinematics_slipplane_opening_instance(p))) kinematic_type => kinematics%get(k) - prm%sdot0 = kinematic_type%get_asFloat('dot_o') - prm%n = kinematic_type%get_asFloat('q') + prm%dot_o = kinematic_type%get_asFloat('dot_o') + prm%q = kinematic_type%get_asFloat('q') N_sl = pl%get_asInts('N_sl') prm%sum_N_sl = sum(abs(N_sl)) @@ -89,15 +89,15 @@ module function kinematics_slipplane_opening_init(kinematics_length) result(myKi prm%P_n(1:3,1:3,i) = math_outer(n(1:3,i), n(1:3,i)) enddo - prm%critLoad = kinematic_type%get_asFloats('g_crit',requiredSize=size(N_sl)) + prm%g_crit = kinematic_type%get_asFloats('g_crit',requiredSize=size(N_sl)) ! expand: family => system - prm%critLoad = math_expand(prm%critLoad,N_sl) + prm%g_crit = math_expand(prm%g_crit,N_sl) ! sanity checks - if (prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' anisoDuctile_n' - if (prm%sdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' anisoDuctile_sdot0' - if (any(prm%critLoad < 0.0_pReal)) extmsg = trim(extmsg)//' anisoDuctile_critLoad' + if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' anisoDuctile_n' + if (prm%dot_o <= 0.0_pReal) extmsg = trim(extmsg)//' anisoDuctile_sdot0' + if (any(prm%g_crit < 0.0_pReal)) extmsg = trim(extmsg)//' anisoDuctile_critLoad' !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range @@ -150,27 +150,27 @@ module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S traction_t = math_tensordot(S,prm%P_t(1:3,1:3,i)) traction_n = math_tensordot(S,prm%P_n(1:3,1:3,i)) - traction_crit = prm%critLoad(i)* damage(homog)%p(damageOffset) ! degrading critical load carrying capacity by damage + traction_crit = prm%g_crit(i)* damage(homog)%p(damageOffset) ! degrading critical load carrying capacity by damage - udotd = sign(1.0_pReal,traction_d)* prm%sdot0* ( abs(traction_d)/traction_crit & - - abs(traction_d)/prm%critLoad(i))**prm%n - udott = sign(1.0_pReal,traction_t)* prm%sdot0* ( abs(traction_t)/traction_crit & - - abs(traction_t)/prm%critLoad(i))**prm%n - udotn = prm%sdot0* ( max(0.0_pReal,traction_n)/traction_crit & - - max(0.0_pReal,traction_n)/prm%critLoad(i))**prm%n + udotd = sign(1.0_pReal,traction_d)* prm%dot_o* ( abs(traction_d)/traction_crit & + - abs(traction_d)/prm%g_crit(i))**prm%q + udott = sign(1.0_pReal,traction_t)* prm%dot_o* ( abs(traction_t)/traction_crit & + - abs(traction_t)/prm%g_crit(i))**prm%q + udotn = prm%dot_o* ( max(0.0_pReal,traction_n)/traction_crit & + - max(0.0_pReal,traction_n)/prm%g_crit(i))**prm%q if (dNeq0(traction_d)) then - dudotd_dt = udotd*prm%n/traction_d + dudotd_dt = udotd*prm%q/traction_d else dudotd_dt = 0.0_pReal endif if (dNeq0(traction_t)) then - dudott_dt = udott*prm%n/traction_t + dudott_dt = udott*prm%q/traction_t else dudott_dt = 0.0_pReal endif if (dNeq0(traction_n)) then - dudotn_dt = udotn*prm%n/traction_n + dudotn_dt = udotn*prm%q/traction_n else dudotn_dt = 0.0_pReal endif diff --git a/src/kinematics_thermal_expansion.f90 b/src/kinematics_thermal_expansion.f90 index 93a48e035..772f5abbf 100644 --- a/src/kinematics_thermal_expansion.f90 +++ b/src/kinematics_thermal_expansion.f90 @@ -11,7 +11,7 @@ submodule(constitutive:constitutive_thermal) kinematics_thermal_expansion real(pReal) :: & T_ref real(pReal), dimension(3,3,3) :: & - expansion = 0.0_pReal + A = 0.0_pReal end type tParameters type(tParameters), dimension(:), allocatable :: param @@ -64,13 +64,13 @@ module function kinematics_thermal_expansion_init(kinematics_length) result(myKi ! read up to three parameters (constant, linear, quadratic with T) temp = kinematic_type%get_asFloats('A_11') - prm%expansion(1,1,1:size(temp)) = temp + prm%A(1,1,1:size(temp)) = temp temp = kinematic_type%get_asFloats('A_22',defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(temp)) - prm%expansion(2,2,1:size(temp)) = temp + prm%A(2,2,1:size(temp)) = temp temp = kinematic_type%get_asFloats('A_33',defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(temp)) - prm%expansion(3,3,1:size(temp)) = temp - do i=1, size(prm%expansion,3) - prm%expansion(1:3,1:3,i) = lattice_applyLatticeSymmetry33(prm%expansion(1:3,1:3,i),& + prm%A(3,3,1:size(temp)) = temp + do i=1, size(prm%A,3) + prm%A(1:3,1:3,i) = lattice_applyLatticeSymmetry33(prm%A(1:3,1:3,i),& phase%get_asString('lattice')) enddo @@ -94,13 +94,13 @@ pure module function kinematics_thermal_expansion_initialStrain(homog,phase,offs offset real(pReal), dimension(3,3) :: & - initialStrain !< initial thermal strain (should be small strain, though) + initialStrain !< initial thermal strain (should be small strain, though) associate(prm => param(kinematics_thermal_expansion_instance(phase))) initialStrain = & - (temperature(homog)%p(offset) - prm%T_ref)**1 / 1. * prm%expansion(1:3,1:3,1) + & ! constant coefficient - (temperature(homog)%p(offset) - prm%T_ref)**2 / 2. * prm%expansion(1:3,1:3,2) + & ! linear coefficient - (temperature(homog)%p(offset) - prm%T_ref)**3 / 3. * prm%expansion(1:3,1:3,3) ! quadratic coefficient + (temperature(homog)%p(offset) - prm%T_ref)**1 / 1. * prm%A(1:3,1:3,1) + & ! constant coefficient + (temperature(homog)%p(offset) - prm%T_ref)**2 / 2. * prm%A(1:3,1:3,2) + & ! linear coefficient + (temperature(homog)%p(offset) - prm%T_ref)**3 / 3. * prm%A(1:3,1:3,3) ! quadratic coefficient end associate end function kinematics_thermal_expansion_initialStrain @@ -133,14 +133,14 @@ module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, i associate(prm => param(kinematics_thermal_expansion_instance(phase))) Li = TDot * ( & - prm%expansion(1:3,1:3,1)*(T - prm%T_ref)**0 & ! constant coefficient - + prm%expansion(1:3,1:3,2)*(T - prm%T_ref)**1 & ! linear coefficient - + prm%expansion(1:3,1:3,3)*(T - prm%T_ref)**2 & ! quadratic coefficient + prm%A(1:3,1:3,1)*(T - prm%T_ref)**0 & ! constant coefficient + + prm%A(1:3,1:3,2)*(T - prm%T_ref)**1 & ! linear coefficient + + prm%A(1:3,1:3,3)*(T - prm%T_ref)**2 & ! quadratic coefficient ) / & (1.0_pReal & - + prm%expansion(1:3,1:3,1)*(T - prm%T_ref)**1 / 1. & - + prm%expansion(1:3,1:3,2)*(T - prm%T_ref)**2 / 2. & - + prm%expansion(1:3,1:3,3)*(T - prm%T_ref)**3 / 3. & + + prm%A(1:3,1:3,1)*(T - prm%T_ref)**1 / 1. & + + prm%A(1:3,1:3,2)*(T - prm%T_ref)**2 / 2. & + + prm%A(1:3,1:3,3)*(T - prm%T_ref)**3 / 3. & ) end associate dLi_dTstar = 0.0_pReal diff --git a/src/lattice.f90 b/src/lattice.f90 index 69305c839..78b3894c0 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -394,13 +394,13 @@ module lattice ! SHOULD NOT BE PART OF LATTICE BEGIN real(pReal), dimension(:), allocatable, public, protected :: & lattice_mu, lattice_nu, & - lattice_damageMobility, & - lattice_massDensity, & - lattice_specificHeat + lattice_M, & + lattice_rho, & + lattice_c_p real(pReal), dimension(:,:,:), allocatable, public, protected :: & lattice_C66, & - lattice_thermalConductivity, & - lattice_damageDiffusion + lattice_K, & + lattice_D integer(kind(lattice_UNDEFINED_ID)), dimension(:), allocatable, public, protected :: & lattice_structure ! SHOULD NOT BE PART OF LATTICE END @@ -465,11 +465,11 @@ subroutine lattice_init allocate(lattice_structure(Nphases),source = lattice_UNDEFINED_ID) allocate(lattice_C66(6,6,Nphases), source=0.0_pReal) - allocate(lattice_thermalConductivity (3,3,Nphases), source=0.0_pReal) - allocate(lattice_damageDiffusion (3,3,Nphases), source=0.0_pReal) + allocate(lattice_K (3,3,Nphases), source=0.0_pReal) + allocate(lattice_D (3,3,Nphases), source=0.0_pReal) - allocate(lattice_damageMobility,& - lattice_massDensity,lattice_specificHeat, & + allocate(lattice_M,& + lattice_rho,lattice_c_p, & lattice_mu, lattice_nu,& source=[(0.0_pReal,i=1,Nphases)]) @@ -517,22 +517,22 @@ subroutine lattice_init ! SHOULD NOT BE PART OF LATTICE BEGIN - lattice_thermalConductivity(1,1,p) = phase%get_asFloat('K_11',defaultVal=0.0_pReal) - lattice_thermalConductivity(2,2,p) = phase%get_asFloat('K_22',defaultVal=0.0_pReal) - lattice_thermalConductivity(3,3,p) = phase%get_asFloat('K_33',defaultVal=0.0_pReal) - lattice_thermalConductivity(1:3,1:3,p) = lattice_applyLatticeSymmetry33(lattice_thermalConductivity(1:3,1:3,p), & - phase%get_asString('lattice')) + lattice_K(1,1,p) = phase%get_asFloat('K_11',defaultVal=0.0_pReal) + lattice_K(2,2,p) = phase%get_asFloat('K_22',defaultVal=0.0_pReal) + lattice_K(3,3,p) = phase%get_asFloat('K_33',defaultVal=0.0_pReal) + lattice_K(1:3,1:3,p) = lattice_applyLatticeSymmetry33(lattice_K(1:3,1:3,p), & + phase%get_asString('lattice')) - lattice_specificHeat(p) = phase%get_asFloat('c_p',defaultVal=0.0_pReal) - lattice_massDensity(p) = phase%get_asFloat('rho', defaultVal=0.0_pReal) + lattice_c_p(p) = phase%get_asFloat('c_p', defaultVal=0.0_pReal) + lattice_rho(p) = phase%get_asFloat('rho', defaultVal=0.0_pReal) - lattice_DamageDiffusion(1,1,p) = phase%get_asFloat('D_11',defaultVal=0.0_pReal) - lattice_DamageDiffusion(2,2,p) = phase%get_asFloat('D_22',defaultVal=0.0_pReal) - lattice_DamageDiffusion(3,3,p) = phase%get_asFloat('D_33',defaultVal=0.0_pReal) - lattice_DamageDiffusion(1:3,1:3,p) = lattice_applyLatticeSymmetry33(lattice_DamageDiffusion(1:3,1:3,p), & - phase%get_asString('lattice')) + lattice_D(1,1,p) = phase%get_asFloat('D_11',defaultVal=0.0_pReal) + lattice_D(2,2,p) = phase%get_asFloat('D_22',defaultVal=0.0_pReal) + lattice_D(3,3,p) = phase%get_asFloat('D_33',defaultVal=0.0_pReal) + lattice_D(1:3,1:3,p) = lattice_applyLatticeSymmetry33(lattice_D(1:3,1:3,p), & + phase%get_asString('lattice')) - lattice_DamageMobility(p) = phase%get_asFloat('M',defaultVal=0.0_pReal) + lattice_M(p) = phase%get_asFloat('M',defaultVal=0.0_pReal) ! SHOULD NOT BE PART OF LATTICE END call selfTest diff --git a/src/source_damage_anisoBrittle.f90 b/src/source_damage_anisoBrittle.f90 index 7911d6d0a..1954e6ea3 100644 --- a/src/source_damage_anisoBrittle.f90 +++ b/src/source_damage_anisoBrittle.f90 @@ -12,11 +12,11 @@ submodule (constitutive:constitutive_damage) source_damage_anisoBrittle type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & - sdot_0, & !< opening rate of cleavage planes - n !< damage rate sensitivity + dot_o, & !< opening rate of cleavage planes + q !< damage rate sensitivity real(pReal), dimension(:), allocatable :: & - critDisp, & !< critical displacement - critLoad !< critical load + s_crit, & !< critical displacement + g_crit !< critical load real(pReal), dimension(:,:,:,:), allocatable :: & cleavage_systems integer :: & @@ -75,18 +75,18 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources) N_cl = src%get_asInts('N_cl',defaultVal=emptyIntArray) prm%sum_N_cl = sum(abs(N_cl)) - prm%n = src%get_asFloat('q') - prm%sdot_0 = src%get_asFloat('dot_o') + prm%q = src%get_asFloat('q') + prm%dot_o = src%get_asFloat('dot_o') - prm%critDisp = src%get_asFloats('s_crit', requiredSize=size(N_cl)) - prm%critLoad = src%get_asFloats('g_crit', requiredSize=size(N_cl)) + prm%s_crit = src%get_asFloats('s_crit', requiredSize=size(N_cl)) + prm%g_crit = src%get_asFloats('g_crit', requiredSize=size(N_cl)) prm%cleavage_systems = lattice_SchmidMatrix_cleavage(N_cl,phase%get_asString('lattice'),& phase%get_asFloat('c/a',defaultVal=0.0_pReal)) ! expand: family => system - prm%critDisp = math_expand(prm%critDisp,N_cl) - prm%critLoad = math_expand(prm%critLoad,N_cl) + prm%s_crit = math_expand(prm%s_crit,N_cl) + prm%g_crit = math_expand(prm%g_crit,N_cl) #if defined (__GFORTRAN__) prm%output = output_asStrings(src) @@ -95,10 +95,10 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources) #endif ! sanity checks - if (prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' q' - if (prm%sdot_0 <= 0.0_pReal) extmsg = trim(extmsg)//' dot_o' - if (any(prm%critLoad < 0.0_pReal)) extmsg = trim(extmsg)//' g_crit' - if (any(prm%critDisp < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit' + if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' + if (prm%dot_o <= 0.0_pReal) extmsg = trim(extmsg)//' dot_o' + if (any(prm%g_crit < 0.0_pReal)) extmsg = trim(extmsg)//' g_crit' + if (any(prm%s_crit < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit' NipcMyPhase = count(material_phaseAt==p) * discretization_nIP call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0) @@ -152,14 +152,14 @@ module subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i)) traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i)) - traction_crit = prm%critLoad(i)*damage(homog)%p(damageOffset)**2.0_pReal + traction_crit = prm%g_crit(i)*damage(homog)%p(damageOffset)**2.0_pReal sourceState(phase)%p(sourceOffset)%dotState(1,constituent) & = sourceState(phase)%p(sourceOffset)%dotState(1,constituent) & - + prm%sdot_0 / prm%critDisp(i) & - * ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%n + & - (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%n + & - (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**prm%n) + + prm%dot_o / prm%s_crit(i) & + * ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%q + & + (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%q + & + (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**prm%q) enddo end associate diff --git a/src/source_damage_anisoDuctile.f90 b/src/source_damage_anisoDuctile.f90 index 52189c839..3b8eb2428 100644 --- a/src/source_damage_anisoDuctile.f90 +++ b/src/source_damage_anisoDuctile.f90 @@ -12,9 +12,9 @@ submodule(constitutive:constitutive_damage) source_damage_anisoDuctile type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & - n !< damage rate sensitivity + q !< damage rate sensitivity real(pReal), dimension(:), allocatable :: & - critPlasticStrain !< critical plastic strain per slip system + gamma_crit !< critical plastic strain per slip system character(len=pStringLen), allocatable, dimension(:) :: & output end type tParameters @@ -67,12 +67,12 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources) associate(prm => param(source_damage_anisoDuctile_instance(p))) src => sources%get(sourceOffset) - N_sl = pl%get_asInts('N_sl',defaultVal=emptyIntArray) - prm%n = src%get_asFloat('q') - prm%critPlasticStrain = src%get_asFloats('gamma_crit',requiredSize=size(N_sl)) + N_sl = pl%get_asInts('N_sl',defaultVal=emptyIntArray) + prm%q = src%get_asFloat('q') + prm%gamma_crit = src%get_asFloats('gamma_crit',requiredSize=size(N_sl)) ! expand: family => system - prm%critPlasticStrain = math_expand(prm%critPlasticStrain,N_sl) + prm%gamma_crit = math_expand(prm%gamma_crit,N_sl) #if defined (__GFORTRAN__) prm%output = output_asStrings(src) @@ -81,8 +81,8 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources) #endif ! sanity checks - if (prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' q' - if (any(prm%critPlasticStrain < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit' + if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' + if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit' NipcMyPhase=count(material_phaseAt==p) * discretization_nIP call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0) @@ -127,7 +127,7 @@ module subroutine source_damage_anisoDuctile_dotState(ipc, ip, el) associate(prm => param(source_damage_anisoDuctile_instance(phase))) sourceState(phase)%p(sourceOffset)%dotState(1,constituent) & - = sum(plasticState(phase)%slipRate(:,constituent)/(damage(homog)%p(damageOffset)**prm%n)/prm%critPlasticStrain) + = sum(plasticState(phase)%slipRate(:,constituent)/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit) end associate end subroutine source_damage_anisoDuctile_dotState diff --git a/src/source_damage_isoDuctile.f90 b/src/source_damage_isoDuctile.f90 index 493183d75..d958aed6a 100644 --- a/src/source_damage_isoDuctile.f90 +++ b/src/source_damage_isoDuctile.f90 @@ -12,8 +12,8 @@ submodule (constitutive:constitutive_damage) source_damage_isoDuctile type:: tParameters !< container type for internal constitutive parameters real(pReal) :: & - critPlasticStrain, & !< critical plastic strain - N + gamma_crit, & !< critical plastic strain + q character(len=pStringLen), allocatable, dimension(:) :: & output end type tParameters @@ -64,8 +64,8 @@ module function source_damage_isoDuctile_init(source_length) result(mySources) associate(prm => param(source_damage_isoDuctile_instance(p))) src => sources%get(sourceOffset) - prm%N = src%get_asFloat('q') - prm%critPlasticStrain = src%get_asFloat('gamma_crit') + prm%q = src%get_asFloat('q') + prm%gamma_crit = src%get_asFloat('gamma_crit') #if defined (__GFORTRAN__) prm%output = output_asStrings(src) @@ -74,8 +74,8 @@ module function source_damage_isoDuctile_init(source_length) result(mySources) #endif ! sanity checks - if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' q' - if (prm%critPlasticStrain <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit' + if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' + if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit' NipcMyPhase=count(material_phaseAt==p) * discretization_nIP call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0) @@ -120,7 +120,7 @@ module subroutine source_damage_isoDuctile_dotState(ipc, ip, el) associate(prm => param(source_damage_isoDuctile_instance(phase))) sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = & - sum(plasticState(phase)%slipRate(:,constituent))/(damage(homog)%p(damageOffset)**prm%N)/prm%critPlasticStrain + sum(plasticState(phase)%slipRate(:,constituent))/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit end associate end subroutine source_damage_isoDuctile_dotState diff --git a/src/source_thermal_externalheat.f90 b/src/source_thermal_externalheat.f90 index 4a644f53b..2eeeb47df 100644 --- a/src/source_thermal_externalheat.f90 +++ b/src/source_thermal_externalheat.f90 @@ -13,8 +13,8 @@ submodule(constitutive:constitutive_thermal) source_thermal_externalheat type :: tParameters !< container type for internal constitutive parameters real(pReal), dimension(:), allocatable :: & - time, & - heat_rate + t_n, & + f_T integer :: & nIntervals end type tParameters @@ -64,10 +64,10 @@ module function source_thermal_externalheat_init(source_length) result(mySources associate(prm => param(source_thermal_externalheat_instance(p))) src => sources%get(sourceOffset) - prm%time = src%get_asFloats('t_n') - prm%nIntervals = size(prm%time) - 1 + prm%t_n = src%get_asFloats('t_n') + prm%nIntervals = size(prm%t_n) - 1 - prm%heat_rate = src%get_asFloats('f_T',requiredSize = size(prm%time)) + prm%f_T = src%get_asFloats('f_T',requiredSize = size(prm%t_n)) NipcMyPhase = count(material_phaseAt==p) * discretization_nIP call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0) @@ -121,13 +121,13 @@ module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_d associate(prm => param(source_thermal_externalheat_instance(phase))) do interval = 1, prm%nIntervals ! scan through all rate segments - frac_time = (sourceState(phase)%p(sourceOffset)%state(1,of) - prm%time(interval)) & - / (prm%time(interval+1) - prm%time(interval)) ! fractional time within segment + frac_time = (sourceState(phase)%p(sourceOffset)%state(1,of) - prm%t_n(interval)) & + / (prm%t_n(interval+1) - prm%t_n(interval)) ! fractional time within segment if ( (frac_time < 0.0_pReal .and. interval == 1) & .or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) & .or. (frac_time >= 0.0_pReal .and. frac_time < 1.0_pReal) ) & - TDot = prm%heat_rate(interval ) * (1.0_pReal - frac_time) + & - prm%heat_rate(interval+1) * frac_time ! interpolate heat rate between segment boundaries... + TDot = prm%f_T(interval ) * (1.0_pReal - frac_time) + & + prm%f_T(interval+1) * frac_time ! interpolate heat rate between segment boundaries... ! ...or extrapolate if outside of bounds enddo dTDot_dT = 0.0 diff --git a/src/thermal_adiabatic.f90 b/src/thermal_adiabatic.f90 index 1601a9b4f..63deb3cd5 100644 --- a/src/thermal_adiabatic.f90 +++ b/src/thermal_adiabatic.f90 @@ -169,7 +169,7 @@ function thermal_adiabatic_getSpecificHeat(ip,el) do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) thermal_adiabatic_getSpecificHeat = thermal_adiabatic_getSpecificHeat & - + lattice_specificHeat(material_phaseAt(grain,el)) + + lattice_c_p(material_phaseAt(grain,el)) enddo thermal_adiabatic_getSpecificHeat = thermal_adiabatic_getSpecificHeat & @@ -195,7 +195,7 @@ function thermal_adiabatic_getMassDensity(ip,el) do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) thermal_adiabatic_getMassDensity = thermal_adiabatic_getMassDensity & - + lattice_massDensity(material_phaseAt(grain,el)) + + lattice_rho(material_phaseAt(grain,el)) enddo thermal_adiabatic_getMassDensity = thermal_adiabatic_getMassDensity & diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index 5556a54e1..69ce8025e 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -127,7 +127,7 @@ function thermal_conduction_getConductivity(ip,el) thermal_conduction_getConductivity = 0.0_pReal do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) thermal_conduction_getConductivity = thermal_conduction_getConductivity + & - crystallite_push33ToRef(grain,ip,el,lattice_thermalConductivity(:,:,material_phaseAt(grain,el))) + crystallite_push33ToRef(grain,ip,el,lattice_K(:,:,material_phaseAt(grain,el))) enddo thermal_conduction_getConductivity = thermal_conduction_getConductivity & @@ -153,7 +153,7 @@ function thermal_conduction_getSpecificHeat(ip,el) do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat & - + lattice_specificHeat(material_phaseAt(grain,el)) + + lattice_c_p(material_phaseAt(grain,el)) enddo thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat & @@ -180,7 +180,7 @@ function thermal_conduction_getMassDensity(ip,el) do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) thermal_conduction_getMassDensity = thermal_conduction_getMassDensity & - + lattice_massDensity(material_phaseAt(grain,el)) + + lattice_rho(material_phaseAt(grain,el)) enddo thermal_conduction_getMassDensity = thermal_conduction_getMassDensity &