diff --git a/PRIVATE b/PRIVATE index f730bcb8d..ee4b1a03b 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit f730bcb8ddd7224011e70c57d0a5c03068532d2d +Subproject commit ee4b1a03b443a2c497a13c45c9498313442d731e diff --git a/install/MarcMentat/apply_DAMASK_modifications.py b/install/MarcMentat/apply_DAMASK_modifications.py index 43bdce26e..3eb6825c1 100755 --- a/install/MarcMentat/apply_DAMASK_modifications.py +++ b/install/MarcMentat/apply_DAMASK_modifications.py @@ -8,7 +8,7 @@ from pathlib import Path import damask -def copy_and_patch(patch,orig,msc_root,editor): +def copy_and_patch(patch,orig,marc_root,editor): try: shutil.copyfile(orig,orig.parent/patch.stem) except shutil.SameFileError: @@ -17,31 +17,31 @@ def copy_and_patch(patch,orig,msc_root,editor): with open(orig.parent/patch.stem) as f_in: content = f_in.read() with open(orig.parent/patch.stem,'w') as f_out: - f_out.write(content.replace('%INSTALLDIR%',msc_root).replace('%EDITOR%',editor)) + f_out.write(content.replace('%INSTALLDIR%',marc_root).replace('%EDITOR%',editor)) parser = argparse.ArgumentParser( - description='Apply DAMASK modification to MSC.Marc/Mentat', + description='Apply DAMASK modification to MSC Marc/Mentat', prog = Path(__file__).name, formatter_class=argparse.ArgumentDefaultsHelpFormatter) parser.add_argument('--editor', dest='editor', metavar='string', default='vi', - help='Name of the editor for MSC.Mentat (executable)') -parser.add_argument('--msc-root', dest='msc_root', metavar='string', - default=damask.solver._marc._msc_root, - help='MSC.Marc/Mentat root directory') -parser.add_argument('--msc-version', dest='msc_version', type=float, metavar='float', - default=damask.solver._marc._msc_version, - help='MSC.Marc/Mentat version') + help='Name of the editor for Marc Mentat (executable)') +parser.add_argument('--marc-root', dest='marc_root', metavar='string', + default=damask.solver._marc._marc_root, + help='Marc root directory') +parser.add_argument('--marc-version', dest='marc_version', type=float, metavar='float', + default=damask.solver._marc._marc_version, + help='Marc version') parser.add_argument('--damask-root', dest='damask_root', metavar = 'string', default=damask.solver._marc._damask_root, help='DAMASK root directory') args = parser.parse_args() -msc_root = Path(args.msc_root).expanduser() +marc_root = Path(args.marc_root).expanduser() damask_root = Path(args.damask_root).expanduser() -msc_version = int(args.msc_version) if str(args.msc_version).split('.')[1] == '0' else \ - args.msc_version +marc_version = int(args.marc_version) if str(args.marc_version).split('.')[1] == '0' else \ + args.marc_version matches = {'Marc_tools': [['comp_user','comp_damask_*mp'], ['run_marc','run_damask_*mp'], @@ -54,15 +54,15 @@ matches = {'Marc_tools': [['comp_user','comp_damask_*mp'], print('patching files...\n') -for directory in glob.glob(str(damask_root/f'install/MarcMentat/{msc_version}/*')): +for directory in glob.glob(str(damask_root/f'install/MarcMentat/{marc_version}/*')): for orig, mods in matches[Path(directory).name]: - product,subfolder = (msc_root/Path(directory)).name.split('_') - orig = msc_root/f'{product.lower()}{msc_version}/{subfolder}/{orig}' + product,subfolder = (marc_root/Path(directory)).name.split('_') + orig = marc_root/f'{product.lower()}{marc_version}/{subfolder}/{orig}' for patch in glob.glob(f'{directory}/{mods}.patch'): - copy_and_patch(Path(patch),orig,msc_root,args.editor) + copy_and_patch(Path(patch),orig,marc_root,args.editor) print('compiling Mentat menu binaries...') -executable = msc_root/f'mentat{msc_version}/bin/mentat' -menu_file = msc_root/f'mentat{msc_version}/menus/linux64/main.msb' +executable = marc_root/f'mentat{marc_version}/bin/mentat' +menu_file = marc_root/f'mentat{marc_version}/menus/linux64/main.msb' os.system(f'xvfb-run -a {executable} -compile {menu_file}') diff --git a/python/damask/VERSION b/python/damask/VERSION index f3e930998..72af4a57b 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-90-gbb6655045 +v3.0.0-alpha5-105-g020759996 diff --git a/python/damask/_configmaterial.py b/python/damask/_configmaterial.py index 375d3e2ce..d1091aa76 100644 --- a/python/damask/_configmaterial.py +++ b/python/damask/_configmaterial.py @@ -440,9 +440,10 @@ class ConfigMaterial(Config): """ N,n,shaped = 1,1,{} + map_dim = {'O':-1,'F_i':-2} for k,v in kwargs.items(): shaped[k] = np.array(v) - s = shaped[k].shape[:-1] if k=='O' else shaped[k].shape + s = shaped[k].shape[:map_dim.get(k,None)] N = max(N,s[0]) if len(s)>0 else N n = max(n,s[1]) if len(s)>1 else n @@ -451,11 +452,12 @@ class ConfigMaterial(Config): if 'v' not in kwargs: shaped['v'] = np.broadcast_to(1/n,(N,n)) + map_shape = {'O':(N,n,4),'F_i':(N,n,3,3)} for k,v in shaped.items(): - target = (N,n,4) if k=='O' else (N,n) + target = map_shape.get(k,(N,n)) obj = np.broadcast_to(v.reshape(util.shapeshifter(v.shape,target,mode='right')),target) for i in range(N): - if k in ['phase','O','v']: + if k in ['phase','O','v','F_i']: for j in range(n): mat[i]['constituents'][j][k] = obj[i,j].item() if isinstance(obj[i,j],np.generic) else obj[i,j] else: diff --git a/python/damask/solver/_marc.py b/python/damask/solver/_marc.py index c2de52dc6..1de16915d 100644 --- a/python/damask/solver/_marc.py +++ b/python/damask/solver/_marc.py @@ -3,14 +3,14 @@ import shlex import re from pathlib import Path -_msc_version = 2021.2 -_msc_root = '/opt/msc' +_marc_version = 2021.2 +_marc_root = '/opt/msc' _damask_root = str(Path(__file__).parents[3]) class Marc: - """Wrapper to run DAMASK with MSC.Marc.""" + """Wrapper to run DAMASK with MSC Marc.""" - def __init__(self,msc_version=_msc_version,msc_root=_msc_root,damask_root=_damask_root): + def __init__(self,marc_version=_marc_version,marc_root=_marc_root,damask_root=_damask_root): """ Create a Marc solver object. @@ -20,14 +20,14 @@ class Marc: Marc version """ - self.msc_version = msc_version - self.msc_root = Path(msc_root) + self.marc_version = marc_version + self.marc_root = Path(marc_root) self.damask_root = Path(damask_root) @property def library_path(self): - path_lib = self.msc_root/f'mentat{self.msc_version}/shlib/linux64' + path_lib = self.marc_root/f'mentat{self.marc_version}/shlib/linux64' if not path_lib.is_dir(): raise FileNotFoundError(f'library path "{path_lib}" not found') @@ -37,7 +37,7 @@ class Marc: @property def tools_path(self): - path_tools = self.msc_root/f'marc{self.msc_version}/tools' + path_tools = self.marc_root/f'marc{self.marc_version}/tools' if not path_tools.is_dir(): raise FileNotFoundError(f'tools path "{path_tools}" not found') diff --git a/python/tests/test_ConfigMaterial.py b/python/tests/test_ConfigMaterial.py index c968f429a..a713b636f 100644 --- a/python/tests/test_ConfigMaterial.py +++ b/python/tests/test_ConfigMaterial.py @@ -99,12 +99,15 @@ class TestConfigMaterial: @pytest.mark.parametrize('N,n,kw',[ (1,1,{'phase':'Gold', 'O':[1,0,0,0], + 'F_i':np.eye(3), 'homogenization':'SX'}), (3,1,{'phase':'Gold', 'O':Rotation.from_random(3), + 'F_i':np.broadcast_to(np.eye(3),(3,3,3)), 'homogenization':'SX'}), (2,3,{'phase':np.broadcast_to(['a','b','c'],(2,3)), 'O':Rotation.from_random((2,3)), + 'F_i':np.broadcast_to(np.eye(3),(2,3,3,3)), 'homogenization':['SX','PX']}), ]) def test_material_add(self,kw,N,n): diff --git a/src/material.f90 b/src/material.f90 index c83cd468c..52940ee4a 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -26,6 +26,7 @@ module material type(tRotationContainer), dimension(:), allocatable :: material_O_0 + type(tTensorContainer), dimension(:), allocatable :: material_F_i_0 integer, dimension(:), allocatable, public, protected :: & homogenization_Nconstituents !< number of grains in each homogenization @@ -48,6 +49,7 @@ module material public :: & tTensorContainer, & tRotationContainer, & + material_F_i_0, & material_O_0, & material_init @@ -159,16 +161,19 @@ subroutine parse() end do allocate(material_O_0(materials%length)) + allocate(material_F_i_0(materials%length)) do ma = 1, materials%length material => materials%get(ma) constituents => material%get('constituents') allocate(material_O_0(ma)%data(constituents%length)) + allocate(material_F_i_0(ma)%data(1:3,1:3,constituents%length)) do co = 1, constituents%length constituent => constituents%get(co) call material_O_0(ma)%data(co)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4)) - end do - end do + material_F_i_0(ma)%data(1:3,1:3,co) = constituent%get_as2dFloat('F_i',defaultVal=math_I3,requiredShape=[3,3]) + enddo + enddo end subroutine parse diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 4186b11cf..a74361198 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -205,6 +205,9 @@ module subroutine mechanical_init(phases) phases integer :: & + ce, & + co, & + ma, & ph, & en, & Nmembers @@ -261,15 +264,21 @@ module subroutine mechanical_init(phases) #endif enddo + do ce = 1, size(material_phaseID,2) + ma = discretization_materialAt((ce-1)/discretization_nIPs+1) + do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) + ph = material_phaseID(co,ce) + phase_mechanical_Fi0(ph)%data(1:3,1:3,material_phaseEntry(co,ce)) = material_F_i_0(ma)%data(1:3,1:3,co) + enddo + enddo + do ph = 1, phases%length do en = 1, count(material_phaseID == ph) phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = phase_O_0(ph)%data(en)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005) phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) & / math_det33(phase_mechanical_Fp0(ph)%data(1:3,1:3,en))**(1.0_pReal/3.0_pReal) - phase_mechanical_Fi0(ph)%data(1:3,1:3,en) = math_I3 - phase_mechanical_F0(ph)%data(1:3,1:3,en) = math_I3 - + phase_mechanical_F0(ph)%data(1:3,1:3,en) = math_I3 phase_mechanical_Fe(ph)%data(1:3,1:3,en) = math_inv33(matmul(phase_mechanical_Fi0(ph)%data(1:3,1:3,en), & phase_mechanical_Fp0(ph)%data(1:3,1:3,en))) ! assuming that euler angles are given in internal strain free configuration enddo diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index d0cc3a261..b2a9bdcc4 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -262,40 +262,43 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & i, j - Mp = matmul(matmul(transpose(Fi),Fi),S) + if (phase_plasticity(ph) == PLASTICITY_NONE_ID) then + Lp = 0.0_pReal + dLp_dFi = 0.0_pReal + dLp_dS = 0.0_pReal + else + Mp = matmul(matmul(transpose(Fi),Fi),S) - plasticType: select case (phase_plasticity(ph)) + plasticType: select case (phase_plasticity(ph)) - case (PLASTICITY_NONE_ID) plasticType - Lp = 0.0_pReal - dLp_dMp = 0.0_pReal + case (PLASTICITY_ISOTROPIC_ID) plasticType + call isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) - case (PLASTICITY_ISOTROPIC_ID) plasticType - call isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) + case (PLASTICITY_PHENOPOWERLAW_ID) plasticType + call phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) - case (PLASTICITY_PHENOPOWERLAW_ID) plasticType - call phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) + case (PLASTICITY_KINEHARDENING_ID) plasticType + call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) - case (PLASTICITY_KINEHARDENING_ID) plasticType - call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) + case (PLASTICITY_NONLOCAL_ID) plasticType + call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en) - case (PLASTICITY_NONLOCAL_ID) plasticType - call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en) + case (PLASTICITY_DISLOTWIN_ID) plasticType + call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en) - case (PLASTICITY_DISLOTWIN_ID) plasticType - call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en) + case (PLASTICITY_DISLOTUNGSTEN_ID) plasticType + call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en) - case (PLASTICITY_DISLOTUNGSTEN_ID) plasticType - call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en) + end select plasticType - end select plasticType + do i=1,3; do j=1,3 + dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + & + matmul(matmul(Fi,dLp_dMp(i,j,1:3,1:3)),S) + dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi) + end do; end do - do i=1,3; do j=1,3 - dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + & - matmul(matmul(Fi,dLp_dMp(i,j,1:3,1:3)),S) - dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi) - end do; end do + end if end subroutine plastic_LpAndItsTangents @@ -318,32 +321,34 @@ module function plastic_dotState(subdt,co,ip,el,ph,en) result(broken) logical :: broken - Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,en)),& - phase_mechanical_Fi(ph)%data(1:3,1:3,en)),phase_mechanical_S(ph)%data(1:3,1:3,en)) + if (phase_plasticity(ph) /= PLASTICITY_NONE_ID) then + Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,en)),& + phase_mechanical_Fi(ph)%data(1:3,1:3,en)),phase_mechanical_S(ph)%data(1:3,1:3,en)) - plasticType: select case (phase_plasticity(ph)) + plasticType: select case (phase_plasticity(ph)) - case (PLASTICITY_ISOTROPIC_ID) plasticType - call isotropic_dotState(Mp,ph,en) + case (PLASTICITY_ISOTROPIC_ID) plasticType + call isotropic_dotState(Mp,ph,en) - case (PLASTICITY_PHENOPOWERLAW_ID) plasticType - call phenopowerlaw_dotState(Mp,ph,en) + case (PLASTICITY_PHENOPOWERLAW_ID) plasticType + call phenopowerlaw_dotState(Mp,ph,en) - case (PLASTICITY_KINEHARDENING_ID) plasticType - call plastic_kinehardening_dotState(Mp,ph,en) + case (PLASTICITY_KINEHARDENING_ID) plasticType + call plastic_kinehardening_dotState(Mp,ph,en) - case (PLASTICITY_DISLOTWIN_ID) plasticType - call dislotwin_dotState(Mp,thermal_T(ph,en),ph,en) + case (PLASTICITY_DISLOTWIN_ID) plasticType + call dislotwin_dotState(Mp,thermal_T(ph,en),ph,en) - case (PLASTICITY_DISLOTUNGSTEN_ID) plasticType - call dislotungsten_dotState(Mp,thermal_T(ph,en),ph,en) + case (PLASTICITY_DISLOTUNGSTEN_ID) plasticType + call dislotungsten_dotState(Mp,thermal_T(ph,en),ph,en) + + case (PLASTICITY_NONLOCAL_ID) plasticType + call nonlocal_dotState(Mp,thermal_T(ph,en),subdt,ph,en,ip,el) + end select plasticType + end if - case (PLASTICITY_NONLOCAL_ID) plasticType - call nonlocal_dotState(Mp,thermal_T(ph,en),subdt,ph,en,ip,el) - end select plasticType broken = any(IEEE_is_NaN(plasticState(ph)%dotState(:,en))) - end function plastic_dotState @@ -390,8 +395,7 @@ module function plastic_deltaState(ph, en) result(broken) integer, intent(in) :: & ph, & en - logical :: & - broken + logical :: broken real(pReal), dimension(3,3) :: & Mp @@ -399,35 +403,34 @@ module function plastic_deltaState(ph, en) result(broken) myOffset, & mySize + broken = .false. - Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,en)),& - phase_mechanical_Fi(ph)%data(1:3,1:3,en)),phase_mechanical_S(ph)%data(1:3,1:3,en)) + select case (phase_plasticity(ph)) + case (PLASTICITY_NONLOCAL_ID,PLASTICITY_KINEHARDENING_ID) - plasticType: select case (phase_plasticity(ph)) + Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,en)),& + phase_mechanical_Fi(ph)%data(1:3,1:3,en)),& + phase_mechanical_S(ph)%data(1:3,1:3,en)) + + plasticType: select case (phase_plasticity(ph)) + + case (PLASTICITY_KINEHARDENING_ID) plasticType + call plastic_kinehardening_deltaState(Mp,ph,en) + + case (PLASTICITY_NONLOCAL_ID) plasticType + call plastic_nonlocal_deltaState(Mp,ph,en) + + end select plasticType - case (PLASTICITY_KINEHARDENING_ID) plasticType - call plastic_kinehardening_deltaState(Mp,ph,en) broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,en))) - - case (PLASTICITY_NONLOCAL_ID) plasticType - call plastic_nonlocal_deltaState(Mp,ph,en) - broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,en))) - - case default - broken = .false. - - end select plasticType - - if (.not. broken) then - select case(phase_plasticity(ph)) - case (PLASTICITY_NONLOCAL_ID,PLASTICITY_KINEHARDENING_ID) - + if (.not. broken) then myOffset = plasticState(ph)%offsetDeltaState mySize = plasticState(ph)%sizeDeltaState plasticState(ph)%state(myOffset + 1:myOffset + mySize,en) = & plasticState(ph)%state(myOffset + 1:myOffset + mySize,en) + plasticState(ph)%deltaState(1:mySize,en) - end select - end if + end if + + end select end function plastic_deltaState @@ -435,7 +438,7 @@ end function plastic_deltaState !-------------------------------------------------------------------------------------------------- !> @brief checks if a plastic module is active or not !-------------------------------------------------------------------------------------------------- -function plastic_active(plastic_label) result(active_plastic) +function plastic_active(plastic_label) result(active_plastic) character(len=*), intent(in) :: plastic_label !< type of plasticity model logical, dimension(:), allocatable :: active_plastic @@ -453,8 +456,8 @@ function plastic_active(plastic_label) result(active_plastic) phase => phases%get(ph) mech => phase%get('mechanical') pl => mech%get('plastic',defaultVal = emptyDict) - if (pl%get_asString('type',defaultVal='none') == plastic_label) active_plastic(ph) = .true. - end do + active_plastic(ph) = pl%get_asString('type',defaultVal='none') == plastic_label + enddo end function plastic_active