223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Merge remote-tracking branch 'origin/development' into restructure-numerics

This commit is contained in:
Sharan 2023-09-06 12:17:13 +02:00
parent 252128a863 a0335aa0c4
commit eacdf38482
23 changed files with 117 additions and 94 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
VERSION
View File

@ -1 +1 @@
3.0.0-alpha7-753-gaa28fe677
3.0.0-alpha7-784-g69239661c

24
python/damask/_crystal.py
View File

@ -2,11 +2,11 @@ from typing import Optional, Union, Dict, List, Tuple
import numpy as np
from ._typehints import FloatSequence, CrystalFamily, CrystalLattice, CrystalKinematics
from ._typehints import FloatSequence, CrystalFamily, BravaisLattice, CrystalKinematics
from . import util
from . import Rotation
lattice_symmetries: Dict[CrystalLattice, CrystalFamily] = {
lattice_symmetries: Dict[BravaisLattice, CrystalFamily] = {
'aP': 'triclinic',
'mP': 'monoclinic',
@ -27,7 +27,7 @@ lattice_symmetries: Dict[CrystalLattice, CrystalFamily] = {
'cF': 'cubic',
}
orientation_relationships: Dict[str, Dict[CrystalLattice,np.ndarray]] = {
orientation_relationships: Dict[str, Dict[BravaisLattice,np.ndarray]] = {
'KS': {
'cF': np.array([
[[-1, 0, 1],[ 1, 1, 1]],
@ -323,7 +323,7 @@ class Crystal():
def __init__(self, *,
family: Optional[CrystalFamily] = None,
lattice: Optional[CrystalLattice] = None,
lattice: Optional[BravaisLattice] = None,
a: Optional[float] = None, b: Optional[float] = None, c: Optional[float] = None,
alpha: Optional[float] = None, beta: Optional[float] = None, gamma: Optional[float] = None,
degrees: bool = False):
@ -548,7 +548,17 @@ class Crystal():
@property
def symmetry_operations(self) -> Rotation:
"""Symmetry operations as Rotations."""
"""
Return symmetry operations.
References
----------
U.F. Kocks et al.,
Texture and Anisotropy:
Preferred Orientations in Polycrystals and their Effect on Materials Properties.
Cambridge University Press 1998. Table II
"""
_symmetry_operations: Dict[CrystalFamily, List] = {
'cubic': [
[ 1.0, 0.0, 0.0, 0.0 ],
@ -772,7 +782,7 @@ class Crystal():
Directions and planes of deformation mode families.
"""
_kinematics: Dict[CrystalLattice, Dict[CrystalKinematics, List[np.ndarray]]] = {
_kinematics: Dict[BravaisLattice, Dict[CrystalKinematics, List[np.ndarray]]] = {
'cF': {
'slip': [np.array([
[+0,+1,-1, +1,+1,+1],
@ -1015,7 +1025,7 @@ class Crystal():
def relation_operations(self,
model: str) -> Tuple[CrystalLattice, Rotation]:
model: str) -> Tuple[BravaisLattice, Rotation]:
"""
Crystallographic orientation relationships for phase transformations.

4
python/damask/_orientation.py
View File

@ -3,7 +3,7 @@ from typing import Optional, Union, TypeVar
import numpy as np
from ._typehints import FloatSequence, IntSequence, CrystalFamily, CrystalLattice
from ._typehints import FloatSequence, IntSequence, CrystalFamily, BravaisLattice
from . import Rotation
from . import Crystal
from . import util
@ -73,7 +73,7 @@ class Orientation(Rotation,Crystal):
rotation: Union[FloatSequence, Rotation] = np.array([1.,0.,0.,0.]),
*,
family: Optional[CrystalFamily] = None,
lattice: Optional[CrystalLattice] = None,
lattice: Optional[BravaisLattice] = None,
a: Optional[float] = None, b: Optional[float] = None, c: Optional[float] = None,
alpha: Optional[float] = None, beta: Optional[float] = None, gamma: Optional[float] = None,
degrees: bool = False):

58
python/damask/_rotation.py
View File

@ -307,7 +307,8 @@ class Rotation:
p_m = self.quaternion[...,1:]
q_o = other.quaternion[...,0:1]
p_o = other.quaternion[...,1:]
q = (q_m*q_o - np.einsum('...i,...i',p_m,p_o).reshape(self.shape+(1,)))
qmo = q_m*q_o
q = (qmo - np.einsum('...i,...i',p_m,p_o).reshape(qmo.shape))
p = q_m*p_o + q_o*p_m + _P * np.cross(p_m,p_o)
return self.copy(Rotation(np.block([q,p]))._standardize())
else:
@ -1323,28 +1324,41 @@ class Rotation:
@staticmethod
def _qu2eu(qu: np.ndarray) -> np.ndarray:
"""Quaternion to Bunge Euler angles."""
q02 = qu[...,0:1]*qu[...,2:3]
q13 = qu[...,1:2]*qu[...,3:4]
q01 = qu[...,0:1]*qu[...,1:2]
q23 = qu[...,2:3]*qu[...,3:4]
q03_s = qu[...,0:1]**2+qu[...,3:4]**2
q12_s = qu[...,1:2]**2+qu[...,2:3]**2
chi = np.sqrt(q03_s*q12_s)
"""
Quaternion to Bunge Euler angles.
eu = np.where(np.abs(q12_s) < 1.e-8,
np.block([np.arctan2(-_P*2.*qu[...,0:1]*qu[...,3:4],qu[...,0:1]**2-qu[...,3:4]**2),
np.zeros(qu.shape[:-1]+(2,))]),
np.where(np.abs(q03_s) < 1.e-8,
np.block([np.arctan2( 2.*qu[...,1:2]*qu[...,2:3],qu[...,1:2]**2-qu[...,2:3]**2),
np.broadcast_to(np.pi,qu[...,0:1].shape),
np.zeros(qu.shape[:-1]+(1,))]),
np.block([np.arctan2((-_P*q02+q13)*chi, (-_P*q01-q23)*chi),
np.arctan2( 2.*chi, q03_s-q12_s ),
np.arctan2(( _P*q02+q13)*chi, (-_P*q01+q23)*chi)])
)
)
eu[np.abs(eu) < 1.e-6] = 0.
References
----------
E. Bernardes and S. Viollet, PLoS ONE 17(11):e0276302, 2022
https://doi.org/10.1371/journal.pone.0276302
"""
a = qu[...,0:1]
b = -_P*qu[...,3:4]
c = -_P*qu[...,1:2]
d = -_P*qu[...,2:3]
eu = np.block([
np.arctan2(b,a),
np.arccos(2*(a**2+b**2)/(a**2+b**2+c**2+d**2)-1),
np.arctan2(-d,c),
])
eu_sum = eu[...,0] + eu[...,2]
eu_diff = eu[...,0] - eu[...,2]
is_zero = np.isclose(eu[...,1],0.0)
is_pi = np.isclose(eu[...,1],np.pi)
is_ok = ~np.logical_or(is_zero,is_pi)
eu[...,0][is_zero] = 2*eu[...,0][is_zero]
eu[...,0][is_pi] = -2*eu[...,2][is_pi]
eu[...,2][~is_ok] = 0.0
eu[...,0][is_ok] = eu_diff[is_ok]
eu[...,2][is_ok] = eu_sum [is_ok]
eu[np.logical_or(np.abs(eu) < 1.e-6,
np.abs(eu-2*np.pi) < 1.e-6)] = 0.
return np.where(eu < 0., eu%(np.pi*np.array([2.,1.,2.])),eu)
@staticmethod

2
python/damask/_typehints.py
View File

@ -11,7 +11,7 @@ IntSequence = Union[np.ndarray,Sequence[int]]
StrSequence = Union[np.ndarray,Sequence[str]]
FileHandle = Union[TextIO, str, Path]
CrystalFamily = Union[None,Literal['triclinic', 'monoclinic', 'orthorhombic', 'tetragonal', 'hexagonal', 'cubic']]
CrystalLattice = Union[None,Literal['aP', 'mP', 'mS', 'oP', 'oS', 'oI', 'oF', 'tP', 'tI', 'hP', 'cP', 'cI', 'cF']]
BravaisLattice = Union[None,Literal['aP', 'mP', 'mS', 'oP', 'oS', 'oI', 'oF', 'tP', 'tI', 'hP', 'cP', 'cI', 'cF']]
CrystalKinematics = Literal['slip', 'twin']
NumpyRngSeed = Union[int, IntSequence, np.random.SeedSequence, np.random.Generator]
# BitGenerator does not exists in older numpy versions

7
python/tests/test_Rotation.py
View File

@ -975,6 +975,13 @@ class TestRotation:
assert np.allclose(rot_broadcast.quaternion[...,i,:], rot.quaternion)
@pytest.mark.parametrize('shape',[(3,2),(4,6)])
def test_broadcastcomposition(self,shape):
a = Rotation.from_random(shape[0])
b = Rotation.from_random(shape[1])
assert (a[:,np.newaxis]*b[np.newaxis,:]).allclose(a.broadcast_to(shape)*b.broadcast_to(shape))
@pytest.mark.parametrize('function,invalid',[(Rotation.from_quaternion, np.array([-1,0,0,0])),
(Rotation.from_quaternion, np.array([1,1,1,0])),
(Rotation.from_Euler_angles, np.array([1,4,0])),

24
src/crystal.f90
View File

@ -1902,7 +1902,7 @@ end function buildInteraction
!> @brief Build a local coordinate system on slip, twin, trans, cleavage systems
!> @details Order: Direction, plane (normal), and common perpendicular
!--------------------------------------------------------------------------------------------------
function buildCoordinateSystem(active,potential,system,lattice,cOverA)
function buildCoordinateSystem(active,potential,system,lattice,cOverA) result(coordinateSystem)
integer, dimension(:), intent(in) :: &
active, & !< # of active systems per family
@ -1914,7 +1914,7 @@ function buildCoordinateSystem(active,potential,system,lattice,cOverA)
real(pREAL), intent(in) :: &
cOverA
real(pREAL), dimension(3,3,sum(active)) :: &
buildCoordinateSystem
coordinateSystem
real(pREAL), dimension(3) :: &
direction, normal
@ -1937,10 +1937,14 @@ function buildCoordinateSystem(active,potential,system,lattice,cOverA)
select case(lattice)
case ('cF','cI','tI')
case ('cF','cI')
direction = system(1:3,p)
normal = system(4:6,p)
case ('tI')
direction = [ system(1,p), system(2,p), system(3,p)*cOverA ]
normal = [ system(4,p), system(5,p), system(6,p)/cOverA ]
case ('hP')
direction = [ system(1,p)*1.5_pREAL, &
(system(1,p)+2.0_pREAL*system(2,p))*sqrt(0.75_pREAL), &
@ -1954,10 +1958,10 @@ function buildCoordinateSystem(active,potential,system,lattice,cOverA)
end select
buildCoordinateSystem(1:3,1,a) = direction/norm2(direction)
buildCoordinateSystem(1:3,2,a) = normal /norm2(normal)
buildCoordinateSystem(1:3,3,a) = math_cross(direction/norm2(direction),&
normal /norm2(normal))
coordinateSystem(1:3,1,a) = direction/norm2(direction)
coordinateSystem(1:3,2,a) = normal /norm2(normal)
coordinateSystem(1:3,3,a) = math_cross(direction/norm2(direction),&
normal /norm2(normal))
end do activeSystems
end do activeFamilies
@ -2245,6 +2249,12 @@ subroutine crystal_selfTest
CoSy = buildCoordinateSystem([1],[1],system,'cF',0.0_pREAL)
if (any(dNeq(CoSy(1:3,1:3,1),math_I3))) error stop 'buildCoordinateSystem'
if (any(dNeq(buildCoordinateSystem(TI_NSLIPSYSTEM,TI_NSLIPSYSTEM,TI_SYSTEMSLIP,'cI',0.0_pReal), &
buildCoordinateSystem(TI_NSLIPSYSTEM,TI_NSLIPSYSTEM,TI_SYSTEMSLIP,'tI',1.0_pReal)))) &
error stop 'cI/tI coordinate system'
if (all(dEq( buildCoordinateSystem(TI_NSLIPSYSTEM,TI_NSLIPSYSTEM,TI_SYSTEMSLIP,'tI',1.1_pReal + r(1)*0.9_pReal), &
buildCoordinateSystem(TI_NSLIPSYSTEM,TI_NSLIPSYSTEM,TI_SYSTEMSLIP,'tI',1.0_pReal)))) &
error stop 'tI coordinate system'
do i = 1, 10
call random_number(C)
C_cF = crystal_symmetrize_C66(C,'cI')

2
src/phase.f90
View File

@ -427,12 +427,12 @@ subroutine phase_init
phase => phases%get_dict(ph)
refs = config_listReferences(phase,indent=3)
if (len(refs) > 0) print'(/,1x,a)', refs
phase_rho(ph) = phase%get_asReal('rho',defaultVal=0.0_pREAL)
phase_lattice(ph) = phase%get_asStr('lattice')
if (all(phase_lattice(ph) /= ['cF','cI','hP','tI'])) &
call IO_error(130,ext_msg='phase_init: '//phase%get_asStr('lattice'))
if (any(phase_lattice(ph) == ['hP','tI'])) &
phase_cOverA(ph) = phase%get_asReal('c/a')
phase_rho(ph) = phase%get_asReal('rho',defaultVal=0.0_pREAL)
allocate(phase_O_0(ph)%data(count(material_ID_phase==ph)))
end do

4
src/phase_damage_anisobrittle.f90
View File

@ -50,7 +50,7 @@ module function anisobrittle_init() result(mySources)
if (count(mySources) == 0) return
print'(/,1x,a)', '<<<+- phase:damage:anisobrittle init -+>>>'
print'(/,a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT)
print'(/,1x,a,1x,i0)', '# phases:',count(mySources); flush(IO_STDOUT)
phases => config_material%get_dict('phase')
@ -64,7 +64,7 @@ module function anisobrittle_init() result(mySources)
associate(prm => param(ph))
print'(/,1x,a,i0,a)', 'phase ',ph,': '//phases%key(ph)
print'(/,1x,a,1x,i0,a)', 'phase',ph,': '//phases%key(ph)
refs = config_listReferences(src,indent=3)
if (len(refs) > 0) print'(/,1x,a)', refs

4
src/phase_damage_isobrittle.f90
View File

@ -48,7 +48,7 @@ module function isobrittle_init() result(mySources)
if (count(mySources) == 0) return
print'(/,1x,a)', '<<<+- phase:damage:isobrittle init -+>>>'
print'(/,a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT)
print'(/,1x,a,1x,i0)', '# phases:',count(mySources); flush(IO_STDOUT)
phases => config_material%get_dict('phase')
@ -66,7 +66,7 @@ module function isobrittle_init() result(mySources)
prm%W_crit = src%get_asReal('G_crit')/src%get_asReal('l_c')
print'(/,1x,a,i0,a)', 'phase ',ph,': '//phases%key(ph)
print'(/,1x,a,1x,i0,a)', 'phase',ph,': '//phases%key(ph)
refs = config_listReferences(src,indent=3)
if (len(refs) > 0) print'(/,1x,a)', refs

28
src/phase_mechanical_eigen.f90
View File

@ -95,34 +95,6 @@ function kinematics_active(kinematics_label,kinematics_length) result(active_ki
end function kinematics_active
!--------------------------------------------------------------------------------------------------
!> @brief Checks if a damage kinematic mechanism is active.
!--------------------------------------------------------------------------------------------------
function kinematics_active2(kinematics_label) result(active_kinematics)
character(len=*), intent(in) :: kinematics_label !< name of kinematic mechanism
logical, dimension(:), allocatable :: active_kinematics
type(tDict), pointer :: &
phases, &
phase, &
kinematics_type
integer :: ph
phases => config_material%get_dict('phase')
allocate(active_kinematics(phases%length), source = .false.)
do ph = 1, phases%length
phase => phases%get_dict(ph)
kinematics_type => phase%get_dict('damage',defaultVal=emptyDict)
active_kinematics(ph) = kinematics_type%get_asStr('type',defaultVal='n/a') == kinematics_label
end do
end function kinematics_active2
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient
! ToDo: MD: S is Mi?

10
src/phase_mechanical_eigen_thermalexpansion.f90
View File

@ -27,7 +27,7 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics)
integer, intent(in) :: kinematics_length
logical, dimension(:,:), allocatable :: myKinematics
integer :: Ninstances, p, k
integer :: p, k
type(tList), pointer :: &
kinematics
type(tDict), pointer :: &
@ -37,15 +37,13 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics)
myKinematics = kinematics_active('thermalexpansion',kinematics_length)
Ninstances = count(myKinematics)
print'(/,a,i2)', ' # phases: ',Ninstances; flush(IO_STDOUT)
if (Ninstances == 0) return
if (count(myKinematics) == 0) return
print'(/,1x,a)', '<<<+- phase:mechanical:eigen:thermalexpansion init -+>>>'
print'(/,1x,a,1x,i0)', '# phases:',count(myKinematics); flush(IO_STDOUT)
phases => config_material%get_dict('phase')
allocate(param(Ninstances))
allocate(param(count(myKinematics)))
allocate(kinematics_thermal_expansion_instance(phases%length), source=0)
do p = 1, phases%length

4
src/phase_mechanical_elastic.f90
View File

@ -34,7 +34,7 @@ module subroutine elastic_init(phases)
print'(/,1x,a)', '<<<+- phase:mechanical:elastic init -+>>>'
print'(/,1x,a)', '<<<+- phase:mechanical:elastic:Hooke init -+>>>'
print'(/,a,i0)', ' # phases: ',phases%length; flush(IO_STDOUT)
print'(/,1x,a,1x,i0)', '# phases:',phases%length; flush(IO_STDOUT)
allocate(param(phases%length))
@ -43,7 +43,7 @@ module subroutine elastic_init(phases)
phase => phases%get_dict(ph)
mech => phase%get_dict('mechanical')
elastic => mech%get_dict('elastic')
print'(/,1x,a,i0,a)', 'phase ',ph,': '//phases%key(ph)
print'(/,1x,a,1x,i0,a)', 'phase',ph,': '//phases%key(ph)
refs = config_listReferences(elastic,indent=3)
if (len(refs) > 0) print'(/,1x,a)', refs
if (elastic%get_asStr('type') /= 'Hooke') call IO_error(200,ext_msg=elastic%get_asStr('type'))

2
src/phase_mechanical_plastic_dislotungsten.f90
View File

@ -108,11 +108,11 @@ module function plastic_dislotungsten_init() result(myPlasticity)
if (count(myPlasticity) == 0) return
print'(/,1x,a)', '<<<+- phase:mechanical:plastic:dislotungsten init -+>>>'
print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
print'(/,1x,a)', 'D. Cereceda et al., International Journal of Plasticity 78:242–256, 2016'
print'( 1x,a)', 'https://doi.org/10.1016/j.ijplas.2015.09.002'
print'(/,1x,a,1x,i0)', '# phases:',count(myPlasticity); flush(IO_STDOUT)
phases => config_material%get_dict('phase')
allocate(param(phases%length))

2
src/phase_mechanical_plastic_dislotwin.f90
View File

@ -155,7 +155,6 @@ module function plastic_dislotwin_init() result(myPlasticity)
if (count(myPlasticity) == 0) return
print'(/,1x,a)', '<<<+- phase:mechanical:plastic:dislotwin init -+>>>'
print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
print'(/,1x,a)', 'A. Ma and F. Roters, Acta Materialia 52(12):3603–3612, 2004'
print'( 1x,a)', 'https://doi.org/10.1016/j.actamat.2004.04.012'
@ -166,6 +165,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
print'(/,1x,a)', 'S.L. Wong et al., Acta Materialia 118:140–151, 2016'
print'( 1x,a)', 'https://doi.org/10.1016/j.actamat.2016.07.032'
print'(/,1x,a,1x,i0)', '# phases:',count(myPlasticity); flush(IO_STDOUT)
phases => config_material%get_dict('phase')
allocate(param(phases%length))

2
src/phase_mechanical_plastic_isotropic.f90
View File

@ -68,11 +68,11 @@ module function plastic_isotropic_init() result(myPlasticity)
if (count(myPlasticity) == 0) return
print'(/,1x,a)', '<<<+- phase:mechanical:plastic:isotropic init -+>>>'
print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
print'(/,1x,a)', 'T. Maiti and P. Eisenlohr, Scripta Materialia 145:37–40, 2018'
print'( 1x,a)', 'https://doi.org/10.1016/j.scriptamat.2017.09.047'
print'(/,1x,a,1x,i0)', '# phases:',count(myPlasticity); flush(IO_STDOUT)
phases => config_material%get_dict('phase')
allocate(param(phases%length))

2
src/phase_mechanical_plastic_kinehardening.f90
View File

@ -100,11 +100,11 @@ module function plastic_kinehardening_init() result(myPlasticity)
if (count(myPlasticity) == 0) return
print'(/,1x,a)', '<<<+- phase:mechanical:plastic:kinehardening init -+>>>'
print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
print'(/,1x,a)', 'J.A. Wollmershauser et al., International Journal of Fatigue 36:181–193, 2012'
print'( 1x,a)', 'https://doi.org/10.1016/j.ijfatigue.2011.07.008'
print'(/,1x,a,1x,i0)', '# phases:',count(myPlasticity); flush(IO_STDOUT)
phases => config_material%get_dict('phase')
allocate(param(phases%length))

2
src/phase_mechanical_plastic_none.f90
View File

@ -25,7 +25,7 @@ module function plastic_none_init() result(myPlasticity)
if (count(myPlasticity) == 0) return
print'(/,1x,a)', '<<<+- phase:mechanical:plastic:none init -+>>>'
print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
print'(/,1x,a,1x,i0)', '# phases:',count(myPlasticity); flush(IO_STDOUT)
phases => config_material%get_dict('phase')

2
src/phase_mechanical_plastic_nonlocal.f90
View File

@ -203,7 +203,6 @@ module function plastic_nonlocal_init() result(myPlasticity)
end if
print'(/,1x,a)', '<<<+- phase:mechanical:plastic:nonlocal init -+>>>'
print'(/,a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT)
print'(/,1x,a)', 'C. Reuber et al., Acta Materialia 71:333–348, 2014'
print'( 1x,a)', 'https://doi.org/10.1016/j.actamat.2014.03.012'
@ -211,6 +210,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
print'(/,1x,a)', 'C. Kords, Dissertation RWTH Aachen, 2014'
print'( 1x,a)', 'http://publications.rwth-aachen.de/record/229993'
print'(/,1x,a,1x,i0)', '# phases:',Ninstances; flush(IO_STDOUT)
phases => config_material%get_dict('phase')
allocate(geom(phases%length))

4
src/phase_mechanical_plastic_phenopowerlaw.f90
View File

@ -105,7 +105,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
if (count(myPlasticity) == 0) return
print'(/,1x,a)', '<<<+- phase:mechanical:plastic:phenopowerlaw init -+>>>'
print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
print'(/,1x,a,1x,i0)', '# phases:',count(myPlasticity); flush(IO_STDOUT)
phases => config_material%get_dict('phase')
@ -124,7 +124,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
mech => phase%get_dict('mechanical')
pl => mech%get_dict('plastic')
print'(/,1x,a,i0,a)', 'phase ',ph,': '//phases%key(ph)
print'(/,1x,a,1x,i0,a)', 'phase',ph,': '//phases%key(ph)
refs = config_listReferences(pl,indent=3)
if (len(refs) > 0) print'(/,1x,a)', refs

14
src/phase_thermal.f90
View File

@ -76,7 +76,9 @@ module subroutine thermal_init(phases)
thermal
type(tList), pointer :: &
sources
character(len=:), allocatable :: refs
character(len=:), allocatable :: &
refs, &
extmsg
integer :: &
ph, so, &
Nmembers
@ -88,6 +90,7 @@ module subroutine thermal_init(phases)
allocate(thermalState(phases%length))
allocate(thermal_Nsources(phases%length),source = 0)
allocate(param(phases%length))
extmsg = ''
do ph = 1, phases%length
Nmembers = count(material_ID_phase == ph)
@ -106,6 +109,15 @@ module subroutine thermal_init(phases)
if (any(phase_lattice(ph) == ['hP','tI'])) param(ph)%K(3,3) = thermal%get_asReal('K_33')
param(ph)%K = crystal_symmetrize_33(param(ph)%K,phase_lattice(ph))
! sanity checks
if ( param(ph)%C_p <= 0.0_pREAL ) extmsg = trim(extmsg)//' C_p'
if (any(param(ph)%K < 0.0_pREAL)) extmsg = trim(extmsg)//' K'
if ( phase_rho(ph) <= 0.0_pREAL ) extmsg = trim(extmsg)//' rho'
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg))
#if defined(__GFORTRAN__)
param(ph)%output = output_as1dStr(thermal)
#else

4
src/phase_thermal_source_dissipation.f90
View File

@ -42,7 +42,7 @@ module function source_dissipation_init(maxNsources) result(isMySource)
if (count(isMySource) == 0) return
print'(/,1x,a)', '<<<+- phase:thermal:source_dissipation init -+>>>'
print'(/,a,i2)', ' # phases: ',count(isMySource); flush(IO_STDOUT)
print'(/,1x,a,1x,i0)', '# phases:',count(isMySource); flush(IO_STDOUT)
phases => config_material%get_dict('phase')
@ -60,7 +60,7 @@ module function source_dissipation_init(maxNsources) result(isMySource)
if (isMySource(so,ph)) then
associate(prm => param(ph))
src => sources%get_dict(so)
print'(1x,a,i0,a,i0)', 'phase ',ph,' source ',so
print'(/,1x,a,1x,i0,1x,a,1x,a,1x,i0)', 'phase',ph,'('//phases%key(ph)//')','source',so
refs = config_listReferences(src,indent=3)
if (len(refs) > 0) print'(/,1x,a)', refs

4
src/phase_thermal_source_externalheat.f90
View File

@ -44,7 +44,7 @@ module function source_externalheat_init(maxNsources) result(isMySource)
if (count(isMySource) == 0) return
print'(/,1x,a)', '<<<+- phase:thermal:source_externalheat init -+>>>'
print'(/,a,i2)', ' # phases: ',count(isMySource); flush(IO_STDOUT)
print'(/,1x,a,1x,i0)', '# phases:',count(isMySource); flush(IO_STDOUT)
phases => config_material%get_dict('phase')
@ -64,7 +64,7 @@ module function source_externalheat_init(maxNsources) result(isMySource)
source_ID(ph) = so
associate(prm => param(ph))
src => sources%get_dict(so)
print'(1x,a,i0,a,i0)', 'phase ',ph,' source ',so
print'(/,1x,a,1x,i0,1x,a,1x,a,1x,i0)', 'phase',ph,'('//phases%key(ph)//')','source',so
refs = config_listReferences(src,indent=3)
if (len(refs) > 0) print'(/,1x,a)', refs